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