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