Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2013-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 "MetricRegister.h"
      23             : #include "RMSDBase.h"
      24             : #include "tools/Matrix.h"
      25             : #include "tools/RMSD.h"
      26             : 
      27             : namespace PLMD {
      28             : 
      29             : class OptimalRMSD : public RMSDBase {
      30             : private:
      31             :   bool fast;
      32             :   RMSD myrmsd;
      33             : public:
      34             :   explicit OptimalRMSD(const ReferenceConfigurationOptions& ro);
      35             :   void read( const PDB& ) override;
      36             :   double calc( const std::vector<Vector>& pos, ReferenceValuePack& myder, const bool& squared ) const override;
      37          45 :   bool pcaIsEnabledForThisReference() override { return true; }
      38         910 :   void setupRMSDObject() override { myrmsd.clear(); myrmsd.set(getAlign(),getDisplace(),getReferencePositions(),"OPTIMAL"); }
      39          72 :   void setupPCAStorage( ReferenceValuePack& mypack ) override {
      40             :     mypack.switchOnPCAOption();
      41          72 :     mypack.centeredpos.resize( getNumberOfAtoms() );
      42          72 :     mypack.displacement.resize( getNumberOfAtoms() );
      43          72 :     mypack.DRotDPos.resize(3,3); mypack.rot.resize(1);
      44          72 :   }
      45             :   void extractAtomicDisplacement( const std::vector<Vector>& pos, std::vector<Vector>& direction ) const override;
      46             :   double projectAtomicDisplacementOnVector( const bool& normalized, const std::vector<Vector>& vecs, ReferenceValuePack& mypack ) const override;
      47             : };
      48             : 
      49       11293 : PLUMED_REGISTER_METRIC(OptimalRMSD,"OPTIMAL")
      50             : 
      51         437 : OptimalRMSD::OptimalRMSD(const ReferenceConfigurationOptions& ro ):
      52             :   ReferenceConfiguration(ro),
      53         437 :   RMSDBase(ro)
      54             : {
      55         437 :   fast=ro.usingFastOption();
      56         437 : }
      57             : 
      58         427 : void OptimalRMSD::read( const PDB& pdb ) {
      59         427 :   readReference( pdb ); setupRMSDObject();
      60         427 : }
      61             : 
      62      218933 : double OptimalRMSD::calc( const std::vector<Vector>& pos, ReferenceValuePack& myder, const bool& squared ) const {
      63             :   double d;
      64      218933 :   if( myder.calcUsingPCAOption() ) {
      65        2786 :     std::vector<Vector> centeredreference( getNumberOfAtoms () );
      66        2786 :     d=myrmsd.calc_PCAelements(pos,myder.getAtomVector(),myder.rot[0],myder.DRotDPos,myder.getAtomsDisplacementVector(),myder.centeredpos,centeredreference,squared);
      67        2786 :     unsigned nat = pos.size();
      68       48548 :     for(unsigned i=0; i<nat; ++i) { myder.getAtomsDisplacementVector()[i] -= getReferencePosition(i); myder.getAtomsDisplacementVector()[i] *= getDisplace()[i]; }
      69      216147 :   } else if( fast ) {
      70      192640 :     if( getAlign()==getDisplace() ) d=myrmsd.optimalAlignment<false,true>(getAlign(),getDisplace(),pos,getReferencePositions(),myder.getAtomVector(),squared);
      71           2 :     else d=myrmsd.optimalAlignment<false,false>(getAlign(),getDisplace(),pos,getReferencePositions(),myder.getAtomVector(),squared);
      72             :   } else {
      73       23507 :     if( getAlign()==getDisplace() ) d=myrmsd.optimalAlignment<true,true>(getAlign(),getDisplace(),pos,getReferencePositions(),myder.getAtomVector(),squared);
      74         111 :     else d=myrmsd.optimalAlignment<true,false>(getAlign(),getDisplace(),pos,getReferencePositions(),myder.getAtomVector(),squared);
      75             :   }
      76    17762561 :   myder.clear(); for(unsigned i=0; i<pos.size(); ++i) myder.setAtomDerivatives( i, myder.getAtomVector()[i] );
      77      218933 :   if( !myder.updateComplete() ) myder.updateDynamicLists();
      78      218933 :   return d;
      79             : }
      80             : 
      81        1107 : void OptimalRMSD::extractAtomicDisplacement( const std::vector<Vector>& pos, std::vector<Vector>& direction ) const {
      82        1107 :   std::vector<Tensor> rot(1);  Matrix<std::vector<Vector> > DRotDPos( 3, 3 );
      83        1107 :   std::vector<Vector> centeredreference( getNumberOfAtoms() ), centeredpos( getNumberOfAtoms() ), avector( getNumberOfAtoms() );
      84        1107 :   myrmsd.calc_PCAelements(pos,avector,rot[0],DRotDPos,direction,centeredpos,centeredreference,true);
      85       15498 :   unsigned nat = pos.size(); for(unsigned i=0; i<nat; ++i) direction[i] = getDisplace()[i]*( direction[i] - getReferencePosition(i) );
      86        1107 : }
      87             : 
      88        1158 : double OptimalRMSD::projectAtomicDisplacementOnVector( const bool& normalized, const std::vector<Vector>& vecs, ReferenceValuePack& mypack ) const {
      89             :   plumed_dbg_assert( mypack.calcUsingPCAOption() );
      90             : 
      91        1158 :   double proj=0.0; mypack.clear();
      92       15662 :   for(unsigned i=0; i<vecs.size(); ++i) {
      93       14504 :     proj += dotProduct( mypack.getAtomsDisplacementVector()[i], vecs[i] );
      94             :   }
      95        4632 :   for(unsigned a=0; a<3; a++) {
      96       13896 :     for(unsigned b=0; b<3; b++) {
      97      140958 :       double tmp1=0.; for(unsigned n=0; n<getNumberOfAtoms(); n++) tmp1+=mypack.centeredpos[n][b]*vecs[n][a];
      98             : 
      99      140958 :       for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
     100      130536 :         if( normalized ) mypack.addAtomDerivatives( iat, getDisplace()[iat]*mypack.DRotDPos[a][b][iat]*tmp1 );
     101      127764 :         else mypack.addAtomDerivatives( iat, mypack.DRotDPos[a][b][iat]*tmp1 );
     102             :       }
     103             :     }
     104             :   }
     105        1158 :   Tensor trot=mypack.rot[0].transpose();
     106        1158 :   Vector v1; v1.zero(); double prefactor = 1. / static_cast<double>( getNumberOfAtoms() );
     107       15662 :   for(unsigned n=0; n<getNumberOfAtoms(); n++) v1+=prefactor*matmul(trot,vecs[n]);
     108        1158 :   if( normalized ) {
     109         374 :     for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) mypack.addAtomDerivatives( iat, getDisplace()[iat]*(matmul(trot,vecs[iat]) - v1) );
     110             :   } else {
     111       15288 :     for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) mypack.addAtomDerivatives( iat, (matmul(trot,vecs[iat]) - v1) );
     112             :   }
     113        1158 :   if( !mypack.updateComplete() ) mypack.updateDynamicLists();
     114             : 
     115        1158 :   return proj;
     116             : }
     117             : 
     118             : }