Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-2018 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/ActionWithValue.h"
      23             : #include "core/ActionPilot.h"
      24             : #include "core/ActionWithArguments.h"
      25             : #include "core/ActionRegister.h"
      26             : #include "tools/Matrix.h"
      27             : #include "PathProjectionCalculator.h"
      28             : 
      29             : //+PLUMEDOC ANALYSIS AVERAGE_PATH_DISPLACEMENT
      30             : /*
      31             : Accumulate the distances between the reference frames in the paths and the configurations visited
      32             : 
      33             : \par Examples
      34             : 
      35             : */
      36             : //+ENDPLUMEDOC
      37             : 
      38             : namespace PLMD {
      39             : namespace mapping {
      40             : 
      41             : class PathDisplacements : public ActionWithValue, public ActionPilot, public ActionWithArguments {
      42             : private:
      43             :   bool clearnextstep;
      44             :   unsigned clearstride;
      45             :   double fadefact;
      46             :   std::vector<double> wsum, displace_v;
      47             :   Matrix<double> displacements;
      48             :   PathProjectionCalculator path_projector;
      49             : public:
      50             :   static void registerKeywords( Keywords& keys );
      51             :   explicit PathDisplacements(const ActionOptions&);
      52             :   unsigned getNumberOfDerivatives();
      53         573 :   void clearDerivatives( const bool& force=false ) {}
      54         573 :   void calculate() {}
      55         573 :   void apply() {}
      56             :   void update();
      57             : };
      58             : 
      59             : PLUMED_REGISTER_ACTION(PathDisplacements,"AVERAGE_PATH_DISPLACEMENT")
      60             : 
      61           6 : void PathDisplacements::registerKeywords( Keywords& keys ) {
      62           6 :   Action::registerKeywords( keys );
      63           6 :   ActionWithValue::registerKeywords( keys );
      64           6 :   ActionPilot::registerKeywords( keys );
      65           6 :   ActionWithArguments::registerKeywords( keys );
      66           6 :   keys.use("ARG");
      67           6 :   PathProjectionCalculator::registerKeywords( keys );
      68          12 :   keys.add("compulsory","STRIDE","1","the frequency with which the average displacements should be collected and added to the average displacements");
      69          12 :   keys.add("compulsory","HALFLIFE","-1","the number of MD steps after which a previously measured path distance weighs only 50 percent in the average. This option may increase convergence by allowing to forget the memory of a bad initial guess path. The default is to set this to infinity");
      70          12 :   keys.add("compulsory","CLEAR","0","the frequency with which to clear all the accumulated data.  The default value "
      71             :            "of 0 implies that all the data will be used and that the grid will never be cleared");
      72           6 :   keys.setValueDescription("vector containing the average displacement between the trajectory and each of the landmarks that makes up the path");
      73           6 : }
      74             : 
      75           2 : PathDisplacements::PathDisplacements(const ActionOptions& ao):
      76             :   Action(ao),
      77             :   ActionWithValue(ao),
      78             :   ActionPilot(ao),
      79             :   ActionWithArguments(ao),
      80           2 :   clearnextstep(false),
      81           2 :   path_projector(this) {
      82             :   // Read in clear instructions
      83           2 :   parse("CLEAR",clearstride);
      84           2 :   if( clearstride>0 ) {
      85           2 :     if( clearstride%getStride()!=0 ) {
      86           0 :       error("CLEAR parameter must be a multiple of STRIDE");
      87             :     }
      88           2 :     log.printf("  clearing average every %u steps \n",clearstride);
      89             :   }
      90             :   double halflife;
      91           2 :   parse("HALFLIFE",halflife);
      92           2 :   log.printf("  weight of contribution to frame halves every %f steps \n",halflife);
      93           2 :   if( halflife<0 ) {
      94           2 :     fadefact=1.0;
      95             :   } else {
      96           0 :     fadefact = exp( -0.693147180559945 / static_cast<double>(halflife) );
      97             :   }
      98             :   // Now create the weights vector and displacements matrix
      99           2 :   unsigned nrows = getPntrToArgument(0)->getShape()[0];
     100           2 :   unsigned ncols = getPntrToArgument(0)->getShape()[1];
     101           2 :   wsum.resize( nrows );
     102             :   displacements.resize( nrows, ncols );
     103          64 :   for(unsigned i=0; i<nrows; ++i) {
     104          62 :     wsum[i]=0;
     105        1740 :     for(unsigned j=0; j<ncols; ++j) {
     106        1678 :       displacements(i,j)=0;
     107             :     }
     108             :   }
     109             :   // Add bibliography
     110           4 :   log<<"  Bibliography "<<plumed.cite("Diaz Leines and Ensing, Phys. Rev. Lett. 109, 020601 (2012)")<<"\n";
     111             :   // And create a value to hold the displacements
     112           2 :   std::vector<unsigned> shape(2);
     113           2 :   shape[0]=nrows;
     114           2 :   shape[1]=ncols;
     115           2 :   addValue( shape );
     116           2 :   setNotPeriodic();
     117           2 :   getPntrToComponent(0)->buildDataStore();
     118           2 :   getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
     119           2 : }
     120             : 
     121           0 : unsigned PathDisplacements::getNumberOfDerivatives() {
     122           0 :   return 0;
     123             : }
     124             : 
     125         573 : void PathDisplacements::update() {
     126         573 :   unsigned nrows = getPntrToArgument(0)->getShape()[0];
     127         573 :   unsigned ncols = getPntrToArgument(0)->getShape()[1];
     128             : 
     129         573 :   if( clearnextstep ) {
     130             :     unsigned k=0;
     131        1010 :     for(unsigned i=0; i<nrows; ++i) {
     132       38700 :       for(unsigned j=0; j<ncols; ++j) {
     133       37714 :         displacements(i,j)=0;
     134       37714 :         getPntrToComponent(0)->set(k,0);
     135       37714 :         k++;
     136             :       }
     137             :     }
     138          24 :     clearnextstep=false;
     139             :   }
     140             : 
     141             :   unsigned k=0, iclose1=0, iclose2=0;
     142             :   double v1v1=0, v3v3=0;
     143       22417 :   for(unsigned i=0; i<nrows; ++i) {
     144             :     double dist = 0;
     145      799020 :     for(unsigned j=0; j<ncols; ++j) {
     146      777176 :       double tmp = getPntrToArgument(0)->get(k);
     147      777176 :       dist += tmp*tmp;
     148      777176 :       k++;
     149             :     }
     150       21844 :     if( i==0 ) {
     151             :       v1v1 = dist;
     152             :       iclose1 = 0;
     153       21271 :     } else if( dist<v1v1 ) {
     154             :       v3v3=v1v1;
     155             :       v1v1=dist;
     156             :       iclose2=iclose1;
     157             :       iclose1=i;
     158       10991 :     } else if( i==1 ) {
     159             :       v3v3=dist;
     160             :       iclose2=1;
     161       10991 :     } else if( dist<v3v3 ) {
     162             :       v3v3=dist;
     163             :       iclose2=i;
     164             :     }
     165             :   }
     166             :   // And find third closest point
     167         573 :   int isign = iclose1 - iclose2;
     168             :   if( isign>1 ) {
     169             :     isign=1;
     170             :   } else if( isign<-1 ) {
     171             :     isign=-1;
     172             :   }
     173         573 :   int iclose3 = iclose1 + isign;
     174         573 :   unsigned ifrom=iclose1, ito=iclose3;
     175         573 :   if( iclose3<0 || iclose3>=nrows ) {
     176           0 :     ifrom=iclose2;
     177           0 :     ito=iclose1;
     178             :   }
     179             : 
     180             :   // Calculate the dot product of v1 with v2
     181         573 :   path_projector.getDisplaceVector( ifrom, ito, displace_v );
     182             :   double v2v2=0, v1v2=0;
     183         573 :   unsigned kclose1 = iclose1*ncols;
     184       19183 :   for(unsigned i=0; i<displace_v.size(); ++i) {
     185       18610 :     v2v2 += displace_v[i]*displace_v[i];
     186       18610 :     v1v2 += displace_v[i]*getPntrToArgument(0)->get(kclose1+i);
     187             :   }
     188             : 
     189         573 :   double root = sqrt( v1v2*v1v2 - v2v2 * ( v1v1 - v3v3) );
     190         573 :   double dx = 0.5 * ( (root + v1v2) / v2v2 - 1.);
     191         573 :   double weight2 = -1.* dx;
     192         573 :   double weight1 = 1.0 + dx;
     193         573 :   if( weight1>1.0 ) {
     194             :     weight1=1.0;
     195             :     weight2=0.0;
     196         573 :   } else if( weight2>1.0 ) {
     197             :     weight1=0.0;
     198             :     weight2=1.0;
     199             :   }
     200             : 
     201             :   // Accumulate displacements for path
     202       19183 :   for(unsigned i=0; i<ncols; ++i) {
     203       18610 :     double displace = getPntrToArgument(0)->get(kclose1+i) - dx*displace_v[i];
     204       18610 :     displacements(iclose1,i) += weight1 * displace;
     205       18610 :     displacements(iclose2,i) += weight2 * displace;
     206             :   }
     207             : 
     208             :   // Update weight accumulators
     209         573 :   wsum[iclose1] *= fadefact;
     210         573 :   wsum[iclose2] *= fadefact;
     211         573 :   wsum[iclose1] += weight1;
     212         573 :   wsum[iclose2] += weight2;
     213             : 
     214             :   // Update numbers in values
     215         573 :   if( wsum[iclose1] > epsilon ) {
     216       19183 :     for(unsigned i=0; i<ncols; ++i) {
     217       18610 :       getPntrToComponent(0)->set( kclose1+i, displacements(iclose1,i) / wsum[iclose1] );
     218             :     }
     219             :   }
     220         573 :   if( wsum[iclose2] > epsilon ) {
     221         573 :     unsigned kclose2 = iclose2*ncols;
     222       19183 :     for(unsigned i=0; i<ncols; ++i) {
     223       18610 :       getPntrToComponent(0)->set( kclose2+i, displacements(iclose2,i) / wsum[iclose2] );
     224             :     }
     225             :   }
     226             : 
     227             :   // Clear if required
     228         573 :   if( (getStep()>0 && clearstride>0 && getStep()%clearstride==0) ) {
     229          25 :     clearnextstep=true;
     230             :   }
     231         573 : }
     232             : 
     233             : }
     234             : }