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 "TrigonometricPathVessel.h"
23 : #include "vesselbase/VesselRegister.h"
24 : #include "tools/PDB.h"
25 :
26 : namespace PLMD {
27 : namespace mapping {
28 :
29 13788 : PLUMED_REGISTER_VESSEL(TrigonometricPathVessel,"GPATH")
30 :
31 3 : void TrigonometricPathVessel::registerKeywords( Keywords& keys ) {
32 3 : StoreDataVessel::registerKeywords(keys);
33 3 : }
34 :
35 4595 : void TrigonometricPathVessel::reserveKeyword( Keywords& keys ) {
36 9190 : keys.reserve("vessel","GPATH","calculate the position on the path using trigonometry");
37 9190 : keys.addOutputComponent("gspath","GPATH","the position on the path calculated using trigonometry");
38 9190 : keys.addOutputComponent("gzpath","GPATH","the distance from the path calculated using trigonometry");
39 4595 : }
40 :
41 3 : TrigonometricPathVessel::TrigonometricPathVessel( const vesselbase::VesselOptions& da ):
42 : StoreDataVessel(da),
43 6 : projdir(ReferenceConfigurationOptions("DIRECTION")),
44 3 : mydpack1( 1, getAction()->getNumberOfDerivatives() ),
45 3 : mydpack2( 1, getAction()->getNumberOfDerivatives() ),
46 3 : mydpack3( 1, getAction()->getNumberOfDerivatives() ),
47 3 : mypack1( 0, 0, mydpack1 ),
48 3 : mypack2( 0, 0, mydpack2 ),
49 9 : mypack3( 0, 0, mydpack3 ) {
50 3 : mymap=dynamic_cast<Mapping*>( getAction() );
51 3 : plumed_massert( mymap, "Trigonometric path vessel can only be used with mappings");
52 : // Retrieve the index of the property in the underlying mapping
53 3 : if( mymap->getNumberOfProperties()!=1 ) {
54 0 : error("cannot use trigonometric paths when there are multiple properties");
55 : }
56 :
57 125 : for(unsigned i=0; i<mymap->getFullNumberOfTasks(); ++i) {
58 122 : if( mymap->getTaskCode(i)!=mymap->getPositionInFullTaskList(i) ) {
59 0 : error("mismatched tasks and codes");
60 : }
61 : }
62 3 : mymap->addComponentWithDerivatives("gspath");
63 3 : mymap->componentIsNotPeriodic("gspath");
64 3 : sp=mymap->copyOutput( mymap->getNumberOfComponents()-1 );
65 3 : sp->resizeDerivatives( mymap->getNumberOfDerivatives() );
66 3 : mymap->addComponentWithDerivatives("gzpath");
67 3 : mymap->componentIsNotPeriodic("gzpath");
68 3 : zp=mymap->copyOutput( mymap->getNumberOfComponents()-1 );
69 3 : zp->resizeDerivatives( mymap->getNumberOfDerivatives() );
70 :
71 : // Check we have PCA
72 3 : ReferenceConfiguration* ref0=mymap->getReferenceConfiguration(0);
73 125 : for(unsigned i=0; i<mymap->getFullNumberOfTasks(); ++i) {
74 122 : if( !(mymap->getReferenceConfiguration(i))->pcaIsEnabledForThisReference() ) {
75 0 : error("pca must be implemented in order to use trigometric path");
76 : }
77 244 : if( ref0->getName()!=(mymap->getReferenceConfiguration(i))->getName() ) {
78 0 : error("cannot use mixed metrics");
79 : }
80 122 : if( mymap->getNumberOfAtoms()!=(mymap->getReferenceConfiguration(i))->getNumberOfReferencePositions() ) {
81 0 : error("all frames must use the same set of atoms");
82 : }
83 122 : if( mymap->getNumberOfArguments()!=(mymap->getReferenceConfiguration(i))->getNumberOfReferenceArguments() ) {
84 0 : error("all frames must use the same set of arguments");
85 : }
86 : }
87 :
88 3 : cargs.resize( mymap->getNumberOfArguments() );
89 3 : std::vector<std::string> argument_names( mymap->getNumberOfArguments() );
90 7 : for(unsigned i=0; i<mymap->getNumberOfArguments(); ++i) {
91 8 : argument_names[i] = (mymap->getPntrToArgument(i))->getName();
92 : }
93 3 : PDB mypdb;
94 3 : mypdb.setAtomNumbers( mymap->getAbsoluteIndexes() );
95 3 : mypdb.addBlockEnd( mymap->getAbsoluteIndexes().size() );
96 3 : if( argument_names.size()>0 ) {
97 2 : mypdb.setArgumentNames( argument_names );
98 : }
99 3 : projdir.read( mypdb );
100 3 : mypack1.resize( mymap->getNumberOfArguments(), mymap->getNumberOfAtoms() );
101 3 : ref0->setupPCAStorage( mypack1 );
102 3 : mypack2.resize( mymap->getNumberOfArguments(), mymap->getNumberOfAtoms() );
103 3 : ref0->setupPCAStorage( mypack2 );
104 3 : mypack3.resize( mymap->getNumberOfArguments(), mymap->getNumberOfAtoms() );
105 16 : for(unsigned i=0; i<mymap->getNumberOfAtoms(); ++i) {
106 : mypack1.setAtomIndex(i,i);
107 : mypack2.setAtomIndex(i,i);
108 : mypack3.setAtomIndex(i,i);
109 : }
110 3 : mypack1_stashd_atoms.resize( mymap->getNumberOfAtoms() );
111 3 : mypack1_stashd_args.resize( mymap->getNumberOfArguments() );
112 3 : }
113 :
114 3 : std::string TrigonometricPathVessel::description() {
115 3 : return "values gspath and gzpath contain the position on and distance from the path calculated using trigonometry";
116 : }
117 :
118 7 : void TrigonometricPathVessel::resize() {
119 7 : StoreDataVessel::resize();
120 7 : if( getAction()->derivativesAreRequired() ) {
121 4 : unsigned nderivatives=getAction()->getNumberOfDerivatives();
122 4 : sp->resizeDerivatives( nderivatives );
123 4 : zp->resizeDerivatives( nderivatives );
124 : }
125 7 : }
126 :
127 1193 : void TrigonometricPathVessel::finish( const std::vector<double>& buffer ) {
128 : // Store the data calculated during mpi loop
129 1193 : StoreDataVessel::finish( buffer );
130 : // Get current value of all arguments
131 2487 : for(unsigned i=0; i<cargs.size(); ++i) {
132 1294 : cargs[i]=mymap->getArgument(i);
133 : }
134 :
135 : // Determine closest and second closest point to current position
136 1193 : double lambda=mymap->getLambda();
137 1193 : std::vector<double> dist( getNumberOfComponents() ), dist2( getNumberOfComponents() );;
138 1193 : retrieveSequentialValue( 0, false, dist );
139 1193 : retrieveSequentialValue( 1, false, dist2 );
140 1193 : iclose1=getStoreIndex(0);
141 1193 : iclose2=getStoreIndex(1);
142 1193 : double mindist1=dist[0], mindist2=dist2[0];
143 1193 : if( lambda>0.0 ) {
144 0 : mindist1=-std::log( dist[0] ) / lambda;
145 0 : mindist2=-std::log( dist2[0] ) / lambda;
146 : }
147 1193 : if( mindist2<mindist1 ) {
148 : double tmp=mindist1;
149 : mindist1=mindist2;
150 : mindist2=tmp;
151 866 : iclose1=getStoreIndex(1);
152 866 : iclose2=getStoreIndex(0);
153 : }
154 56519 : for(unsigned i=2; i<getNumberOfStoredValues(); ++i) {
155 55326 : retrieveSequentialValue( i, false, dist );
156 55326 : double ndist=dist[0];
157 55326 : if( lambda>0.0 ) {
158 0 : ndist=-std::log( dist[0] ) / lambda;
159 : }
160 55326 : if( ndist<mindist1 ) {
161 : mindist2=mindist1;
162 19737 : iclose2=iclose1;
163 : mindist1=ndist;
164 19737 : iclose1=getStoreIndex(i);
165 35589 : } else if( ndist<mindist2 ) {
166 : mindist2=ndist;
167 1062 : iclose2=getStoreIndex(i);
168 : }
169 : }
170 : // And find third closest point
171 1193 : int isign = iclose1 - iclose2;
172 : if( isign>1 ) {
173 : isign=1;
174 : } else if( isign<-1 ) {
175 : isign=-1;
176 : }
177 1193 : int iclose3 = iclose1 + isign;
178 : double v2v2;
179 : // We now have to compute vectors connecting the three closest points to the
180 : // new point
181 1193 : double v1v1 = (mymap->getReferenceConfiguration( iclose1 ))->calculate( mymap->getPositions(), mymap->getPbc(), mymap->getArguments(), mypack1, true );
182 1193 : double v3v3 = (mymap->getReferenceConfiguration( iclose2 ))->calculate( mymap->getPositions(), mymap->getPbc(), mymap->getArguments(), mypack3, true );
183 1193 : if( iclose3<0 || iclose3>=mymap->getFullNumberOfTasks() ) {
184 31 : ReferenceConfiguration* conf2=mymap->getReferenceConfiguration( iclose1 );
185 62 : v2v2=(mymap->getReferenceConfiguration( iclose2 ))->calc( conf2->getReferencePositions(), mymap->getPbc(), mymap->getArguments(),
186 31 : conf2->getReferenceArguments(), mypack2, true );
187 62 : (mymap->getReferenceConfiguration( iclose2 ))->extractDisplacementVector( conf2->getReferencePositions(), mymap->getArguments(),
188 62 : conf2->getReferenceArguments(), false, projdir );
189 : } else {
190 1162 : ReferenceConfiguration* conf2=mymap->getReferenceConfiguration( iclose3 );
191 2324 : v2v2=(mymap->getReferenceConfiguration( iclose1 ))->calc( conf2->getReferencePositions(), mymap->getPbc(), mymap->getArguments(),
192 1162 : conf2->getReferenceArguments(), mypack2, true );
193 2324 : (mymap->getReferenceConfiguration( iclose1 ))->extractDisplacementVector( conf2->getReferencePositions(), mymap->getArguments(),
194 2324 : conf2->getReferenceArguments(), false, projdir );
195 : }
196 :
197 : // Stash derivatives of v1v1
198 2487 : for(unsigned i=0; i<mymap->getNumberOfArguments(); ++i) {
199 1294 : mypack1_stashd_args[i]=mypack1.getArgumentDerivative(i);
200 : }
201 1193 : if( mymap->getNumberOfAtoms()>0 ) {
202 546 : ReferenceAtoms* at = dynamic_cast<ReferenceAtoms*>( mymap->getReferenceConfiguration( iclose1 ) );
203 : const std::vector<double> & displace( at->getDisplace() );
204 7644 : for(unsigned i=0; i<mymap->getNumberOfAtoms(); ++i) {
205 7098 : mypack1_stashd_atoms[i]=mypack1.getAtomDerivative(i);
206 7098 : mypack1.getAtomsDisplacementVector()[i] /= displace[i];
207 : }
208 : }
209 : // Calculate the dot product of v1 with v2
210 1193 : double v1v2 = (mymap->getReferenceConfiguration(iclose1))->projectDisplacementOnVector( projdir, mymap->getArguments(), cargs, mypack1 );
211 :
212 : // This computes s value
213 1193 : double spacing = mymap->getPropertyValue( iclose1, (mymap->property.begin())->first ) - mymap->getPropertyValue( iclose2, (mymap->property.begin())->first );
214 1193 : double root = sqrt( v1v2*v1v2 - v2v2 * ( v1v1 - v3v3) );
215 1193 : dx = 0.5 * ( (root + v1v2) / v2v2 - 1.);
216 1193 : double path_s = mymap->getPropertyValue(iclose1, (mymap->property.begin())->first ) + spacing * dx;
217 1193 : sp->set( path_s );
218 1193 : double fact = 0.25*spacing / v2v2;
219 : // Derivative of s wrt arguments
220 2487 : for(unsigned i=0; i<mymap->getNumberOfArguments(); ++i) {
221 1294 : sp->setDerivative( i, fact*( mypack2.getArgumentDerivative(i) + (v2v2 * (-mypack1_stashd_args[i] + mypack3.getArgumentDerivative(i))
222 1294 : + v1v2*mypack2.getArgumentDerivative(i) )/root ) );
223 : }
224 : // Derivative of s wrt atoms
225 1193 : unsigned narg=mymap->getNumberOfArguments();
226 1193 : Tensor vir;
227 1193 : vir.zero();
228 1193 : fact = 0.5*spacing / v2v2;
229 1193 : if( mymap->getNumberOfAtoms()>0 ) {
230 7644 : for(unsigned i=0; i<mymap->getNumberOfAtoms(); ++i) {
231 7098 : Vector ader = fact*(( v1v2*mypack1.getAtomDerivative(i) + 0.5*v2v2*(-mypack1_stashd_atoms[i] + mypack3.getAtomDerivative(i) ) )/root + mypack1.getAtomDerivative(i) );
232 28392 : for(unsigned k=0; k<3; ++k) {
233 21294 : sp->setDerivative( narg+3*i+k, ader[k] );
234 : }
235 7098 : vir-=Tensor( mymap->getPosition(i), ader );
236 : }
237 : // Set the virial
238 546 : unsigned nbase=narg+3*mymap->getNumberOfAtoms();
239 2184 : for(unsigned i=0; i<3; ++i)
240 6552 : for(unsigned j=0; j<3; ++j) {
241 4914 : sp->setDerivative( nbase+3*i+j, vir(i,j) );
242 : }
243 : }
244 : // Now compute z value
245 1193 : ReferenceConfiguration* conf2=mymap->getReferenceConfiguration( iclose1 );
246 2386 : double v4v4=(mymap->getReferenceConfiguration( iclose2 ))->calc( conf2->getReferencePositions(), mymap->getPbc(), mymap->getArguments(),
247 1193 : conf2->getReferenceArguments(), mypack2, true );
248 : // Extract vector connecting frames
249 2386 : (mymap->getReferenceConfiguration( iclose2 ))->extractDisplacementVector( conf2->getReferencePositions(), mymap->getArguments(),
250 1193 : conf2->getReferenceArguments(), false, projdir );
251 : // Calculate projection of vector on line connnecting frames
252 1193 : double proj = (mymap->getReferenceConfiguration(iclose1))->projectDisplacementOnVector( projdir, mymap->getArguments(), cargs, mypack1 );
253 1193 : double path_z = v1v1 + dx*dx*v4v4 - 2*dx*proj;
254 : // Derivatives for z path
255 1193 : path_z = sqrt(path_z);
256 1193 : zp->set( path_z );
257 1193 : vir.zero();
258 2487 : for(unsigned i=0; i<mymap->getNumberOfArguments(); ++i) {
259 1294 : zp->setDerivative( i, (mypack1_stashd_args[i] - 2*dx*mypack1.getArgumentDerivative(i))/(2.0*path_z) );
260 : }
261 : // Derivative wrt atoms
262 1193 : if( mymap->getNumberOfAtoms()>0 ) {
263 7644 : for(unsigned i=0; i<mymap->getNumberOfAtoms(); ++i) {
264 7098 : Vector dxder;
265 28392 : for(unsigned k=0; k<3; ++k) {
266 21294 : dxder[k] = ( 2*v4v4*dx - 2*proj )*spacing*sp->getDerivative( narg + 3*i+k );
267 : }
268 7098 : Vector ader = ( mypack1_stashd_atoms[i] - 2.*dx*mypack1.getAtomDerivative(i) + dxder )/ (2.0*path_z);
269 28392 : for(unsigned k=0; k<3; ++k) {
270 21294 : zp->setDerivative( narg+3*i+k, ader[k] );
271 : }
272 7098 : vir-=Tensor( mymap->getPosition(i), ader );
273 : }
274 : // Set the virial
275 546 : unsigned nbase=narg+3*mymap->getNumberOfAtoms();
276 2184 : for(unsigned i=0; i<3; ++i)
277 6552 : for(unsigned j=0; j<3; ++j) {
278 4914 : zp->setDerivative( nbase+3*i+j, vir(i,j) );
279 : }
280 : }
281 1193 : }
282 :
283 1193 : bool TrigonometricPathVessel::applyForce( std::vector<double>& forces ) {
284 1193 : std::vector<double> tmpforce( forces.size(), 0.0 );
285 1193 : forces.assign(forces.size(),0.0);
286 : bool wasforced=false;
287 1193 : if( sp->applyForce( tmpforce ) ) {
288 : wasforced=true;
289 0 : for(unsigned j=0; j<forces.size(); ++j) {
290 0 : forces[j]+=tmpforce[j];
291 : }
292 : }
293 1193 : tmpforce.assign(forces.size(),0.0);
294 1193 : if( zp->applyForce( tmpforce ) ) {
295 : wasforced=true;
296 0 : for(unsigned j=0; j<forces.size(); ++j) {
297 0 : forces[j]+=tmpforce[j];
298 : }
299 : }
300 1193 : return wasforced;
301 : }
302 :
303 : }
304 : }
|
