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 {
      38          45 :     return true;
      39             :   }
      40         473 :   void setupRMSDObject() override {
      41         473 :     myrmsd.clear();
      42         473 :     myrmsd.set(getAlign(),getDisplace(),getReferencePositions(),"OPTIMAL");
      43         473 :   }
      44          72 :   void setupPCAStorage( ReferenceValuePack& mypack ) override {
      45             :     mypack.switchOnPCAOption();
      46          72 :     mypack.centeredpos.resize( getNumberOfAtoms() );
      47          72 :     mypack.displacement.resize( getNumberOfAtoms() );
      48             :     mypack.DRotDPos.resize(3,3);
      49          72 :     mypack.rot.resize(1);
      50          72 :   }
      51             :   void extractAtomicDisplacement( const std::vector<Vector>& pos, std::vector<Vector>& direction ) const override;
      52             :   double projectAtomicDisplacementOnVector( const bool& normalized, const std::vector<Vector>& vecs, ReferenceValuePack& mypack ) const override;
      53             : };
      54             : 
      55       14695 : PLUMED_REGISTER_METRIC(OptimalRMSD,"OPTIMAL")
      56             : 
      57         455 : OptimalRMSD::OptimalRMSD(const ReferenceConfigurationOptions& ro ):
      58             :   ReferenceConfiguration(ro),
      59         455 :   RMSDBase(ro) {
      60         455 :   fast=ro.usingFastOption();
      61         455 : }
      62             : 
      63         427 : void OptimalRMSD::read( const PDB& pdb ) {
      64         427 :   readReference( pdb );
      65         427 :   setupRMSDObject();
      66         427 : }
      67             : 
      68      273521 : double OptimalRMSD::calc( const std::vector<Vector>& pos, ReferenceValuePack& myder, const bool& squared ) const {
      69             :   double d;
      70      273521 :   if( myder.calcUsingPCAOption() ) {
      71        2786 :     std::vector<Vector> centeredreference( getNumberOfAtoms () );
      72        2786 :     d=myrmsd.calc_PCAelements(pos,myder.getAtomVector(),myder.rot[0],myder.DRotDPos,myder.getAtomsDisplacementVector(),myder.centeredpos,centeredreference,squared);
      73        2786 :     unsigned nat = pos.size();
      74       48548 :     for(unsigned i=0; i<nat; ++i) {
      75       45762 :       myder.getAtomsDisplacementVector()[i] -= getReferencePosition(i);
      76       45762 :       myder.getAtomsDisplacementVector()[i] *= getDisplace()[i];
      77             :     }
      78      270735 :   } else if( fast ) {
      79      219934 :     if( getAlign()==getDisplace() ) {
      80      219932 :       d=myrmsd.optimalAlignment<false,true>(getAlign(),getDisplace(),pos,getReferencePositions(),myder.getAtomVector(),squared);
      81             :     } else {
      82           2 :       d=myrmsd.optimalAlignment<false,false>(getAlign(),getDisplace(),pos,getReferencePositions(),myder.getAtomVector(),squared);
      83             :     }
      84             :   } else {
      85       50801 :     if( getAlign()==getDisplace() ) {
      86       50690 :       d=myrmsd.optimalAlignment<true,true>(getAlign(),getDisplace(),pos,getReferencePositions(),myder.getAtomVector(),squared);
      87             :     } else {
      88         111 :       d=myrmsd.optimalAlignment<true,false>(getAlign(),getDisplace(),pos,getReferencePositions(),myder.getAtomVector(),squared);
      89             :     }
      90             :   }
      91      273521 :   myder.clear();
      92    19454789 :   for(unsigned i=0; i<pos.size(); ++i) {
      93    19181268 :     myder.setAtomDerivatives( i, myder.getAtomVector()[i] );
      94             :   }
      95      273521 :   if( !myder.updateComplete() ) {
      96      273521 :     myder.updateDynamicLists();
      97             :   }
      98      273521 :   return d;
      99             : }
     100             : 
     101        1107 : void OptimalRMSD::extractAtomicDisplacement( const std::vector<Vector>& pos, std::vector<Vector>& direction ) const {
     102        1107 :   std::vector<Tensor> rot(1);
     103             :   Matrix<std::vector<Vector> > DRotDPos( 3, 3 );
     104        1107 :   std::vector<Vector> centeredreference( getNumberOfAtoms() ), centeredpos( getNumberOfAtoms() ), avector( getNumberOfAtoms() );
     105        1107 :   myrmsd.calc_PCAelements(pos,avector,rot[0],DRotDPos,direction,centeredpos,centeredreference,true);
     106        1107 :   unsigned nat = pos.size();
     107       15498 :   for(unsigned i=0; i<nat; ++i) {
     108       14391 :     direction[i] = getDisplace()[i]*( direction[i] - getReferencePosition(i) );
     109             :   }
     110        1107 : }
     111             : 
     112        1158 : double OptimalRMSD::projectAtomicDisplacementOnVector( const bool& normalized, const std::vector<Vector>& vecs, ReferenceValuePack& mypack ) const {
     113             :   plumed_dbg_assert( mypack.calcUsingPCAOption() );
     114             : 
     115             :   double proj=0.0;
     116        1158 :   mypack.clear();
     117       15662 :   for(unsigned i=0; i<vecs.size(); ++i) {
     118       14504 :     proj += dotProduct( mypack.getAtomsDisplacementVector()[i], vecs[i] );
     119             :   }
     120        4632 :   for(unsigned a=0; a<3; a++) {
     121       13896 :     for(unsigned b=0; b<3; b++) {
     122             :       double tmp1=0.;
     123      140958 :       for(unsigned n=0; n<getNumberOfAtoms(); n++) {
     124      130536 :         tmp1+=mypack.centeredpos[n][b]*vecs[n][a];
     125             :       }
     126             : 
     127      140958 :       for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
     128      130536 :         if( normalized ) {
     129        2772 :           mypack.addAtomDerivatives( iat, getDisplace()[iat]*mypack.DRotDPos[a][b][iat]*tmp1 );
     130             :         } else {
     131      127764 :           mypack.addAtomDerivatives( iat, mypack.DRotDPos[a][b][iat]*tmp1 );
     132             :         }
     133             :       }
     134             :     }
     135             :   }
     136        1158 :   Tensor trot=mypack.rot[0].transpose();
     137        1158 :   Vector v1;
     138        1158 :   v1.zero();
     139        1158 :   double prefactor = 1. / static_cast<double>( getNumberOfAtoms() );
     140       15662 :   for(unsigned n=0; n<getNumberOfAtoms(); n++) {
     141       14504 :     v1+=prefactor*matmul(trot,vecs[n]);
     142             :   }
     143        1158 :   if( normalized ) {
     144         374 :     for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
     145         308 :       mypack.addAtomDerivatives( iat, getDisplace()[iat]*(matmul(trot,vecs[iat]) - v1) );
     146             :     }
     147             :   } else {
     148       15288 :     for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
     149       14196 :       mypack.addAtomDerivatives( iat, (matmul(trot,vecs[iat]) - v1) );
     150             :     }
     151             :   }
     152        1158 :   if( !mypack.updateComplete() ) {
     153        1158 :     mypack.updateDynamicLists();
     154             :   }
     155             : 
     156        1158 :   return proj;
     157             : }
     158             : 
     159             : }