Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2012-2020 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "Colvar.h"
      23             : #include "ActionRegister.h"
      24             : 
      25             : #include <string>
      26             : #include <cmath>
      27             : 
      28             : using namespace std;
      29             : 
      30             : namespace PLMD {
      31             : namespace colvar {
      32             : 
      33             : //+PLUMEDOC COLVAR DIPOLE
      34             : /*
      35             : Calculate the dipole moment for a group of atoms.
      36             : 
      37             : When running with periodic boundary conditions, the atoms should be
      38             : in the proper periodic image. This is done automatically since PLUMED 2.5,
      39             : by considering the ordered list of atoms and rebuilding the molecule with a procedure
      40             : that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
      41             : rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
      42             : which actually modifies the coordinates stored in PLUMED.
      43             : 
      44             : In case you want to recover the old behavior you should use the NOPBC flag.
      45             : In that case you need to take care that atoms are in the correct
      46             : periodic image.
      47             : 
      48             : \par Examples
      49             : 
      50             : The following tells plumed to calculate the dipole of the group of atoms containing
      51             : the atoms from 1-10 and print it every 5 steps
      52             : \plumedfile
      53             : d: DIPOLE GROUP=1-10
      54             : PRINT FILE=output STRIDE=5 ARG=d
      55             : \endplumedfile
      56             : 
      57             : \attention
      58             : If the total charge Q of the group in non zero, then a charge Q/N will be subtracted to every atom,
      59             : where N is the number of atoms. This implies that the dipole (which for a charged system depends
      60             : on the position) is computed on the geometric center of the group.
      61             : 
      62             : 
      63             : */
      64             : //+ENDPLUMEDOC
      65             : 
      66         159 : class Dipole : public Colvar {
      67             :   vector<AtomNumber> ga_lista;
      68             :   bool components;
      69             :   bool nopbc;
      70             : public:
      71             :   explicit Dipole(const ActionOptions&);
      72             :   virtual void calculate();
      73             :   static void registerKeywords(Keywords& keys);
      74             : };
      75             : 
      76        7462 : PLUMED_REGISTER_ACTION(Dipole,"DIPOLE")
      77             : 
      78          54 : void Dipole::registerKeywords(Keywords& keys) {
      79          54 :   Colvar::registerKeywords(keys);
      80         216 :   keys.add("atoms","GROUP","the group of atoms we are calculating the dipole moment for");
      81         162 :   keys.addFlag("COMPONENTS",false,"calculate the x, y and z components of the dipole separately and store them as label.x, label.y and label.z");
      82         216 :   keys.addOutputComponent("x","COMPONENTS","the x-component of the dipole");
      83         216 :   keys.addOutputComponent("y","COMPONENTS","the y-component of the dipole");
      84         216 :   keys.addOutputComponent("z","COMPONENTS","the z-component of the dipole");
      85          54 : }
      86             : 
      87          53 : Dipole::Dipole(const ActionOptions&ao):
      88             :   PLUMED_COLVAR_INIT(ao),
      89          53 :   components(false)
      90             : {
      91         106 :   parseAtomList("GROUP",ga_lista);
      92         106 :   parseFlag("COMPONENTS",components);
      93         106 :   parseFlag("NOPBC",nopbc);
      94          53 :   checkRead();
      95          53 :   if(components) {
      96           6 :     addComponentWithDerivatives("x"); componentIsNotPeriodic("x");
      97           6 :     addComponentWithDerivatives("y"); componentIsNotPeriodic("y");
      98           6 :     addComponentWithDerivatives("z"); componentIsNotPeriodic("z");
      99             :   } else {
     100          51 :     addValueWithDerivatives(); setNotPeriodic();
     101             :   }
     102             : 
     103         106 :   log.printf("  of %u atoms\n",static_cast<unsigned>(ga_lista.size()));
     104        1075 :   for(unsigned int i=0; i<ga_lista.size(); ++i) {
     105         646 :     log.printf("  %d", ga_lista[i].serial());
     106             :   }
     107          53 :   log.printf("  \n");
     108          53 :   if(nopbc) log.printf("  without periodic boundary conditions\n");
     109          15 :   else      log.printf("  using periodic boundary conditions\n");
     110             : 
     111          53 :   requestAtoms(ga_lista);
     112          53 : }
     113             : 
     114             : // calculator
     115         918 : void Dipole::calculate()
     116             : {
     117         918 :   if(!nopbc) makeWhole();
     118             :   double ctot=0.;
     119             :   unsigned N=getNumberOfAtoms();
     120         918 :   vector<double> charges(N);
     121         918 :   Vector dipje;
     122             : 
     123       13568 :   for(unsigned i=0; i<N; ++i) {
     124       12650 :     charges[i]=getCharge(i);
     125        6325 :     ctot+=charges[i];
     126             :   }
     127         918 :   ctot/=(double)N;
     128             : 
     129       13568 :   for(unsigned i=0; i<N; ++i) {
     130       12650 :     charges[i]-=ctot;
     131       12650 :     dipje += charges[i]*getPosition(i);
     132             :   }
     133             : 
     134         918 :   if(!components) {
     135         773 :     double dipole = dipje.modulo();
     136         773 :     double idip = 1./dipole;
     137             : 
     138       11683 :     for(unsigned i=0; i<N; i++) {
     139       10910 :       double dfunc=charges[i]*idip;
     140       10910 :       setAtomsDerivatives(i,dfunc*dipje);
     141             :     }
     142         773 :     setBoxDerivativesNoPbc();
     143         773 :     setValue(dipole);
     144             :   } else {
     145         290 :     Value* valuex=getPntrToComponent("x");
     146         290 :     Value* valuey=getPntrToComponent("y");
     147         290 :     Value* valuez=getPntrToComponent("z");
     148        1885 :     for(unsigned i=0; i<N; i++) {
     149        1740 :       setAtomsDerivatives(valuex,i,charges[i]*Vector(1.0,0.0,0.0));
     150        1740 :       setAtomsDerivatives(valuey,i,charges[i]*Vector(0.0,1.0,0.0));
     151        1740 :       setAtomsDerivatives(valuez,i,charges[i]*Vector(0.0,0.0,1.0));
     152             :     }
     153         145 :     setBoxDerivativesNoPbc(valuex);
     154         145 :     setBoxDerivativesNoPbc(valuey);
     155         145 :     setBoxDerivativesNoPbc(valuez);
     156         145 :     valuex->set(dipje[0]);
     157         145 :     valuey->set(dipje[1]);
     158         145 :     valuez->set(dipje[2]);
     159             :   }
     160         918 : }
     161             : 
     162             : }
     163        5517 : }