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