Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-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 "RMSDVector.h"
      23             : #include "core/PlumedMain.h"
      24             : #include "core/ActionRegister.h"
      25             : #include "tools/PDB.h"
      26             : 
      27             : //+PLUMEDOC DCOLVAR RMSD_VECTOR
      28             : /*
      29             : Calculate the RMSD distance between the instaneous configuration and multiple reference configurations
      30             : 
      31             : 
      32             : \par Examples
      33             : 
      34             : */
      35             : //+ENDPLUMEDOC
      36             : 
      37             : namespace PLMD {
      38             : namespace colvar {
      39             : 
      40             : PLUMED_REGISTER_ACTION(RMSDVector,"RMSD_VECTOR")
      41             : 
      42          81 : void RMSDVector::registerKeywords(Keywords& keys) {
      43         162 :   ActionWithVector::registerKeywords(keys); keys.use("ARG"); keys.setDisplayName("RMSD");
      44         162 :   keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed.  Should be OPTIMAL or SIMPLE.");
      45         162 :   keys.add("compulsory","ALIGN","1.0","the weights to use when aligning to the reference structure");
      46         162 :   keys.add("compulsory","DISPLACE","1.0","the weights to use when calculating the displacement from the reference structure");
      47         162 :   keys.addFlag("SQUARED",false," This should be set if you want mean squared displacement instead of RMSD ");
      48         162 :   keys.addFlag("UNORMALIZED",false,"by default the mean sequare deviation or root mean square deviation is calculated.  If this option is given no averaging is done");
      49         162 :   keys.addFlag("DISPLACEMENT",false,"Calculate the vector of displacements instead of the length of this vector");
      50         162 :   keys.addOutputComponent("disp","DISPLACEMENT","the vector of displacements for the atoms");
      51         162 :   keys.addOutputComponent("dist","DISPLACEMENT","the RMSD distance the atoms have moved");
      52          81 :   keys.setValueDescription("a vector containing the RMSD between the instantaneous structure and each of the reference structures that were input");
      53          81 : }
      54             : 
      55          49 : RMSDVector::RMSDVector(const ActionOptions&ao):
      56             :   Action(ao),
      57             :   ActionWithVector(ao),
      58          49 :   firststep(true)
      59             : {
      60          49 :   if( getPntrToArgument(0)->getRank()!=1 ) error("first argument should be vector");
      61          49 :   if( getPntrToArgument(1)->getRank()<1 ) error("second argument should be matrix or a vector");
      62          49 :   if( getPntrToArgument(1)->getRank()==1 ) {
      63           7 :     if( getPntrToArgument(0)->getNumberOfValues()!=getPntrToArgument(1)->getNumberOfValues() ) error("mismatch between sizes of input vectors");
      64          42 :   } else if( getPntrToArgument(1)->getRank()==2 ) {
      65          42 :     if( getPntrToArgument(0)->getNumberOfValues()!=getPntrToArgument(1)->getShape()[1] ) error("mismatch between sizes of input vectors");
      66             :   }
      67          49 :   if( getPntrToArgument(0)->getNumberOfValues()%3!=0 ) error("number of components in input arguments should be multiple of three");
      68             : 
      69          49 :   unsigned natoms = getPntrToArgument(0)->getNumberOfValues() / 3;
      70          98 :   type.assign("SIMPLE"); parse("TYPE",type); parseFlag("SQUARED",squared);
      71          49 :   align.resize( natoms ); parseVector("ALIGN",align);
      72          49 :   displace.resize( natoms ); parseVector("DISPLACE",displace);
      73          98 :   bool unorm=false; parseFlag("UNORMALIZED",unorm); norm_weights=!unorm;
      74          49 :   double wa=0, wd=0; sqrtdisplace.resize( displace.size() );
      75         630 :   for(unsigned i=0; i<align.size(); ++i) { wa+=align[i]; wd+=displace[i]; }
      76             : 
      77          49 :   if( wa>epsilon ) {
      78          49 :     double iwa = 1. / wa;
      79         630 :     for(unsigned i=0; i<align.size(); ++i) align[i] *= iwa;
      80             :   } else {
      81           0 :     double iwa = 1. / natoms;
      82           0 :     for(unsigned i=0; i<align.size(); ++i) align[i] = iwa;
      83             :   }
      84          49 :   if( wd>epsilon ) {
      85          49 :     if( !norm_weights ) { wd = 1; } double iwd = 1. / wd;
      86         630 :     for(unsigned i=0; i<align.size(); ++i) displace[i] *= iwd;
      87             :   } else {
      88           0 :     double iwd = 1. / natoms;
      89           0 :     for(unsigned i=0; i<align.size(); ++i) displace[i] = iwd;
      90             :   }
      91         630 :   for(unsigned i=0; i<align.size(); ++i) sqrtdisplace[i] = sqrt(displace[i]);
      92             : 
      93          49 :   parseFlag("DISPLACEMENT",displacement);
      94          49 :   if( displacement && (getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getShape()[0]<=1) ) {
      95          78 :     addComponentWithDerivatives("dist"); componentIsNotPeriodic("dist");
      96          26 :     std::vector<unsigned> shape( 1, getPntrToArgument(0)->getNumberOfValues() );
      97          78 :     addComponent( "disp", shape ); getPntrToComponent(1)->buildDataStore(); componentIsNotPeriodic("disp");
      98          23 :   } else if( displacement ) {
      99           2 :     std::vector<unsigned> shape( 1, getPntrToArgument(1)->getShape()[0] );
     100           4 :     addComponent( "dist", shape ); getPntrToComponent(0)->buildDataStore();
     101           2 :     componentIsNotPeriodic("dist");
     102           2 :     shape.resize(2); shape[0] = getPntrToArgument(1)->getShape()[0]; shape[1] = getPntrToArgument(0)->getNumberOfValues();
     103           4 :     addComponent( "disp", shape ); getPntrToComponent(1)->buildDataStore(); getPntrToComponent(1)->reshapeMatrixStore( shape[1] );
     104           4 :     componentIsNotPeriodic("disp");
     105          21 :   } else if( (getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getShape()[0]==1) ) {
     106          26 :     addValue(); setNotPeriodic();
     107             :   } else {
     108           8 :     std::vector<unsigned> shape( 1, getPntrToArgument(1)->getShape()[0] );
     109           8 :     addValue( shape ); setNotPeriodic();
     110             :   }
     111          49 :   if( getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getNumberOfValues()==0 ) myrmsd.resize(1);
     112          42 :   else myrmsd.resize( getPntrToArgument(1)->getShape()[0] );
     113             : 
     114          49 :   if( getPntrToArgument(1)->getRank()==1 ) log.printf("  calculating RMSD distance between %d atoms. Distance between the avectors of atoms in %s and %s\n",
     115             :         natoms, getPntrToArgument(0)->getName().c_str(), getPntrToArgument(1)->getName().c_str() );
     116          42 :   else log.printf("  calculating RMSD distance between %d sets of %d atoms. Distance between vector %s of atoms and matrix of configurations in %s\n",
     117             :                     getPntrToArgument(1)->getShape()[0], natoms, getPntrToArgument(0)->getName().c_str(), getPntrToArgument(1)->getName().c_str() );
     118          49 :   log.printf("  method for alignment : %s \n",type.c_str() );
     119          49 :   if(squared)log.printf("  chosen to use SQUARED option for MSD instead of RMSD\n");
     120          21 :   else      log.printf("  using periodic boundary conditions\n");
     121          49 : }
     122             : 
     123         804 : unsigned RMSDVector::getNumberOfDerivatives() {
     124         804 :   return getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues();
     125             : }
     126             : 
     127       10631 : void RMSDVector::setReferenceConfigurations() {
     128       10631 :   unsigned natoms = getPntrToArgument(0)->getShape()[0] / 3;
     129       10631 :   Vector center; std::vector<Vector> pos( natoms );
     130       22578 :   for(unsigned jconf=0; jconf<myrmsd.size(); ++jconf) {
     131       11947 :     center.zero();
     132      172431 :     for(unsigned i=0; i<pos.size(); ++i) {
     133      641936 :       for(unsigned j=0; j<3; ++j) pos[i][j] = getPntrToArgument(1)->get( (3*jconf+j)*pos.size() + i );
     134      160484 :       center+=pos[i]*align[i];
     135             :     }
     136      172431 :     for(unsigned i=0; i<pos.size(); ++i) pos[i] -= center;
     137       11947 :     myrmsd[jconf].clear(); myrmsd[jconf].set(align,displace,pos,type,true,norm_weights);
     138             :   }
     139       10631 : }
     140             : 
     141      192996 : double RMSDVector::calculateRMSD( const unsigned& current, std::vector<Vector>& pos, std::vector<Vector>& der, std::vector<Vector>& direction ) const {
     142      192996 :   unsigned natoms = pos.size();
     143     2704147 :   for(unsigned i=0; i<natoms; ++i) {
     144    10044604 :     for(unsigned j=0; j<3; ++j) pos[i][j] = getPntrToArgument(0)->get( j*natoms + i );
     145             :   }
     146             : 
     147      192996 :   if( displacement && type=="SIMPLE" ) {
     148         646 :     const Value* myval = getConstPntrToComponent(1);
     149         646 :     double r = myrmsd[current].simpleAlignment( align, displace, pos, myrmsd[current].getReference(), der, direction, squared );
     150         646 :     if( !doNotCalculateDerivatives() && myval->forcesWereAdded() ) {
     151          33 :       Vector comforce; comforce.zero();
     152         187 :       for(unsigned i=0; i<natoms; i++) {
     153         616 :         for(unsigned k=0; k<3; ++k) comforce[k] += align[i]*myval->getForce( (3*current+k)*natoms + i);
     154             :       }
     155         187 :       for(unsigned i=0; i<natoms; i++) {
     156         616 :         for(unsigned k=0; k<3; ++k) direction[i][k] = myval->getForce( (3*current+k)*natoms + i ) - comforce[k];
     157             :       }
     158             :     }
     159         646 :     return r;
     160      192350 :   } else if( displacement ) {
     161       55589 :     const Value* myval = getConstPntrToComponent(1);
     162      111178 :     Tensor rot; Matrix<std::vector<Vector> > DRotDPos(3,3); std::vector<Vector> centeredpos( natoms ), centeredreference( natoms );
     163       55589 :     double r = myrmsd[current].calc_PCAelements( pos, der, rot, DRotDPos, direction, centeredpos, centeredreference, squared ); std::vector<Vector> ref( myrmsd[current].getReference() );
     164       55589 :     if( !doNotCalculateDerivatives() && myval->forcesWereAdded() ) {
     165        1671 :       Tensor trot=rot.transpose(); double prefactor = 1 / static_cast<double>( natoms ); Vector v1; v1.zero();
     166       22573 :       for(unsigned n=0; n<natoms; n++) {
     167       83608 :         Vector ff; for(unsigned k=0; k<3; ++k ) ff[k] = myval->getForce( (3*current+k)*natoms + n );
     168       20902 :         v1+=prefactor*matmul(trot,ff);
     169             :       }
     170             :       // Notice that we use centreredreference here to accumulate the true forces
     171       22573 :       for(unsigned n=0; n<natoms; n++) {
     172       83608 :         Vector ff; for(unsigned k=0; k<3; ++k ) ff[k] = myval->getForce( (3*current+k)*natoms + n );
     173       20902 :         centeredreference[n] = sqrtdisplace[n]*( matmul(trot,ff) - v1 );
     174             :       }
     175        6684 :       for(unsigned a=0; a<3; a++) {
     176       20052 :         for(unsigned b=0; b<3; b++) {
     177      203157 :           double tmp1=0.; for(unsigned m=0; m<natoms; m++) tmp1+=centeredpos[m][b]*myval->getForce( (3*current+a)*natoms + m );
     178      203157 :           for(unsigned i=0; i<natoms; i++) centeredreference[i] += sqrtdisplace[i]*tmp1*DRotDPos[a][b][i];
     179             :         }
     180             :       }
     181             :       // Now subtract the current force and add on the true force
     182       22573 :       for(unsigned n=0; n<natoms; n++) {
     183       83608 :         for(unsigned k=0; k<3; ++k) direction[n][k] = centeredreference[n][k];
     184             :       }
     185             :     } else {
     186      759226 :       for(unsigned i=0; i<direction.size(); ++i) direction[i] = sqrtdisplace[i]*( direction[i] - ref[i] );
     187             :     }
     188             :     return r;
     189             :   }
     190      136761 :   return myrmsd[current].calculate( pos, der, squared );
     191             : }
     192             : 
     193             : // calculator
     194       17182 : void RMSDVector::calculate() {
     195       17182 :   if( firststep || !getPntrToArgument(1)->isConstant() ) { setReferenceConfigurations(); firststep=false; }
     196             : 
     197       17182 :   if( getPntrToComponent(0)->getRank()==0 ) {
     198       11796 :     unsigned natoms = getPntrToArgument(0)->getShape()[0] / 3;
     199       11796 :     std::vector<Vector> pos( natoms ), der( natoms ), direction( natoms );
     200       11796 :     double r = calculateRMSD( 0, pos, der, direction );
     201             : 
     202       11796 :     getPntrToComponent(0)->set( r );
     203       11796 :     if( getNumberOfComponents()==2 ) {
     204       11775 :       Value* mydisp = getPntrToComponent(1);
     205      168204 :       for(unsigned i=0; i<natoms; i++) {
     206      625716 :         for(unsigned j=0; j<3; ++j ) mydisp->set( j*natoms+i, direction[i][j] );
     207             :       }
     208             :     }
     209       11796 :     if( doNotCalculateDerivatives() ) return;
     210             : 
     211         612 :     Value* myval = getPntrToComponent(0);
     212        7472 :     for(unsigned i=0; i<natoms; i++) {
     213       27440 :       for(unsigned j=0; j<3; ++j ) myval->setDerivative( j*natoms+i, der[i][j] );
     214             :     }
     215        5386 :   } else runAllTasks();
     216             : }
     217             : 
     218      160524 : bool RMSDVector::checkForTaskForce( const unsigned& itask, const Value* myval ) const {
     219      160524 :   if( myval->getRank()<2 ) return ActionWithVector::checkForTaskForce( itask, myval );
     220       22932 :   unsigned nelements = myval->getShape()[1], startr = itask*nelements;
     221      874692 :   for(unsigned j=0; j<nelements; ++j ) {
     222      852852 :     if( fabs( myval->getForce( startr + j ) )>epsilon ) return true;
     223             :   }
     224             :   return false;
     225             : }
     226             : 
     227       17162 : void RMSDVector::apply() {
     228       17162 :   if( doNotCalculateDerivatives() ) return;
     229             : 
     230        4980 :   if( getPntrToComponent(0)->getRank()==0 ) {
     231         612 :     std::vector<double> forces( getNumberOfDerivatives(), 0 );
     232         612 :     bool wasforced = getPntrToComponent(0)->applyForce( forces );
     233             : 
     234         612 :     if( getNumberOfComponents()==2 && getPntrToComponent(1)->forcesWereAdded() ) {
     235         612 :       unsigned natoms = getPntrToArgument(0)->getShape()[0] / 3;
     236         612 :       std::vector<Vector> pos( natoms ), der( natoms ), direction( natoms );
     237         612 :       double r = calculateRMSD( 0, pos, der, direction );
     238        7472 :       for(unsigned i=0; i<natoms; ++i) {
     239       27440 :         for(unsigned j=0; j<3; ++j ) forces[j*natoms+i] += direction[i][j];
     240             :       }
     241             :       wasforced=true;
     242             :     }
     243         612 :     if( wasforced ) { unsigned ss=0; addForcesOnArguments( 0, forces, ss, getLabel() ); }
     244        4368 :   } else ActionWithVector::apply();
     245             : }
     246             : 
     247      180588 : void RMSDVector::performTask( const unsigned& current, MultiValue& myvals ) const {
     248      180588 :   unsigned natoms = getPntrToArgument(0)->getShape()[0] / 3;
     249             :   std::vector<Vector>& pos( myvals.getFirstAtomVector() );
     250             :   std::vector<std::vector<Vector> > & allder( myvals.getFirstAtomDerivativeVector() );
     251      180588 :   if( allder.size()!=2 ) allder.resize(2);
     252             :   std::vector<Vector>& der( allder[0] ); std::vector<Vector>& direction( allder[1] );
     253      180588 :   if( pos.size()!=natoms ) { pos.resize( natoms ); der.resize( natoms ); direction.resize( natoms ); }
     254     2528232 :   for(unsigned i=0; i<pos.size(); ++i) {
     255     9390576 :     for(unsigned j=0; j<3; ++j) pos[i][j] = getPntrToArgument(0)->get( j*natoms + i );
     256             :   }
     257      180588 :   double r = calculateRMSD( current, pos, der, direction );
     258      180588 :   unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream();
     259      180588 :   myvals.setValue( ostrn, r );
     260             : 
     261      180588 :   if( doNotCalculateDerivatives() ) return;
     262             : 
     263     1929648 :   for(unsigned i=0; i<natoms; i++) {
     264     7167264 :     for(unsigned j=0; j<3; ++j ) { myvals.addDerivative( ostrn, j*natoms+i, der[i][j] ); myvals.updateIndex( ostrn, j*natoms+i ); }
     265             :   }
     266             : }
     267             : 
     268      202902 : void RMSDVector::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
     269             :                                     const unsigned& bufstart, std::vector<double>& buffer ) const {
     270      202902 :   if( getConstPntrToComponent(valindex)->getRank()==1 ) { ActionWithVector::gatherStoredValue( valindex, code, myvals, bufstart, buffer ); return; }
     271             :   const std::vector<Vector>& direction( myvals.getConstFirstAtomDerivativeVector()[1] );
     272       42756 :   unsigned natoms = direction.size(); unsigned vindex = bufstart + 3*code*natoms;
     273      598584 :   for(unsigned i=0; i<natoms; ++i) {
     274     2223312 :     for(unsigned j=0; j<3; ++j ) buffer[vindex + j*natoms + i] += direction[i][j];
     275             :   }
     276             : }
     277             : 
     278       20442 : void RMSDVector::gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
     279       20442 :   if( myval->getRank()==1 ) { ActionWithVector::gatherForcesOnStoredValue( myval, itask, myvals, forces ); return; }
     280        1092 :   const std::vector<Vector>& direction( myvals.getConstFirstAtomDerivativeVector()[1] ); unsigned natoms = direction.size();
     281       15288 :   for(unsigned i=0; i<natoms; ++i) {
     282       56784 :     for(unsigned j=0; j<3; ++j ) forces[j*natoms+i] += direction[i][j];
     283             :   }
     284             : }
     285             : 
     286             : 
     287             : }
     288             : }
     289             : 
     290             : 
     291             :