Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-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 "core/ActionWithVirtualAtom.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/ActionSet.h"
26 : #include "MultiColvarBase.h"
27 : #include "CatomPack.h"
28 : #include "BridgedMultiColvarFunction.h"
29 : #include "vesselbase/StoreDataVessel.h"
30 :
31 : //+PLUMEDOC VATOM CENTER_OF_MULTICOLVAR
32 : /*
33 : Calculate a a weighted average position based on the value of some multicolvar.
34 :
35 : This action calculates the position of a new virtual atom using the following formula:
36 :
37 : \f[
38 : x_\alpha = \frac{1}{2\pi} \arctan \left[ \frac{ \sum_i w_i f_i \sin\left( 2\pi x_{i,\alpha} \right) }{ \sum_i w_i f_i \cos\left( 2\pi x_{i,\alpha} \right) } \right]
39 : \f]
40 :
41 : Where in this expression the \f$w_i\f$ values are a set of weights calculated within a multicolvar
42 : action and the \f$f_i\f$ are the values of the multicolvar functions. The \f$x_{i,\alpha}\f$ values are
43 : the positions (in scaled coordinates) associated with each of the multicolvars calculated.
44 :
45 : \bug The virial contribution for this type of virtual atom is not currently evaluated so do not use in bias functions unless the volume of the cell is fixed
46 :
47 : \par Examples
48 :
49 : Lets suppose that you are examining the formation of liquid droplets from gas. You may want to
50 : determine the center of mass of any of the droplets formed. In doing this calculation you recognize that
51 : the atoms in the liquid droplets will have a higher coordination number than those in the surrounding gas.
52 : As you want to calculate the position of the droplets you thus recognize that these atoms with high coordination
53 : numbers should have a high weight in the weighted average you are using to calculate the position of the droplet.
54 : You can thus calculate the position of the droplet using an input like the one shown below:
55 :
56 : \plumedfile
57 : c1: COORDINATIONNUMBER LOWMEM SPECIES=1-512 SWITCH={EXP D_0=4.0 R_0=0.5}
58 : cc: CENTER_OF_MULTICOLVAR DATA=c1
59 : \endplumedfile
60 :
61 : The first line here calculates the coordination numbers of all the atoms in the system. The virtual atom then uses the values
62 : of the coordination numbers calculated by the action labelled c1 when it calculates the Berry Phase average described above.
63 : (N.B. the \f$w_i\f$ in the above expression are all set equal to 1 in this case)
64 :
65 : The above input is fine we can, however, refine this somewhat by making use of a multicolvar transform action as shown below:
66 :
67 : \plumedfile
68 : c1: COORDINATIONNUMBER SPECIES=1-512 SWITCH={EXP D_0=4.0 R_0=0.5}
69 : cf: MTRANSFORM_MORE DATA=c1 SWITCH={RATIONAL D_0=2.0 R_0=0.1} LOWMEM
70 : cc: CENTER_OF_MULTICOLVAR DATA=cf
71 : \endplumedfile
72 :
73 : This input once again calculates the coordination numbers of all the atoms in the system. The middle line then transforms these
74 : coordination numbers to numbers between 0 and 1. Essentially any atom with a coordination number larger than 2.0 is given a weight
75 : of one and below this value the transformed value decays to zero. It is these transformed coordination numbers that are used to calculate
76 : the Berry phase average described in the previous section.
77 :
78 : */
79 : //+ENDPLUMEDOC
80 :
81 : namespace PLMD {
82 : namespace multicolvar {
83 :
84 :
85 : class CenterOfMultiColvar : public ActionWithVirtualAtom {
86 : private:
87 : unsigned comp;
88 : vesselbase::StoreDataVessel* mystash;
89 : MultiColvarBase* mycolv;
90 : public:
91 : static void registerKeywords( Keywords& keys );
92 : explicit CenterOfMultiColvar(const ActionOptions&ao);
93 : void calculate() override;
94 : };
95 :
96 13791 : PLUMED_REGISTER_ACTION(CenterOfMultiColvar,"CENTER_OF_MULTICOLVAR")
97 :
98 7 : void CenterOfMultiColvar::registerKeywords(Keywords& keys) {
99 7 : ActionWithVirtualAtom::registerKeywords(keys);
100 14 : keys.add("compulsory","DATA","find the average value for a multicolvar");
101 14 : keys.add("optional","COMPONENT","if your input multicolvar is a vector then specify which component you would like to use in calculating the weight");
102 7 : }
103 :
104 3 : CenterOfMultiColvar::CenterOfMultiColvar(const ActionOptions&ao):
105 : Action(ao),
106 3 : ActionWithVirtualAtom(ao) {
107 : std::string mlab;
108 3 : parse("DATA",mlab);
109 3 : mycolv= plumed.getActionSet().selectWithLabel<MultiColvarBase*>(mlab);
110 3 : if(!mycolv) {
111 0 : error("action labelled " + mlab + " does not exist or does not have vessels");
112 : }
113 : // Copy the atoms from the input multicolvar
114 3 : BridgedMultiColvarFunction* mybr=dynamic_cast<BridgedMultiColvarFunction*>( mycolv );
115 3 : if( mybr ) {
116 2 : requestAtoms( (mybr->getPntrToMultiColvar())->getAbsoluteIndexes() );
117 2 : comp=1;
118 : } else {
119 1 : if( mycolv->getNumberOfQuantities()>5 ) {
120 0 : int incomp=-1;
121 0 : parse("COMPONENT",incomp);
122 0 : if( incomp<0 ) {
123 0 : error("vector input but component was not specified");
124 : }
125 0 : comp=incomp;
126 : } else {
127 1 : comp=1;
128 : }
129 1 : requestAtoms( mycolv->getAbsoluteIndexes () );
130 : }
131 : // We need the derivatives
132 3 : mycolv->turnOnDerivatives();
133 3 : addDependency(mycolv);
134 3 : mystash = mycolv->buildDataStashes( NULL );
135 3 : log.printf(" building center of mass based on weights calculated in multicolvar action named %s \n",mycolv->getLabel().c_str() );
136 3 : }
137 :
138 4 : void CenterOfMultiColvar::calculate() {
139 : // Retrieve the periodic boundary conditions
140 4 : const Pbc& pbc=mycolv->getPbc();
141 4 : if( !pbc.isOrthorombic() ) {
142 0 : error("Berry phase does not work for non orthorhombic cells");
143 : }
144 :
145 : // Create a multivalue to store the derivatives
146 4 : MultiValue myvals( 7, mycolv->getNumberOfDerivatives() );
147 4 : myvals.clearAll();
148 4 : MultiValue tvals( mycolv->getNumberOfQuantities(), mycolv->getNumberOfDerivatives() );
149 4 : tvals.clearAll();
150 :
151 : // Now loop over all active multicolvars
152 4 : Vector stmp, ctmp, scom, ccom, sder, cder;
153 4 : scom.zero();
154 4 : ccom.zero();
155 : double norm=0;
156 4 : std::vector<double> cvals( mycolv->getNumberOfQuantities() );
157 1818 : for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
158 : // Retrieve value and derivatives
159 1814 : mystash->retrieveSequentialValue( i, false, cvals );
160 1814 : mystash->retrieveDerivatives( mycolv->getPositionInFullTaskList(i), false, tvals );
161 : // Convert position into fractionals
162 1814 : Vector fpos = pbc.realToScaled( mycolv->getCentralAtomPos( mycolv->getPositionInFullTaskList(i) ) );
163 : // Now accumulate Berry phase averages
164 7256 : for(unsigned j=0; j<3; ++j) {
165 5442 : stmp[j] = std::sin( 2*pi*fpos[j] );
166 5442 : ctmp[j] = std::cos( 2*pi*fpos[j] );
167 5442 : scom[j] += cvals[0]*cvals[comp]*stmp[j];
168 5442 : ccom[j] += cvals[0]*cvals[comp]*ctmp[j];
169 5442 : double icell = 1.0 / getPbc().getBox().getRow(j).modulo();
170 5442 : sder[j] = 2*pi*icell*cvals[0]*cvals[comp]*std::cos( 2*pi*fpos[j] );
171 5442 : cder[j]=-2*pi*icell*cvals[0]*cvals[comp]*std::sin( 2*pi*fpos[j] );
172 : }
173 : // Now accumulate derivatives
174 1623560 : for(unsigned k=0; k<tvals.getNumberActive(); ++k) {
175 1621746 : unsigned icomp=tvals.getActiveIndex(k);
176 1621746 : myvals.addDerivative( 0, icomp, cvals[0]*tvals.getDerivative( comp, icomp ) + cvals[comp]*tvals.getDerivative( 0, icomp ) );
177 6486984 : for(unsigned k=0; k<3; ++k) {
178 4865238 : myvals.addDerivative( 1+k, icomp, stmp[k]*( cvals[0]*tvals.getDerivative( comp, icomp ) +
179 4865238 : cvals[comp]*tvals.getDerivative( 0, icomp ) ) );
180 4865238 : myvals.addDerivative( 4+k, icomp, ctmp[k]*( cvals[0]*tvals.getDerivative( comp, icomp ) +
181 4865238 : cvals[comp]*tvals.getDerivative( 0, icomp ) ) );
182 : }
183 : }
184 : // Get the central atom pack
185 : CatomPack mypack;
186 1814 : mycolv->getCentralAtomPack( 0, mycolv->getPositionInFullTaskList(i), mypack );
187 3628 : for(unsigned j=0; j<mypack.getNumberOfAtomsWithDerivatives(); ++j) {
188 1814 : unsigned jder=3*mypack.getIndex(j);
189 : // Derivatives of sine
190 1814 : myvals.addDerivative( 1, jder+0, mypack.getDerivative(j, 0, sder) );
191 1814 : myvals.addDerivative( 2, jder+1, mypack.getDerivative(j, 1, sder) );
192 1814 : myvals.addDerivative( 3, jder+2, mypack.getDerivative(j, 2, sder) );
193 : // Derivatives of cosine
194 1814 : myvals.addDerivative( 4, jder+0, mypack.getDerivative(j, 0, cder) );
195 1814 : myvals.addDerivative( 5, jder+1, mypack.getDerivative(j, 1, cder) );
196 1814 : myvals.addDerivative( 6, jder+2, mypack.getDerivative(j, 2, cder) );
197 : }
198 1814 : norm += cvals[0]*cvals[comp];
199 1814 : tvals.clearAll();
200 : }
201 :
202 : // And now finish Berry phase average
203 4 : scom /= norm;
204 4 : ccom /=norm;
205 4 : Vector cpos;
206 16 : for(unsigned j=0; j<3; ++j) {
207 12 : cpos[j] = std::atan2( scom[j], ccom[j] ) / (2*pi);
208 : }
209 4 : Vector cart_pos = pbc.scaledToReal( cpos );
210 : setPosition(cart_pos);
211 : setMass(1.0); // This could be much cleverer but not without changing many things in PLMED
212 :
213 : // And derivatives
214 4 : Vector tander;
215 : myvals.updateDynamicList();
216 4 : double inv_weight = 1.0 / norm;
217 16 : for(unsigned j=0; j<3; ++j) {
218 12 : double tmp = scom[j] / ccom[j];
219 12 : tander[j] = getPbc().getBox().getRow(j).modulo() / (2*pi*( 1 + tmp*tmp ));
220 : }
221 5578 : for(unsigned i=0; i<myvals.getNumberActive(); ++i) {
222 5574 : unsigned ider=myvals.getActiveIndex(i);
223 22296 : for(unsigned j=0; j<3; ++j) {
224 16722 : double sderv = inv_weight*myvals.getDerivative(1+j,ider) - inv_weight*scom[j]*myvals.getDerivative(0,ider);
225 16722 : double cderv = inv_weight*myvals.getDerivative(4+j,ider) - inv_weight*ccom[j]*myvals.getDerivative(0,ider);
226 16722 : myvals.setDerivative( 1+j, ider, tander[j]*(sderv/ccom[j] - scom[j]*cderv/(ccom[j]*ccom[j])) );
227 : //if( j==2 ) std::printf("DERIV %d %10.4f %10.4f %10.4f %10.4f \n",i,myvals.getDerivative(0,ider),sderv,cderv,myvals.getDerivative(1+j,ider ) );
228 : }
229 : }
230 :
231 : // Atom derivatives
232 4 : std::vector<Tensor> fderiv( getNumberOfAtoms() );
233 2052 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
234 8192 : for(unsigned k=0; k<3; ++k) {
235 6144 : if( myvals.isActive(3*j+k) )
236 22152 : for(unsigned n=0; n<3; ++n) {
237 16614 : fderiv[j](k,n) = myvals.getDerivative( 1+n, 3*j+k );
238 : } else
239 2424 : for(unsigned n=0; n<3; ++n) {
240 1818 : fderiv[j](k,n) = 0;
241 : }
242 : }
243 : }
244 : setAtomsDerivatives( fderiv );
245 : // Box derivatives?
246 4 : }
247 :
248 : }
249 : }
|
