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 "TrigonometricPathVessel.h"
23 : #include "vesselbase/VesselRegister.h"
24 :
25 : namespace PLMD {
26 : namespace mapping {
27 :
28 6455 : PLUMED_REGISTER_VESSEL(TrigonometricPathVessel,"GPATH")
29 :
30 3 : void TrigonometricPathVessel::registerKeywords( Keywords& keys ) {
31 3 : StoreDataVessel::registerKeywords(keys);
32 3 : }
33 :
34 1613 : void TrigonometricPathVessel::reserveKeyword( Keywords& keys ) {
35 6452 : keys.reserve("vessel","GPATH","calculate the position on the path using trigonometry");
36 6452 : keys.addOutputComponent("gspath","GPATH","the position on the path calculated using trigonometry");
37 6452 : keys.addOutputComponent("gzpath","GPATH","the distance from the path calculated using trigonometry");
38 1613 : }
39 :
40 3 : TrigonometricPathVessel::TrigonometricPathVessel( const vesselbase::VesselOptions& da ):
41 : StoreDataVessel(da),
42 9 : projdir(ReferenceConfigurationOptions("DIRECTION")),
43 6 : mydpack1( 1, getAction()->getNumberOfDerivatives() ),
44 6 : mydpack2( 1, getAction()->getNumberOfDerivatives() ),
45 6 : mydpack3( 1, getAction()->getNumberOfDerivatives() ),
46 : mypack1( 0, 0, mydpack1 ),
47 : mypack2( 0, 0, mydpack2 ),
48 36 : mypack3( 0, 0, mydpack3 )
49 : {
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 ) error("cannot use trigonometric paths when there are multiple properties");
54 :
55 372 : for(unsigned i=0; i<mymap->getFullNumberOfTasks(); ++i) {
56 122 : if( mymap->getTaskCode(i)!=mymap->getPositionInFullTaskList(i) ) error("mismatched tasks and codes");
57 : }
58 9 : mymap->addComponentWithDerivatives("gspath"); mymap->componentIsNotPeriodic("gspath");
59 6 : sp=mymap->copyOutput( mymap->getNumberOfComponents()-1 ); sp->resizeDerivatives( mymap->getNumberOfDerivatives() );
60 9 : mymap->addComponentWithDerivatives("gzpath"); mymap->componentIsNotPeriodic("gzpath");
61 6 : zp=mymap->copyOutput( mymap->getNumberOfComponents()-1 ); zp->resizeDerivatives( mymap->getNumberOfDerivatives() );
62 :
63 : // Check we have PCA
64 3 : ReferenceConfiguration* ref0=mymap->getReferenceConfiguration(0);
65 250 : for(unsigned i=0; i<mymap->getFullNumberOfTasks(); ++i) {
66 122 : if( !(mymap->getReferenceConfiguration(i))->pcaIsEnabledForThisReference() ) error("pca must be implemented in order to use trigometric path");
67 366 : if( ref0->getName()!=(mymap->getReferenceConfiguration(i))->getName() ) error("cannot use mixed metrics");
68 244 : if( mymap->getNumberOfAtoms()!=(mymap->getReferenceConfiguration(i))->getNumberOfReferencePositions() ) error("all frames must use the same set of atoms");
69 122 : if( mymap->getNumberOfArguments()!=(mymap->getReferenceConfiguration(i))->getNumberOfReferenceArguments() ) error("all frames must use the same set of arguments");
70 : }
71 :
72 6 : cargs.resize( mymap->getNumberOfArguments() ); std::vector<std::string> argument_names( mymap->getNumberOfArguments() );
73 11 : for(unsigned i=0; i<mymap->getNumberOfArguments(); ++i) argument_names[i] = (mymap->getPntrToArgument(i))->getName();
74 3 : projdir.setNamesAndAtomNumbers( mymap->getAbsoluteIndexes(), argument_names );
75 6 : mypack1.resize( mymap->getNumberOfArguments(), mymap->getNumberOfAtoms() ); ref0->setupPCAStorage( mypack1 );
76 6 : mypack2.resize( mymap->getNumberOfArguments(), mymap->getNumberOfAtoms() ); ref0->setupPCAStorage( mypack2 );
77 6 : mypack3.resize( mymap->getNumberOfArguments(), mymap->getNumberOfAtoms() );
78 45 : for(unsigned i=0; i<mymap->getNumberOfAtoms(); ++i) { mypack1.setAtomIndex(i,i); mypack2.setAtomIndex(i,i); mypack3.setAtomIndex(i,i); }
79 3 : mypack1_stashd_atoms.resize( mymap->getNumberOfAtoms() ); mypack1_stashd_args.resize( mymap->getNumberOfArguments() );
80 3 : }
81 :
82 3 : std::string TrigonometricPathVessel::description() {
83 3 : return "values gspath and gzpath contain the position on and distance from the path calculated using trigonometry";
84 : }
85 :
86 7 : void TrigonometricPathVessel::resize() {
87 7 : StoreDataVessel::resize();
88 7 : if( getAction()->derivativesAreRequired() ) {
89 4 : unsigned nderivatives=getAction()->getNumberOfDerivatives();
90 8 : sp->resizeDerivatives( nderivatives ); zp->resizeDerivatives( nderivatives );
91 : }
92 7 : }
93 :
94 1193 : void TrigonometricPathVessel::finish( const std::vector<double>& buffer ) {
95 : // Store the data calculated during mpi loop
96 1193 : StoreDataVessel::finish( buffer );
97 : // Get current value of all arguments
98 6268 : for(unsigned i=0; i<cargs.size(); ++i) cargs[i]=mymap->getArgument(i);
99 :
100 : // Determine closest and second closest point to current position
101 1193 : double lambda=mymap->getLambda();
102 2386 : std::vector<double> dist( getNumberOfComponents() ), dist2( getNumberOfComponents() );;
103 1193 : retrieveSequentialValue( 0, false, dist );
104 1193 : retrieveSequentialValue( 1, false, dist2 );
105 1193 : iclose1=getStoreIndex(0); iclose2=getStoreIndex(1);
106 2386 : double mindist1=dist[0], mindist2=dist2[0];
107 1193 : if( lambda>0.0 ) {
108 0 : mindist1=-std::log( dist[0] ) / lambda;
109 0 : mindist2=-std::log( dist2[0] ) / lambda;
110 : }
111 1193 : if( mindist2<mindist1 ) {
112 : double tmp=mindist1; mindist1=mindist2; mindist2=tmp;
113 866 : iclose1=getStoreIndex(1); iclose2=getStoreIndex(0);
114 : }
115 56519 : for(unsigned i=2; i<getNumberOfStoredValues(); ++i) {
116 55326 : retrieveSequentialValue( i, false, dist );
117 55326 : double ndist=dist[0];
118 55326 : if( lambda>0.0 ) ndist=-std::log( dist[0] ) / lambda;
119 55326 : if( ndist<mindist1 ) {
120 19737 : mindist2=mindist1; iclose2=iclose1;
121 19737 : mindist1=ndist; iclose1=getStoreIndex(i);
122 35589 : } else if( ndist<mindist2 ) {
123 1062 : mindist2=ndist; iclose2=getStoreIndex(i);
124 : }
125 : }
126 : // And find third closest point
127 1193 : int isign = iclose1 - iclose2;
128 1193 : if( isign>1 ) isign=1; else if( isign<-1 ) isign=-1;
129 1193 : int iclose3 = iclose1 + isign; double v2v2;
130 : // We now have to compute vectors connecting the three closest points to the
131 : // new point
132 2386 : double v1v1 = (mymap->getReferenceConfiguration( iclose1 ))->calculate( mymap->getPositions(), mymap->getPbc(), mymap->getArguments(), mypack1, true );
133 2386 : double v3v3 = (mymap->getReferenceConfiguration( iclose2 ))->calculate( mymap->getPositions(), mymap->getPbc(), mymap->getArguments(), mypack3, true );
134 2369 : if( iclose3<0 || iclose3>=mymap->getFullNumberOfTasks() ) {
135 31 : ReferenceConfiguration* conf2=mymap->getReferenceConfiguration( iclose1 );
136 124 : v2v2=(mymap->getReferenceConfiguration( iclose2 ))->calc( conf2->getReferencePositions(), mymap->getPbc(), mymap->getArguments(),
137 62 : conf2->getReferenceArguments(), mypack2, true );
138 93 : (mymap->getReferenceConfiguration( iclose2 ))->extractDisplacementVector( conf2->getReferencePositions(), mymap->getArguments(),
139 31 : conf2->getReferenceArguments(), false, projdir );
140 : } else {
141 2324 : ReferenceConfiguration* conf2=mymap->getReferenceConfiguration( iclose3 );
142 4648 : v2v2=(mymap->getReferenceConfiguration( iclose1 ))->calc( conf2->getReferencePositions(), mymap->getPbc(), mymap->getArguments(),
143 2324 : conf2->getReferenceArguments(), mypack2, true );
144 3486 : (mymap->getReferenceConfiguration( iclose1 ))->extractDisplacementVector( conf2->getReferencePositions(), mymap->getArguments(),
145 1162 : conf2->getReferenceArguments(), false, projdir );
146 : }
147 :
148 : // Stash derivatives of v1v1
149 3781 : for(unsigned i=0; i<mymap->getNumberOfArguments(); ++i) mypack1_stashd_args[i]=mypack1.getArgumentDerivative(i);
150 2386 : if( mymap->getNumberOfAtoms()>0 ) {
151 546 : ReferenceAtoms* at = dynamic_cast<ReferenceAtoms*>( mymap->getReferenceConfiguration( iclose1 ) );
152 : const std::vector<double> & displace( at->getDisplace() );
153 15288 : for(unsigned i=0; i<mymap->getNumberOfAtoms(); ++i) {
154 28392 : mypack1_stashd_atoms[i]=mypack1.getAtomDerivative(i); mypack1.getAtomsDisplacementVector()[i] /= displace[i];
155 : }
156 : }
157 : // Calculate the dot product of v1 with v2
158 2386 : double v1v2 = (mymap->getReferenceConfiguration(iclose1))->projectDisplacementOnVector( projdir, mymap->getArguments(), cargs, mypack1 );
159 :
160 : // This computes s value
161 2386 : double spacing = mymap->getPropertyValue( iclose1, 0 ) - mymap->getPropertyValue( iclose2, 0 );
162 1193 : double root = sqrt( v1v2*v1v2 - v2v2 * ( v1v1 - v3v3) );
163 1193 : dx = 0.5 * ( (root + v1v2) / v2v2 - 1.);
164 2386 : double path_s = mymap->getPropertyValue(iclose1, 0) + spacing * dx; sp->set( path_s );
165 1193 : double fact = 0.25*spacing / v2v2;
166 : // Derivative of s wrt arguments
167 3781 : for(unsigned i=0; i<mymap->getNumberOfArguments(); ++i) {
168 5176 : sp->setDerivative( i, fact*( mypack2.getArgumentDerivative(i) + (v2v2 * (-mypack1_stashd_args[i] + mypack3.getArgumentDerivative(i))
169 1294 : + v1v2*mypack2.getArgumentDerivative(i) )/root ) );
170 : }
171 : // Derivative of s wrt atoms
172 1193 : unsigned narg=mymap->getNumberOfArguments(); Tensor vir; vir.zero(); fact = 0.5*spacing / v2v2;
173 2386 : if( mymap->getNumberOfAtoms()>0 ) {
174 15288 : for(unsigned i=0; i<mymap->getNumberOfAtoms(); ++i) {
175 14196 : Vector ader = fact*(( v1v2*mypack1.getAtomDerivative(i) + 0.5*v2v2*(-mypack1_stashd_atoms[i] + mypack3.getAtomDerivative(i) ) )/root + mypack1.getAtomDerivative(i) );
176 49686 : for(unsigned k=0; k<3; ++k) sp->setDerivative( narg+3*i+k, ader[k] );
177 14196 : vir-=Tensor( mymap->getPosition(i), ader );
178 : }
179 : // Set the virial
180 546 : unsigned nbase=narg+3*mymap->getNumberOfAtoms();
181 7098 : for(unsigned i=0; i<3; ++i) for(unsigned j=0; j<3; ++j) sp->setDerivative( nbase+3*i+j, vir(i,j) );
182 : }
183 :
184 : // Now compute z value
185 1193 : ReferenceConfiguration* conf2=mymap->getReferenceConfiguration( iclose1 );
186 5965 : double v4v4=(mymap->getReferenceConfiguration( iclose2 ))->calc( conf2->getReferencePositions(), mymap->getPbc(), mymap->getArguments(),
187 3579 : conf2->getReferenceArguments(), mypack2, true );
188 : // Extract vector connecting frames
189 3579 : (mymap->getReferenceConfiguration( iclose2 ))->extractDisplacementVector( conf2->getReferencePositions(), mymap->getArguments(),
190 1193 : conf2->getReferenceArguments(), false, projdir );
191 : // Calculate projection of vector on line connnecting frames
192 2386 : double proj = (mymap->getReferenceConfiguration(iclose1))->projectDisplacementOnVector( projdir, mymap->getArguments(), cargs, mypack1 );
193 1193 : double path_z = v1v1 + dx*dx*v4v4 - 2*dx*proj;
194 :
195 : // Derivatives for z path
196 2386 : path_z = sqrt(path_z); zp->set( path_z ); vir.zero();
197 6369 : for(unsigned i=0; i<mymap->getNumberOfArguments(); ++i) zp->setDerivative( i, (mypack1_stashd_args[i] - 2*dx*mypack1.getArgumentDerivative(i))/(2.0*path_z) );
198 : // Derivative wrt atoms
199 2386 : if( mymap->getNumberOfAtoms()>0 ) {
200 15288 : for(unsigned i=0; i<mymap->getNumberOfAtoms(); ++i) {
201 28392 : Vector dxder; for(unsigned k=0; k<3; ++k) dxder[k] = ( 2*v4v4*dx - 2*proj )*spacing*sp->getDerivative( narg + 3*i+k );
202 14196 : Vector ader = ( mypack1_stashd_atoms[i] - 2.*dx*mypack1.getAtomDerivative(i) + dxder )/ (2.0*path_z);
203 49686 : for(unsigned k=0; k<3; ++k) zp->setDerivative( narg+3*i+k, ader[k] );
204 14196 : vir-=Tensor( mymap->getPosition(i), ader );
205 : }
206 : // Set the virial
207 546 : unsigned nbase=narg+3*mymap->getNumberOfAtoms();
208 7098 : for(unsigned i=0; i<3; ++i) for(unsigned j=0; j<3; ++j) zp->setDerivative( nbase+3*i+j, vir(i,j) );
209 : }
210 1193 : }
211 :
212 1193 : bool TrigonometricPathVessel::applyForce( std::vector<double>& forces ) {
213 2386 : std::vector<double> tmpforce( forces.size(), 0.0 );
214 2386 : forces.assign(forces.size(),0.0); bool wasforced=false;
215 1193 : if( sp->applyForce( tmpforce ) ) {
216 : wasforced=true;
217 0 : for(unsigned j=0; j<forces.size(); ++j) forces[j]+=tmpforce[j];
218 : }
219 2386 : tmpforce.assign(forces.size(),0.0);
220 1193 : if( zp->applyForce( tmpforce ) ) {
221 : wasforced=true;
222 0 : for(unsigned j=0; j<forces.size(); ++j) forces[j]+=tmpforce[j];
223 : }
224 1193 : return wasforced;
225 : }
226 :
227 : }
228 4839 : }
|
