Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2023 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 "ColvarShortcut.h"
24 : #include "MultiColvarTemplate.h"
25 : #include "core/ActionRegister.h"
26 :
27 : namespace PLMD {
28 : namespace colvar {
29 :
30 : //+PLUMEDOC COLVAR DIPOLE
31 : /*
32 : Calculate the dipole moment for a group of atoms.
33 :
34 : The following input tells plumed to calculate the dipole for the group of atoms containing
35 : the atoms from 1-10 and print it every 5 steps
36 :
37 : ```plumed
38 : d: DIPOLE GROUP=1-10
39 : PRINT FILE=output STRIDE=5 ARG=d
40 : ```
41 :
42 : The output value from this input is a scalar that tells you the total magnitude of the dipole vector. If you would
43 : like to access the dipole vector directly you can use the command:
44 :
45 : ```plumed
46 : d: DIPOLE GROUP=1-10 COMPONENTS
47 : PRINT ARG=d.* FILE=output
48 : ```
49 :
50 : This command will output three values d.x, d.y and d.z, which are the x, y and z components of the dipole respectively.
51 :
52 : You can calculate three instinguishable dipoles using a single DIPOLE command by using an input like the one below:
53 :
54 : ```plumed
55 : d: DIPOLE GROUP1=1-10 GROUP2=11-20 GROUP3=21-30
56 : PRINT ARG=d FILE=output
57 : ```
58 :
59 : The output, d, here is a three dimensional vector. The first element of this vector is the magnitude of the dipole for
60 : atoms 1-10, the second is the magnitude of the dipole for atoms 11-20 and the third is the magnitude of the dipole for
61 : atoms 21-30. You can also obtain vector components for the three dipoles above by using the following input:
62 :
63 : ```plumed
64 : d: DIPOLE COMPONENTS GROUP1=1-10 GROUP2=11-20 GROUP3=21-30
65 : PRINT ARG=d.x,d.y,d.z FILE=output
66 : ```
67 :
68 : The output from the DIPOLE command now consists of three three dimensional vectors called d.x, d.y and d.z that contain the
69 : x, y and z components of the three dipoles respectively.
70 :
71 : When running with periodic boundary conditions, the atoms in every group should be
72 : in the proper periodic image. This is done automatically since PLUMED 2.5,
73 : by considering the ordered list of atoms and rebuilding the molecule with a procedure
74 : that is equivalent to that done in [WHOLEMOLECULES](WHOLEMOLECULES.md). Notice that
75 : rebuilding is local to this action. This is different from [WHOLEMOLECULES](WHOLEMOLECULES.md)
76 : which actually modifies the coordinates stored in PLUMED. If you want to recover the old behavior
77 : you should use the NOPBC flag. In that case you need to take care that atoms are in the correct
78 : periodic image.
79 :
80 : > [!CAUTION]
81 : > If the total charge Q of any of the specified groups is non zero, then a charge Q/N will be subtracted from every atom,
82 : > where N is the number of atoms. This implies that the dipole (which for a charged system depends
83 : > on the position) is computed on the geometric center of the group.
84 :
85 : */
86 : //+ENDPLUMEDOC
87 :
88 : //+PLUMEDOC COLVAR DIPOLE_SCALAR
89 : /*
90 : Calculate the dipole moment for a group of atoms.
91 :
92 : \par Examples
93 :
94 : */
95 : //+ENDPLUMEDOC
96 :
97 : //+PLUMEDOC MCOLVAR DIPOLE_VECTOR
98 : /*
99 : Calculate a vector of dipole moments for a set of groups of atoms.
100 :
101 : \par Examples
102 :
103 : */
104 : //+ENDPLUMEDOC
105 :
106 : class Dipole : public Colvar {
107 : std::vector<AtomNumber> ga_lista;
108 : bool components;
109 : bool nopbc;
110 : std::vector<double> value, masses, charges;
111 : std::vector<std::vector<Vector> > derivs;
112 : std::vector<Tensor> virial;
113 : Value* valuex=nullptr;
114 : Value* valuey=nullptr;
115 : Value* valuez=nullptr;
116 : public:
117 : explicit Dipole(const ActionOptions&);
118 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
119 : static unsigned getModeAndSetupValues( ActionWithValue* av );
120 : void calculate() override;
121 : static void registerKeywords(Keywords& keys);
122 : static void calculateCV( const unsigned& mode, const std::vector<double>& masses, std::vector<double>& charges,
123 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
124 : std::vector<Tensor>& virial, const ActionAtomistic* aa );
125 : };
126 :
127 : typedef ColvarShortcut<Dipole> DipoleShortcut;
128 : PLUMED_REGISTER_ACTION(DipoleShortcut,"DIPOLE")
129 : PLUMED_REGISTER_ACTION(Dipole,"DIPOLE_SCALAR")
130 : typedef MultiColvarTemplate<Dipole> DipoleMulti;
131 : PLUMED_REGISTER_ACTION(DipoleMulti,"DIPOLE_VECTOR")
132 :
133 192 : void Dipole::registerKeywords(Keywords& keys) {
134 192 : Colvar::registerKeywords(keys);
135 384 : keys.setDisplayName("DIPOLE");
136 192 : keys.add("atoms","GROUP","the group of atoms we are calculating the dipole moment for");
137 192 : 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");
138 384 : keys.addOutputComponent("x","COMPONENTS","scalar/vector","the x-component of the dipole");
139 384 : keys.addOutputComponent("y","COMPONENTS","scalar/vector","the y-component of the dipole");
140 384 : keys.addOutputComponent("z","COMPONENTS","scalar/vector","the z-component of the dipole");
141 192 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
142 384 : keys.setValueDescription("scalar/vector","the DIPOLE for these atoms");
143 192 : }
144 :
145 58 : Dipole::Dipole(const ActionOptions&ao):
146 : PLUMED_COLVAR_INIT(ao),
147 58 : components(false),
148 58 : value(1),
149 0 : derivs(1),
150 58 : virial(1) {
151 58 : parseAtomList(-1,ga_lista,this);
152 58 : charges.resize(ga_lista.size());
153 58 : components=(getModeAndSetupValues(this)==1);
154 58 : if( components ) {
155 4 : value.resize(3);
156 4 : derivs.resize(3);
157 4 : virial.resize(3);
158 4 : valuex=getPntrToComponent("x");
159 4 : valuey=getPntrToComponent("y");
160 4 : valuez=getPntrToComponent("z");
161 : }
162 124 : for(unsigned i=0; i<derivs.size(); ++i) {
163 66 : derivs[i].resize( ga_lista.size() );
164 : }
165 58 : parseFlag("NOPBC",nopbc);
166 58 : checkRead();
167 :
168 58 : if(nopbc) {
169 38 : log.printf(" without periodic boundary conditions\n");
170 : } else {
171 20 : log.printf(" using periodic boundary conditions\n");
172 : }
173 :
174 58 : requestAtoms(ga_lista);
175 58 : }
176 :
177 576 : void Dipole::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
178 1152 : aa->parseAtomList("GROUP",num,t);
179 576 : if( t.size()>0 ) {
180 573 : aa->log.printf(" of %u atoms\n",static_cast<unsigned>(t.size()));
181 2479 : for(unsigned int i=0; i<t.size(); ++i) {
182 1906 : aa->log.printf(" %d", t[i].serial());
183 : }
184 573 : aa->log.printf(" \n");
185 : }
186 576 : }
187 :
188 61 : unsigned Dipole::getModeAndSetupValues( ActionWithValue* av ) {
189 : bool c;
190 61 : av->parseFlag("COMPONENTS",c);
191 61 : if( c ) {
192 12 : av->addComponentWithDerivatives("x");
193 6 : av->componentIsNotPeriodic("x");
194 12 : av->addComponentWithDerivatives("y");
195 6 : av->componentIsNotPeriodic("y");
196 12 : av->addComponentWithDerivatives("z");
197 6 : av->componentIsNotPeriodic("z");
198 6 : return 1;
199 : }
200 55 : av->addValueWithDerivatives();
201 55 : av->setNotPeriodic();
202 : return 0;
203 : }
204 :
205 : // calculator
206 943 : void Dipole::calculate() {
207 943 : if(!nopbc) {
208 561 : makeWhole();
209 : }
210 : unsigned N=getNumberOfAtoms();
211 7398 : for(unsigned i=0; i<N; ++i) {
212 6455 : charges[i]=getCharge(i);
213 : }
214 :
215 943 : if(!components) {
216 788 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
217 6313 : for(unsigned i=0; i<N; i++) {
218 5525 : setAtomsDerivatives(i,derivs[0][i]);
219 : }
220 788 : setBoxDerivatives(virial[0]);
221 788 : setValue(value[0]);
222 : } else {
223 155 : calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this );
224 1085 : for(unsigned i=0; i<N; i++) {
225 930 : setAtomsDerivatives(valuex,i,derivs[0][i]);
226 930 : setAtomsDerivatives(valuey,i,derivs[1][i]);
227 930 : setAtomsDerivatives(valuez,i,derivs[2][i]);
228 : }
229 155 : setBoxDerivatives(valuex,virial[0]);
230 155 : setBoxDerivatives(valuey,virial[1]);
231 155 : setBoxDerivatives(valuez,virial[2]);
232 155 : valuex->set(value[0]);
233 155 : valuey->set(value[1]);
234 155 : valuez->set(value[2]);
235 : }
236 943 : }
237 :
238 4029 : void Dipole::calculateCV( const unsigned& mode, const std::vector<double>& masses, std::vector<double>& charges,
239 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
240 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
241 4029 : unsigned N=pos.size();
242 : double ctot=0.;
243 19802 : for(unsigned i=0; i<N; ++i) {
244 15773 : ctot += charges[i];
245 : }
246 4029 : ctot/=(double)N;
247 :
248 4029 : Vector dipje;
249 19802 : for(unsigned i=0; i<N; ++i) {
250 15773 : charges[i]-=ctot;
251 15773 : dipje += charges[i]*pos[i];
252 : }
253 :
254 4029 : if( mode==1 ) {
255 13419 : for(unsigned i=0; i<N; i++) {
256 10188 : derivs[0][i]=charges[i]*Vector(1.0,0.0,0.0);
257 10188 : derivs[1][i]=charges[i]*Vector(0.0,1.0,0.0);
258 10188 : derivs[2][i]=charges[i]*Vector(0.0,0.0,1.0);
259 : }
260 12924 : for(unsigned i=0; i<3; ++i ) {
261 9693 : vals[i] = dipje[i];
262 : }
263 : } else {
264 798 : vals[0] = dipje.modulo();
265 798 : double idip = 1./vals[0];
266 6383 : for(unsigned i=0; i<N; i++) {
267 5585 : double dfunc=charges[i]*idip;
268 5585 : derivs[0][i] = dfunc*dipje;
269 : }
270 : }
271 4029 : setBoxDerivativesNoPbc( pos, derivs, virial );
272 4029 : }
273 :
274 : }
275 : }
|
