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 192 : void Dipole::registerKeywords(Keywords& keys) {
109 192 : Colvar::registerKeywords(keys);
110 192 : keys.setDisplayName("DIPOLE");
111 384 : keys.add("atoms","GROUP","the group of atoms we are calculating the dipole moment for");
112 384 : 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 384 : keys.addOutputComponent("x","COMPONENTS","the x-component of the dipole");
114 384 : keys.addOutputComponent("y","COMPONENTS","the y-component of the dipole");
115 384 : keys.addOutputComponent("z","COMPONENTS","the z-component of the dipole");
116 384 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
117 192 : keys.setValueDescription("the DIPOLE for these atoms");
118 192 : }
119 :
120 58 : Dipole::Dipole(const ActionOptions&ao):
121 : PLUMED_COLVAR_INIT(ao),
122 58 : components(false),
123 58 : value(1),
124 0 : derivs(1),
125 58 : virial(1) {
126 58 : parseAtomList(-1,ga_lista,this);
127 58 : charges.resize(ga_lista.size());
128 58 : components=(getModeAndSetupValues(this)==1);
129 58 : if( components ) {
130 4 : value.resize(3);
131 4 : derivs.resize(3);
132 4 : virial.resize(3);
133 4 : valuex=getPntrToComponent("x");
134 4 : valuey=getPntrToComponent("y");
135 4 : valuez=getPntrToComponent("z");
136 : }
137 124 : for(unsigned i=0; i<derivs.size(); ++i) {
138 66 : derivs[i].resize( ga_lista.size() );
139 : }
140 58 : parseFlag("NOPBC",nopbc);
141 58 : checkRead();
142 :
143 58 : if(nopbc) {
144 38 : log.printf(" without periodic boundary conditions\n");
145 : } else {
146 20 : log.printf(" using periodic boundary conditions\n");
147 : }
148 :
149 58 : requestAtoms(ga_lista);
150 58 : }
151 :
152 576 : void Dipole::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
153 1152 : aa->parseAtomList("GROUP",num,t);
154 576 : if( t.size()>0 ) {
155 573 : aa->log.printf(" of %u atoms\n",static_cast<unsigned>(t.size()));
156 2479 : for(unsigned int i=0; i<t.size(); ++i) {
157 1906 : aa->log.printf(" %d", t[i].serial());
158 : }
159 573 : aa->log.printf(" \n");
160 : }
161 576 : }
162 :
163 61 : unsigned Dipole::getModeAndSetupValues( ActionWithValue* av ) {
164 : bool c;
165 61 : av->parseFlag("COMPONENTS",c);
166 61 : if( c ) {
167 12 : av->addComponentWithDerivatives("x");
168 6 : av->componentIsNotPeriodic("x");
169 12 : av->addComponentWithDerivatives("y");
170 6 : av->componentIsNotPeriodic("y");
171 12 : av->addComponentWithDerivatives("z");
172 6 : av->componentIsNotPeriodic("z");
173 6 : return 1;
174 : }
175 55 : av->addValueWithDerivatives();
176 55 : av->setNotPeriodic();
177 : return 0;
178 : }
179 :
180 : // calculator
181 943 : void Dipole::calculate() {
182 943 : if(!nopbc) {
183 561 : makeWhole();
184 : }
185 : unsigned N=getNumberOfAtoms();
186 7398 : for(unsigned i=0; i<N; ++i) {
187 6455 : charges[i]=getCharge(i);
188 : }
189 :
190 943 : if(!components) {
191 788 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
192 6313 : for(unsigned i=0; i<N; i++) {
193 5525 : setAtomsDerivatives(i,derivs[0][i]);
194 : }
195 788 : setBoxDerivatives(virial[0]);
196 788 : setValue(value[0]);
197 : } else {
198 155 : calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this );
199 1085 : for(unsigned i=0; i<N; i++) {
200 930 : setAtomsDerivatives(valuex,i,derivs[0][i]);
201 930 : setAtomsDerivatives(valuey,i,derivs[1][i]);
202 930 : setAtomsDerivatives(valuez,i,derivs[2][i]);
203 : }
204 155 : setBoxDerivatives(valuex,virial[0]);
205 155 : setBoxDerivatives(valuey,virial[1]);
206 155 : setBoxDerivatives(valuez,virial[2]);
207 155 : valuex->set(value[0]);
208 155 : valuey->set(value[1]);
209 155 : valuez->set(value[2]);
210 : }
211 943 : }
212 :
213 4029 : void Dipole::calculateCV( const unsigned& mode, const std::vector<double>& masses, std::vector<double>& charges,
214 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
215 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
216 4029 : unsigned N=pos.size();
217 : double ctot=0.;
218 19802 : for(unsigned i=0; i<N; ++i) {
219 15773 : ctot += charges[i];
220 : }
221 4029 : ctot/=(double)N;
222 :
223 4029 : Vector dipje;
224 19802 : for(unsigned i=0; i<N; ++i) {
225 15773 : charges[i]-=ctot;
226 15773 : dipje += charges[i]*pos[i];
227 : }
228 :
229 4029 : if( mode==1 ) {
230 13419 : for(unsigned i=0; i<N; i++) {
231 10188 : derivs[0][i]=charges[i]*Vector(1.0,0.0,0.0);
232 10188 : derivs[1][i]=charges[i]*Vector(0.0,1.0,0.0);
233 10188 : derivs[2][i]=charges[i]*Vector(0.0,0.0,1.0);
234 : }
235 12924 : for(unsigned i=0; i<3; ++i ) {
236 9693 : vals[i] = dipje[i];
237 : }
238 : } else {
239 798 : vals[0] = dipje.modulo();
240 798 : double idip = 1./vals[0];
241 6383 : for(unsigned i=0; i<N; i++) {
242 5585 : double dfunc=charges[i]*idip;
243 5585 : derivs[0][i] = dfunc*dipje;
244 : }
245 : }
246 4029 : setBoxDerivativesNoPbc( pos, derivs, virial );
247 4029 : }
248 :
249 : }
250 : }
|
