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 "DRMSD.h"
      23             : #include "MetricRegister.h"
      24             : #include "tools/Pbc.h"
      25             : 
      26             : namespace PLMD {
      27             : 
      28       13855 : PLUMED_REGISTER_METRIC(DRMSD,"DRMSD")
      29             : 
      30          37 : DRMSD::DRMSD( const ReferenceConfigurationOptions& ro ):
      31             :   ReferenceConfiguration( ro ),
      32             :   SingleDomainRMSD( ro ),
      33          37 :   nopbc(true),
      34          37 :   bounds_were_set(false),
      35          37 :   lower(0),
      36          37 :   upper(std::numeric_limits<double>::max( )) {
      37          37 : }
      38             : 
      39          37 : void DRMSD::setBoundsOnDistances( bool dopbc, double lbound, double ubound ) {
      40          37 :   bounds_were_set=true;
      41          37 :   nopbc=!dopbc;
      42          37 :   lower=lbound;
      43          37 :   upper=ubound;
      44          37 : }
      45             : 
      46          11 : void DRMSD::readBounds( const PDB& pdb ) {
      47          11 :   if( bounds_were_set ) {
      48          11 :     return;
      49             :   }
      50             :   double tmp;
      51           0 :   nopbc=pdb.hasFlag("NOPBC");
      52           0 :   if( pdb.getArgumentValue("LOWER_CUTOFF",tmp) ) {
      53           0 :     lower=tmp;
      54             :   }
      55           0 :   if( pdb.getArgumentValue("UPPER_CUTOFF",tmp) ) {
      56           0 :     upper=tmp;
      57             :   }
      58           0 :   bounds_were_set=true;
      59             : }
      60             : 
      61           9 : void DRMSD::read( const PDB& pdb ) {
      62           9 :   readAtomsFromPDB( pdb );
      63           9 :   readBounds( pdb );
      64           9 :   setup_targets();
      65           9 : }
      66             : 
      67          26 : void DRMSD::setReferenceAtoms( const std::vector<Vector>& conf, const std::vector<double>& align_in, const std::vector<double>& displace_in ) {
      68          26 :   SingleDomainRMSD::setReferenceAtoms( conf, align_in, displace_in );
      69          26 :   setup_targets();
      70          26 : }
      71             : 
      72          35 : void DRMSD::setup_targets() {
      73          35 :   plumed_massert( bounds_were_set, "I am missing a call to DRMSD::setBoundsOnDistances");
      74             : 
      75          35 :   unsigned natoms = getNumberOfReferencePositions();
      76         839 :   for(unsigned i=0; i<natoms-1; ++i) {
      77       11955 :     for(unsigned j=i+1; j<natoms; ++j) {
      78       11151 :       double distance = delta( getReferencePosition(i), getReferencePosition(j) ).modulo();
      79       11151 :       if(distance < upper && distance > lower ) {
      80       10283 :         targets[std::make_pair(i,j)] = distance;
      81             :       }
      82             :     }
      83             :   }
      84          35 :   if( targets.empty() ) {
      85           0 :     error("drmsd will compare no distances - check upper and lower bounds are sensible");
      86             :   }
      87          35 : }
      88             : 
      89       38582 : double DRMSD::calc( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceValuePack& myder, const bool& squared ) const {
      90             :   plumed_dbg_assert(!targets.empty());
      91             : 
      92       38582 :   Vector distance;
      93       38582 :   myder.clear();
      94             :   double drmsd=0.;
      95    15553320 :   for(const auto & it : targets) {
      96             : 
      97    15514738 :     const unsigned i=getAtomIndex( it.first.first );
      98    15514738 :     const unsigned j=getAtomIndex( it.first.second );
      99             : 
     100    15514738 :     if(nopbc) {
     101       29400 :       distance=delta( pos[i], pos[j] );
     102             :     } else {
     103    15485338 :       distance=pbc.distance( pos[i], pos[j] );
     104             :     }
     105             : 
     106    15514738 :     const double len = distance.modulo();
     107    15514738 :     const double diff = len - it.second;
     108    15514738 :     const double der = diff / len;
     109             : 
     110    15514738 :     drmsd += diff * diff;
     111    15514738 :     myder.addAtomDerivatives( i, -der * distance );
     112    15514738 :     myder.addAtomDerivatives( j,  der * distance );
     113    15514738 :     myder.addBoxDerivatives( - der * Tensor(distance,distance) );
     114             :   }
     115             : 
     116       38582 :   const double inpairs = 1./static_cast<double>(targets.size());
     117             :   double idrmsd;
     118             : 
     119       38582 :   if(squared) {
     120          10 :     drmsd = drmsd * inpairs;
     121          10 :     idrmsd = 2.0 * inpairs;
     122             :   } else {
     123       38572 :     drmsd = std::sqrt( drmsd * inpairs );
     124       38572 :     idrmsd = inpairs / drmsd ;
     125             :   }
     126             : 
     127       38582 :   myder.scaleAllDerivatives( idrmsd );
     128             : 
     129       38582 :   return drmsd;
     130             : }
     131             : 
     132             : }