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