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 10425 : PLUMED_REGISTER_ACTION(CenterOfMultiColvar,"CENTER_OF_MULTICOLVAR")
97 :
98 4 : void CenterOfMultiColvar::registerKeywords(Keywords& keys) {
99 4 : ActionWithVirtualAtom::registerKeywords(keys);
100 8 : keys.add("compulsory","DATA","find the average value for a multicolvar");
101 8 : 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 4 : }
103 :
104 3 : CenterOfMultiColvar::CenterOfMultiColvar(const ActionOptions&ao):
105 : Action(ao),
106 3 : ActionWithVirtualAtom(ao)
107 : {
108 3 : std::string mlab; parse("DATA",mlab);
109 3 : mycolv= plumed.getActionSet().selectWithLabel<MultiColvarBase*>(mlab);
110 3 : if(!mycolv) error("action labelled " + mlab + " does not exist or does not have vessels");
111 : // Copy the atoms from the input multicolvar
112 3 : BridgedMultiColvarFunction* mybr=dynamic_cast<BridgedMultiColvarFunction*>( mycolv );
113 3 : if( mybr ) {
114 2 : requestAtoms( (mybr->getPntrToMultiColvar())->getAbsoluteIndexes() ); comp=1;
115 : } else {
116 1 : if( mycolv->getNumberOfQuantities()>5 ) {
117 0 : int incomp=-1; parse("COMPONENT",incomp);
118 0 : if( incomp<0 ) error("vector input but component was not specified");
119 0 : comp=incomp;
120 : } else {
121 1 : comp=1;
122 : }
123 1 : requestAtoms( mycolv->getAbsoluteIndexes () );
124 : }
125 : // We need the derivatives
126 3 : mycolv->turnOnDerivatives(); addDependency(mycolv);
127 3 : mystash = mycolv->buildDataStashes( NULL );
128 3 : log.printf(" building center of mass based on weights calculated in multicolvar action named %s \n",mycolv->getLabel().c_str() );
129 3 : }
130 :
131 4 : void CenterOfMultiColvar::calculate() {
132 : // Retrieve the periodic boundary conditions
133 4 : const Pbc& pbc=mycolv->getPbc();
134 4 : if( !pbc.isOrthorombic() ) error("Berry phase does not work for non orthorhombic cells");
135 :
136 : // Create a multivalue to store the derivatives
137 4 : MultiValue myvals( 7, mycolv->getNumberOfDerivatives() ); myvals.clearAll();
138 4 : MultiValue tvals( mycolv->getNumberOfQuantities(), mycolv->getNumberOfDerivatives() );
139 4 : tvals.clearAll();
140 :
141 : // Now loop over all active multicolvars
142 4 : Vector stmp, ctmp, scom, ccom, sder, cder;
143 4 : scom.zero(); ccom.zero(); double norm=0;
144 4 : std::vector<double> cvals( mycolv->getNumberOfQuantities() );
145 1818 : for(unsigned i=0; i<mystash->getNumberOfStoredValues(); ++i) {
146 : // Retrieve value and derivatives
147 1814 : mystash->retrieveSequentialValue( i, false, cvals );
148 1814 : mystash->retrieveDerivatives( mycolv->getPositionInFullTaskList(i), false, tvals );
149 : // Convert position into fractionals
150 1814 : Vector fpos = pbc.realToScaled( mycolv->getCentralAtomPos( mycolv->getPositionInFullTaskList(i) ) );
151 : // Now accumulate Berry phase averages
152 7256 : for(unsigned j=0; j<3; ++j) {
153 5442 : stmp[j] = std::sin( 2*pi*fpos[j] ); ctmp[j] = std::cos( 2*pi*fpos[j] );
154 5442 : scom[j] += cvals[0]*cvals[comp]*stmp[j]; ccom[j] += cvals[0]*cvals[comp]*ctmp[j];
155 5442 : double icell = 1.0 / getPbc().getBox().getRow(j).modulo();
156 5442 : sder[j] = 2*pi*icell*cvals[0]*cvals[comp]*std::cos( 2*pi*fpos[j] );
157 5442 : cder[j]=-2*pi*icell*cvals[0]*cvals[comp]*std::sin( 2*pi*fpos[j] );
158 : }
159 : // Now accumulate derivatives
160 1623560 : for(unsigned k=0; k<tvals.getNumberActive(); ++k) {
161 1621746 : unsigned icomp=tvals.getActiveIndex(k);
162 1621746 : myvals.addDerivative( 0, icomp, cvals[0]*tvals.getDerivative( comp, icomp ) + cvals[comp]*tvals.getDerivative( 0, icomp ) );
163 6486984 : for(unsigned k=0; k<3; ++k) {
164 4865238 : myvals.addDerivative( 1+k, icomp, stmp[k]*( cvals[0]*tvals.getDerivative( comp, icomp ) +
165 4865238 : cvals[comp]*tvals.getDerivative( 0, icomp ) ) );
166 4865238 : myvals.addDerivative( 4+k, icomp, ctmp[k]*( cvals[0]*tvals.getDerivative( comp, icomp ) +
167 4865238 : cvals[comp]*tvals.getDerivative( 0, icomp ) ) );
168 : }
169 : }
170 : // Get the central atom pack
171 1814 : CatomPack mypack; mycolv->getCentralAtomPack( 0, mycolv->getPositionInFullTaskList(i), mypack );
172 3628 : for(unsigned j=0; j<mypack.getNumberOfAtomsWithDerivatives(); ++j) {
173 1814 : unsigned jder=3*mypack.getIndex(j);
174 : // Derivatives of sine
175 1814 : myvals.addDerivative( 1, jder+0, mypack.getDerivative(j, 0, sder) );
176 1814 : myvals.addDerivative( 2, jder+1, mypack.getDerivative(j, 1, sder) );
177 1814 : myvals.addDerivative( 3, jder+2, mypack.getDerivative(j, 2, sder) );
178 : // Derivatives of cosine
179 1814 : myvals.addDerivative( 4, jder+0, mypack.getDerivative(j, 0, cder) );
180 1814 : myvals.addDerivative( 5, jder+1, mypack.getDerivative(j, 1, cder) );
181 1814 : myvals.addDerivative( 6, jder+2, mypack.getDerivative(j, 2, cder) );
182 : }
183 1814 : norm += cvals[0]*cvals[comp]; tvals.clearAll();
184 1814 : }
185 :
186 : // And now finish Berry phase average
187 4 : scom /= norm; ccom /=norm; Vector cpos;
188 16 : for(unsigned j=0; j<3; ++j) cpos[j] = std::atan2( scom[j], ccom[j] ) / (2*pi);
189 4 : Vector cart_pos = pbc.scaledToReal( cpos );
190 : setPosition(cart_pos); setMass(1.0); // This could be much cleverer but not without changing many things in PLMED
191 :
192 : // And derivatives
193 4 : Vector tander; myvals.updateDynamicList(); double inv_weight = 1.0 / norm;
194 16 : for(unsigned j=0; j<3; ++j) {
195 12 : double tmp = scom[j] / ccom[j];
196 12 : tander[j] = getPbc().getBox().getRow(j).modulo() / (2*pi*( 1 + tmp*tmp ));
197 : }
198 5578 : for(unsigned i=0; i<myvals.getNumberActive(); ++i) {
199 5574 : unsigned ider=myvals.getActiveIndex(i);
200 22296 : for(unsigned j=0; j<3; ++j) {
201 16722 : double sderv = inv_weight*myvals.getDerivative(1+j,ider) - inv_weight*scom[j]*myvals.getDerivative(0,ider);
202 16722 : double cderv = inv_weight*myvals.getDerivative(4+j,ider) - inv_weight*ccom[j]*myvals.getDerivative(0,ider);
203 16722 : myvals.setDerivative( 1+j, ider, tander[j]*(sderv/ccom[j] - scom[j]*cderv/(ccom[j]*ccom[j])) );
204 : //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 ) );
205 : }
206 : }
207 :
208 : // Atom derivatives
209 4 : std::vector<Tensor> fderiv( getNumberOfAtoms() );
210 2052 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
211 8192 : for(unsigned k=0; k<3; ++k) {
212 22758 : if( myvals.isActive(3*j+k) ) for(unsigned n=0; n<3; ++n) fderiv[j](k,n) = myvals.getDerivative( 1+n, 3*j+k );
213 2424 : else for(unsigned n=0; n<3; ++n) fderiv[j](k,n) = 0;
214 : }
215 : }
216 : setAtomsDerivatives( fderiv );
217 : // Box derivatives?
218 4 : }
219 :
220 : }
221 : }
|
