Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-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 "Colvar.h"
      23             : #include "core/PlumedMain.h"
      24             : #include "core/ActionRegister.h"
      25             : #include "tools/RMSD.h"
      26             : #include "tools/PDB.h"
      27             : 
      28             : namespace PLMD {
      29             : namespace colvar {
      30             : 
      31             : class RMSD : public Colvar {
      32             : 
      33             :   bool squared;
      34             :   bool nopbc;
      35             :   PLMD::RMSD myrmsd;
      36             :   std::vector<Vector> der;
      37             : public:
      38             :   explicit RMSD(const ActionOptions&);
      39             :   void calculate() override;
      40             :   static void registerKeywords(Keywords& keys);
      41             : };
      42             : 
      43             : //+PLUMEDOC DCOLVAR RMSD_SCALAR
      44             : /*
      45             : Calculate the RMSD with respect to a reference structure.
      46             : 
      47             : \par Examples
      48             : 
      49             : */
      50             : //+ENDPLUMEDOC
      51             : 
      52             : //+PLUMEDOC DCOLVAR RMSD
      53             : /*
      54             : Calculate the RMSD with respect to a reference structure.
      55             : 
      56             : The aim with this colvar it to calculate something like:
      57             : 
      58             : \f[
      59             : d(X,X') = \vert X-X' \vert
      60             : \f]
      61             : 
      62             : where \f$ X \f$ is the instantaneous position of all the atoms in the system and
      63             : \f$ X' \f$ is the positions of the atoms in some reference structure provided as input.
      64             : \f$ d(X,X') \f$ thus measures the distance all the atoms have moved away from this reference configuration.
      65             : Oftentimes, it is only the internal motions of the structure - i.e. not the translations of the center of
      66             : mass or the rotations of the reference frame - that are interesting.  Hence, when calculating the
      67             : the root-mean-square deviation between the atoms in two configurations
      68             : you must first superimpose the two structures in some way. At present PLUMED provides two distinct ways
      69             : of performing this superposition.  The first method is applied when you use TYPE=SIMPLE in the input
      70             : line.  This instruction tells PLUMED that the root mean square deviation is to be calculated after the
      71             : positions of the geometric centers in the reference and instantaneous configurations are aligned.  In
      72             : other words \f$d(X,x')\f$ is to be calculated using:
      73             : 
      74             : \f[
      75             :  d(X,X') = \sqrt{ \sum_i \sum_\alpha^{x,y,z}  \frac{w_i}{\sum_j w_j}( X_{i,\alpha}-com_\alpha(X)-{X'}_{i,\alpha}+com_\alpha(X') )^2 }
      76             : \f]
      77             : with
      78             : \f[
      79             : com_\alpha(X)= \sum_i  \frac{w'_{i}}{\sum_j w'_j}X_{i,\alpha}
      80             : \f]
      81             : and
      82             : \f[
      83             : com_\alpha(X')= \sum_i  \frac{w'_{i}}{\sum_j w'_j}X'_{i,\alpha}
      84             : \f]
      85             : Obviously, \f$ com_\alpha(X) \f$ and  \f$ com_\alpha(X') \f$  represent the positions of the center of mass in the reference
      86             : and instantaneous configurations if the weights $w'$ are set equal to the atomic masses.  If the weights are all set equal to
      87             : one, however, \f$com_\alpha(X) \f$ and  \f$ com_\alpha(X') \f$ are the positions of the geometric centers.
      88             : Notice that there are sets of weights:  \f$ w' \f$ and  \f$ w \f$. The first is used to calculate the position of the center of mass
      89             : (so it determines how the atoms are \e aligned).  Meanwhile, the second is used when calculating how far the atoms have actually been
      90             : \e displaced.  These weights are assigned in the reference configuration that you provide as input (i.e. the appear in the input file
      91             : to this action that you set using REFERENCE=whatever.pdb).  This input reference configuration consists of a simple pdb file
      92             : containing the set of atoms for which you want to calculate the RMSD displacement and their positions in the reference configuration.
      93             : It is important to note that the indices in this pdb need to be set correctly.  The indices in this file determine the indices of the
      94             : instantaneous atomic positions that are used by PLUMED when calculating this colvar.  As such if you want to calculate the RMSD distance
      95             : moved by the first, fourth, sixth and twenty eighth atoms in the MD codes input file then the indices of the corresponding reference positions in this pdb
      96             : file should be set equal to 1, 4, 6 and 28.
      97             : 
      98             : The pdb input file should also contain the values of \f$w\f$ and \f$w'\f$. In particular, the OCCUPANCY column (the first column after the coordinates)
      99             : is used provides the values of \f$ w'\f$ that are used to calculate the position of the center of mass.  The BETA column (the second column
     100             : after the Cartesian coordinates) is used to provide the \f$ w \f$ values which are used in the the calculation of the displacement.
     101             : Please note that it is possible to use fractional values for beta and for the occupancy. However, we recommend you only do this when
     102             : you really know what you are doing however as the results can be rather strange.
     103             : 
     104             : In PDB files the atomic coordinates and box lengths should be in Angstroms unless
     105             : you are working with natural units.  If you are working with natural units then the coordinates
     106             : should be in your natural length unit.  For more details on the PDB file format visit http://www.wwpdb.org/docs.html.
     107             : Make sure your PDB file is correctly formatted as explained \ref pdbreader "in this page".
     108             : 
     109             : A different method is used to calculate the RMSD distance when you use TYPE=OPTIMAL on the input line.  In this case  the root mean square
     110             : deviation is calculated after the positions of geometric centers in the reference and instantaneous configurations are aligned AND after
     111             : an optimal alignment of the two frames is performed so that motion due to rotation of the reference frame between the two structures is
     112             : removed.  The equation for \f$d(X,X')\f$ in this case reads:
     113             : 
     114             : \f[
     115             : d(X,X') = \sqrt{ \sum_i \sum_\alpha^{x,y,z}  \frac{w_i}{\sum_j w_j}[ X_{i,\alpha}-com_\alpha(X)- \sum_\beta M(X,X',w')_{\alpha,\beta}({X'}_{i,\beta}-com_\beta(X')) ]^2 }
     116             : \f]
     117             : 
     118             : where \f$ M(X,X',w') \f$ is the optimal alignment matrix which is calculated using the Kearsley \cite kearsley algorithm.  Again different sets of
     119             : weights are used for the alignment (\f$w'\f$) and for the displacement calculations (\f$w\f$).
     120             : This gives a great deal of flexibility as it allows you to use a different sets of atoms (which may or may not overlap) for the alignment and displacement
     121             : parts of the calculation. This may be very useful when you want to calculate how a ligand moves about in a protein cavity as you can use the protein as a reference
     122             : system and do no alignment of the ligand.
     123             : 
     124             : (Note: when this form of RMSD is used to calculate the secondary structure variables (\ref ALPHARMSD, \ref ANTIBETARMSD and \ref PARABETARMSD
     125             : all the atoms in the segment are assumed to be part of both the alignment and displacement sets and all weights are set equal to one)
     126             : 
     127             : Please note that there are a number of other methods for calculating the distance between the instantaneous configuration and a reference configuration
     128             : that are available in plumed.  More information on these various methods can be found in the section of the manual on \ref dists.
     129             : 
     130             : When running with periodic boundary conditions, the atoms should be
     131             : in the proper periodic image. This is done automatically since PLUMED 2.5,
     132             : by considering the ordered list of atoms and rebuilding molecules using a procedure
     133             : that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
     134             : rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
     135             : which actually modifies the coordinates stored in PLUMED.
     136             : 
     137             : In case you want to recover the old behavior you should use the NOPBC flag.
     138             : In that case you need to take care that atoms are in the correct
     139             : periodic image.
     140             : 
     141             : \par Examples
     142             : 
     143             : The following tells plumed to calculate the RMSD distance between
     144             : the positions of the atoms in the reference file and their instantaneous
     145             : position.  The Kearsley algorithm is used so this is done optimally.
     146             : 
     147             : \plumedfile
     148             : RMSD REFERENCE=file.pdb TYPE=OPTIMAL
     149             : \endplumedfile
     150             : 
     151             : The reference configuration is specified in a pdb file that will have a format similar to the one shown below:
     152             : 
     153             : \auxfile{file.pdb}
     154             : ATOM      1  CL  ALA     1      -3.171   0.295   2.045  1.00  1.00
     155             : ATOM      5  CLP ALA     1      -1.819  -0.143   1.679  1.00  1.00
     156             : ATOM      6  OL  ALA     1      -1.177  -0.889   2.401  1.00  1.00
     157             : ATOM      7  NL  ALA     1      -1.313   0.341   0.529  1.00  1.00
     158             : ATOM      8  HL  ALA     1      -1.845   0.961  -0.011  1.00  1.00
     159             : END
     160             : \endauxfile
     161             : 
     162             : ...
     163             : 
     164             : */
     165             : //+ENDPLUMEDOC
     166             : 
     167             : PLUMED_REGISTER_ACTION(RMSD,"RMSD_SCALAR")
     168             : 
     169         163 : void RMSD::registerKeywords(Keywords& keys) {
     170         163 :   Colvar::registerKeywords(keys); keys.setDisplayName("RMSD");
     171         326 :   keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
     172         326 :   keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed.  Should be OPTIMAL or SIMPLE.");
     173         326 :   keys.addFlag("SQUARED",false," This should be set if you want mean squared displacement instead of RMSD ");
     174         163 :   keys.setValueDescription("the RMSD between the instantaneous structure and the reference structure that was input");
     175         163 : }
     176             : 
     177          81 : RMSD::RMSD(const ActionOptions&ao):
     178             :   PLUMED_COLVAR_INIT(ao),
     179          81 :   squared(false),
     180          81 :   nopbc(false)
     181             : {
     182             :   std::string reference;
     183         163 :   parse("REFERENCE",reference);
     184             :   std::string type;
     185          81 :   type.assign("SIMPLE");
     186          81 :   parse("TYPE",type);
     187          81 :   parseFlag("SQUARED",squared);
     188          81 :   parseFlag("NOPBC",nopbc);
     189          81 :   checkRead();
     190             : 
     191         163 :   addValueWithDerivatives(); setNotPeriodic();
     192          81 :   PDB pdb;
     193             : 
     194             :   // read everything in ang and transform to nm if we are not in natural units
     195          81 :   if( !pdb.read(reference,usingNaturalUnits(),0.1/getUnits().getLength()) )
     196           0 :     error("missing input file " + reference );
     197          81 :   myrmsd.set( pdb, type, true, true );
     198             : 
     199          81 :   std::vector<AtomNumber> atoms( pdb.getAtomNumbers() );
     200          81 :   requestAtoms( atoms ); der.resize( atoms.size() );
     201             : 
     202          80 :   log.printf("  reference from file %s\n",reference.c_str());
     203          80 :   log.printf("  which contains %d atoms\n",getNumberOfAtoms());
     204          80 :   log.printf("  with indices : ");
     205        3867 :   for(unsigned i=0; i<atoms.size(); ++i) {
     206        3787 :     if(i%25==0) log<<"\n";
     207        3787 :     log.printf("%d ",atoms[i].serial());
     208             :   }
     209          80 :   log.printf("\n");
     210          80 :   log.printf("  method for alignment : %s \n",type.c_str() );
     211          80 :   if(squared)log.printf("  chosen to use SQUARED option for MSD instead of RMSD\n");
     212          80 :   if(nopbc) log.printf("  without periodic boundary conditions\n");
     213          19 :   else      log.printf("  using periodic boundary conditions\n");
     214         164 : }
     215             : 
     216             : 
     217             : // calculator
     218       40516 : void RMSD::calculate() {
     219       40516 :   if(!nopbc) makeWhole();
     220       40516 :   double r=myrmsd.calculate( getPositions(), der, squared );
     221             : 
     222       40516 :   setValue(r);
     223    14617897 :   for(unsigned i=0; i<getNumberOfAtoms(); i++) setAtomsDerivatives( i, der[i] );
     224       40516 :   setBoxDerivativesNoPbc();
     225       40516 : }
     226             : 
     227             : }
     228             : }
     229             : 
     230             : 
     231             :