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