1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2012-2021 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 class Dipole : public Colvar {
67 vector<AtomNumber> ga_lista;
68 bool components;
69 bool nopbc;
70 public:
71 explicit Dipole(const ActionOptions&);
72 void calculate() override;
73 static void registerKeywords(Keywords& keys);
74 };
75
76 PLUMED_REGISTER_ACTION(Dipole,"DIPOLE")
77
registerKeywords(Keywords & keys)78 void Dipole::registerKeywords(Keywords& keys) {
79 Colvar::registerKeywords(keys);
80 keys.add("atoms","GROUP","the group of atoms we are calculating the dipole moment for");
81 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 keys.addOutputComponent("x","COMPONENTS","the x-component of the dipole");
83 keys.addOutputComponent("y","COMPONENTS","the y-component of the dipole");
84 keys.addOutputComponent("z","COMPONENTS","the z-component of the dipole");
85 }
86
Dipole(const ActionOptions & ao)87 Dipole::Dipole(const ActionOptions&ao):
88 PLUMED_COLVAR_INIT(ao),
89 components(false)
90 {
91 parseAtomList("GROUP",ga_lista);
92 parseFlag("COMPONENTS",components);
93 parseFlag("NOPBC",nopbc);
94 checkRead();
95 if(components) {
96 addComponentWithDerivatives("x"); componentIsNotPeriodic("x");
97 addComponentWithDerivatives("y"); componentIsNotPeriodic("y");
98 addComponentWithDerivatives("z"); componentIsNotPeriodic("z");
99 } else {
100 addValueWithDerivatives(); setNotPeriodic();
101 }
102
103 log.printf(" of %u atoms\n",static_cast<unsigned>(ga_lista.size()));
104 for(unsigned int i=0; i<ga_lista.size(); ++i) {
105 log.printf(" %d", ga_lista[i].serial());
106 }
107 log.printf(" \n");
108 if(nopbc) log.printf(" without periodic boundary conditions\n");
109 else log.printf(" using periodic boundary conditions\n");
110
111 requestAtoms(ga_lista);
112 }
113
114 // calculator
calculate()115 void Dipole::calculate()
116 {
117 if(!nopbc) makeWhole();
118 double ctot=0.;
119 unsigned N=getNumberOfAtoms();
120 vector<double> charges(N);
121 Vector dipje;
122
123 for(unsigned i=0; i<N; ++i) {
124 charges[i]=getCharge(i);
125 ctot+=charges[i];
126 }
127 ctot/=(double)N;
128
129 for(unsigned i=0; i<N; ++i) {
130 charges[i]-=ctot;
131 dipje += charges[i]*getPosition(i);
132 }
133
134 if(!components) {
135 double dipole = dipje.modulo();
136 double idip = 1./dipole;
137
138 for(unsigned i=0; i<N; i++) {
139 double dfunc=charges[i]*idip;
140 setAtomsDerivatives(i,dfunc*dipje);
141 }
142 setBoxDerivativesNoPbc();
143 setValue(dipole);
144 } else {
145 Value* valuex=getPntrToComponent("x");
146 Value* valuey=getPntrToComponent("y");
147 Value* valuez=getPntrToComponent("z");
148 for(unsigned i=0; i<N; i++) {
149 setAtomsDerivatives(valuex,i,charges[i]*Vector(1.0,0.0,0.0));
150 setAtomsDerivatives(valuey,i,charges[i]*Vector(0.0,1.0,0.0));
151 setAtomsDerivatives(valuez,i,charges[i]*Vector(0.0,0.0,1.0));
152 }
153 setBoxDerivativesNoPbc(valuex);
154 setBoxDerivativesNoPbc(valuey);
155 setBoxDerivativesNoPbc(valuez);
156 valuex->set(dipje[0]);
157 valuey->set(dipje[1]);
158 valuez->set(dipje[2]);
159 }
160 }
161
162 }
163 }
164