Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2013-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 "MetricRegister.h"
      23             : #include "RMSDBase.h"
      24             : #include "tools/Matrix.h"
      25             : #include "tools/RMSD.h"
      26             : 
      27             : namespace PLMD {
      28             : 
      29        1836 : 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& );
      36             :   double calc( const std::vector<Vector>& pos, ReferenceValuePack& myder, const bool& squared ) const ;
      37           3 :   bool pcaIsEnabledForThisReference() { return true; }
      38         923 :   void setupRMSDObject() { myrmsd.clear(); myrmsd.set(getAlign(),getDisplace(),getReferencePositions(),"OPTIMAL"); }
      39          68 :   void setupPCAStorage( ReferenceValuePack& mypack ) {
      40          68 :     mypack.switchOnPCAOption();
      41          68 :     mypack.centeredpos.resize( getNumberOfAtoms() );
      42          68 :     mypack.displacement.resize( getNumberOfAtoms() );
      43          68 :     mypack.DRotDPos.resize(3,3); mypack.rot.resize(1);
      44          68 :   }
      45             :   double projectAtomicDisplacementOnVector( const unsigned& iv, const Matrix<Vector>& vecs, const std::vector<Vector>& pos, ReferenceValuePack& mypack ) const ;
      46             : };
      47             : 
      48        3441 : PLUMED_REGISTER_METRIC(OptimalRMSD,"OPTIMAL")
      49             : 
      50         918 : OptimalRMSD::OptimalRMSD(const ReferenceConfigurationOptions& ro ):
      51             :   ReferenceConfiguration(ro),
      52         918 :   RMSDBase(ro)
      53             : {
      54         918 :   fast=ro.usingFastOption();
      55         918 : }
      56             : 
      57         344 : void OptimalRMSD::read( const PDB& pdb ) {
      58         344 :   readReference( pdb ); setupRMSDObject();
      59         344 : }
      60             : 
      61      191601 : double OptimalRMSD::calc( const std::vector<Vector>& pos, ReferenceValuePack& myder, const bool& squared ) const {
      62             :   double d;
      63      191601 :   if( myder.calcUsingPCAOption() ) {
      64        1124 :     std::vector<Vector> centeredreference( getNumberOfAtoms () );
      65        1124 :     d=myrmsd.calc_PCAelements(pos,myder.getAtomVector(),myder.rot[0],myder.DRotDPos,myder.getAtomsDisplacementVector(),myder.centeredpos,centeredreference,squared);
      66        1124 :     unsigned nat = pos.size(); for(unsigned i=0; i<nat; ++i) myder.getAtomsDisplacementVector()[i] -= getReferencePosition(i);
      67      190871 :   } else if( fast ) {
      68      167397 :     if( getAlign()==getDisplace() ) d=myrmsd.optimalAlignment<false,true>(getAlign(),getDisplace(),pos,getReferencePositions(),myder.getAtomVector(),squared);
      69           2 :     else d=myrmsd.optimalAlignment<false,false>(getAlign(),getDisplace(),pos,getReferencePositions(),myder.getAtomVector(),squared);
      70             :   } else {
      71       23474 :     if( getAlign()==getDisplace() ) d=myrmsd.optimalAlignment<true,true>(getAlign(),getDisplace(),pos,getReferencePositions(),myder.getAtomVector(),squared);
      72         111 :     else d=myrmsd.optimalAlignment<true,false>(getAlign(),getDisplace(),pos,getReferencePositions(),myder.getAtomVector(),squared);
      73             :   }
      74      193314 :   myder.clear(); for(unsigned i=0; i<pos.size(); ++i) myder.setAtomDerivatives( i, myder.getAtomVector()[i] );
      75      193735 :   if( !myder.updateComplete() ) myder.updateDynamicLists();
      76      193763 :   return d;
      77             : }
      78             : 
      79          66 : double OptimalRMSD::projectAtomicDisplacementOnVector( const unsigned& iv, const Matrix<Vector>& vecs, const std::vector<Vector>& pos, ReferenceValuePack& mypack ) const {
      80             :   plumed_dbg_assert( mypack.calcUsingPCAOption() );
      81             : 
      82          66 :   double proj=0.0; mypack.clear();
      83         374 :   for(unsigned i=0; i<pos.size(); ++i) {
      84         308 :     proj += dotProduct( mypack.getAtomsDisplacementVector()[i], vecs(iv,i) );
      85             :   }
      86         264 :   for(unsigned a=0; a<3; a++) {
      87         792 :     for(unsigned b=0; b<3; b++) {
      88        3366 :       for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
      89        2772 :         double tmp1=0.;
      90        2772 :         for(unsigned n=0; n<getNumberOfAtoms(); n++) tmp1+=mypack.centeredpos[n][b]*vecs(iv,n)[a];
      91        2772 :         mypack.addAtomDerivatives( iat, mypack.DRotDPos[a][b][iat]*tmp1 );
      92             :       }
      93             :     }
      94             :   }
      95          66 :   Tensor trot=mypack.rot[0].transpose();
      96          66 :   Vector v1; v1.zero(); double prefactor = 1. / static_cast<double>( getNumberOfAtoms() );
      97          66 :   for(unsigned n=0; n<getNumberOfAtoms(); n++) v1+=prefactor*matmul(trot,vecs(iv,n));
      98          66 :   for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) mypack.addAtomDerivatives( iat, matmul(trot,vecs(iv,iat))-v1 );
      99          66 :   if( !mypack.updateComplete() ) mypack.updateDynamicLists();
     100             : 
     101          66 :   return proj;
     102             : }
     103             : 
     104        2523 : }