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