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 "Colvar.h"
      23             : #include "core/PlumedMain.h"
      24             : #include "ActionRegister.h"
      25             : #include "tools/PDB.h"
      26             : #include "reference/MultiDomainRMSD.h"
      27             : #include "reference/MetricRegister.h"
      28             : #include "core/Atoms.h"
      29             : 
      30             : namespace PLMD {
      31             : namespace colvar {
      32             : 
      33             : class MultiRMSD : public Colvar {
      34             : 
      35             :   std::unique_ptr<PLMD::MultiDomainRMSD> rmsd;
      36             :   bool squared;
      37             :   MultiValue myvals;
      38             :   ReferenceValuePack mypack;
      39             :   bool nopbc;
      40             : 
      41             : public:
      42             :   explicit MultiRMSD(const ActionOptions&);
      43             :   void calculate() override;
      44             :   static void registerKeywords(Keywords& keys);
      45             : };
      46             : 
      47             : //+PLUMEDOC DCOLVAR MULTI_RMSD
      48             : /*
      49             : Calculate the RMSD distance moved by a number of separated domains from their positions in a reference structure.
      50             : 
      51             : 
      52             : When you have large proteins the calculation of the root mean squared deviation between all the atoms in a reference
      53             : structure and the instantaneous configuration becomes prohibitively expensive.  You may thus instead want to calculate
      54             : the RMSD between the atoms in a set of domains of your protein and your reference structure.  That is to say:
      55             : 
      56             : \f[
      57             : d(X,X_r) = \sqrt{ \sum_{i} w_i\vert X_i - X_i' \vert^2 }
      58             : \f]
      59             : 
      60             : where here the sum is over the domains of the protein, \f$X_i\f$ represents the positions of the atoms in domain \f$i\f$
      61             : in the instantaneous configuration and \f$X_i'\f$ is the positions of the atoms in domain \f$i\f$ in the reference
      62             : configuration.  \f$w_i\f$ is an optional weight.
      63             : 
      64             : The distances for each of the domains in the above sum can be calculated using the \ref DRMSD or \ref RMSD measures or
      65             : using a combination of these distance.  The reference configuration is specified in a pdb file like the one below:
      66             : 
      67             : \auxfile{file1.pdb}
      68             : ATOM      2  O   ALA     2      -0.926  -2.447  -0.497  1.00  1.00      DIA  O
      69             : ATOM      4  HNT ALA     2       0.533  -0.396   1.184  1.00  1.00      DIA  H
      70             : ATOM      6  HT1 ALA     2      -0.216  -2.590   1.371  1.00  1.00      DIA  H
      71             : ATOM      7  HT2 ALA     2      -0.309  -1.255   2.315  1.00  1.00      DIA  H
      72             : ATOM      8  HT3 ALA     2      -1.480  -1.560   1.212  1.00  1.00      DIA  H
      73             : ATOM      9  CAY ALA     2      -0.096   2.144  -0.669  1.00  1.00      DIA  C
      74             : ATOM     10  HY1 ALA     2       0.871   2.385  -0.588  1.00  1.00      DIA  H
      75             : TER
      76             : ATOM     12  HY3 ALA     2      -0.520   2.679  -1.400  1.00  1.00      DIA  H
      77             : ATOM     14  OY  ALA     2      -1.139   0.931  -0.973  1.00  1.00      DIA  O
      78             : ATOM     16  HN  ALA     2       1.713   1.021  -0.873  1.00  1.00      DIA  H
      79             : ATOM     18  HA  ALA     2       0.099  -0.774  -2.218  1.00  1.00      DIA  H
      80             : ATOM     19  CB  ALA     2       2.063  -1.223  -1.276  1.00  1.00      DIA  C
      81             : ATOM     20  HB1 ALA     2       2.670  -0.716  -2.057  1.00  1.00      DIA  H
      82             : ATOM     21  HB2 ALA     2       2.556  -1.051  -0.295  1.00  1.00      DIA  H
      83             : ATOM     22  HB3 ALA     2       2.070  -2.314  -1.490  1.00  1.00      DIA  H
      84             : END
      85             : \endauxfile
      86             : 
      87             : with the TER keyword being used to separate the various domains in you protein.
      88             : 
      89             : 
      90             : \par Examples
      91             : 
      92             : The following tells plumed to calculate the RMSD distance between
      93             : the positions of the atoms in the reference file and their instantaneous
      94             : position.  The Kearsley algorithm for each of the domains.
      95             : 
      96             : \plumedfile
      97             : MULTI_RMSD REFERENCE=file1.pdb TYPE=MULTI-OPTIMAL
      98             : \endplumedfile
      99             : 
     100             : The following tells plumed to calculate the RMSD distance between the positions of
     101             : the atoms in the domains of reference the reference structure and their instantaneous
     102             : positions.  Here distances are calculated using the \ref DRMSD measure.
     103             : 
     104             : \plumedfile
     105             : MULTI_RMSD REFERENCE=file1.pdb TYPE=MULTI-DRMSD
     106             : \endplumedfile
     107             : 
     108             : in this case it is possible to use the following DRMSD options in the pdb file using the REMARK syntax:
     109             : \verbatim
     110             : NOPBC to calculate distances without PBC
     111             : LOWER_CUTOFF=# only pairs of atoms further than LOWER_CUTOFF are considered in the calculation
     112             : UPPER_CUTOFF=# only pairs of atoms further than UPPER_CUTOFF are considered in the calculation
     113             : \endverbatim
     114             : as shown in the following example
     115             : 
     116             : \auxfile{file2.pdb}
     117             : REMARK NOPBC
     118             : REMARK LOWER_CUTOFF=0.1
     119             : REMARK UPPER_CUTOFF=0.8
     120             : ATOM      2  O   ALA     2      -0.926  -2.447  -0.497  1.00  1.00      DIA  O
     121             : ATOM      4  HNT ALA     2       0.533  -0.396   1.184  1.00  1.00      DIA  H
     122             : ATOM      6  HT1 ALA     2      -0.216  -2.590   1.371  1.00  1.00      DIA  H
     123             : ATOM      7  HT2 ALA     2      -0.309  -1.255   2.315  1.00  1.00      DIA  H
     124             : ATOM      8  HT3 ALA     2      -1.480  -1.560   1.212  1.00  1.00      DIA  H
     125             : ATOM      9  CAY ALA     2      -0.096   2.144  -0.669  1.00  1.00      DIA  C
     126             : ATOM     10  HY1 ALA     2       0.871   2.385  -0.588  1.00  1.00      DIA  H
     127             : TER
     128             : ATOM     12  HY3 ALA     2      -0.520   2.679  -1.400  1.00  1.00      DIA  H
     129             : ATOM     14  OY  ALA     2      -1.139   0.931  -0.973  1.00  1.00      DIA  O
     130             : ATOM     16  HN  ALA     2       1.713   1.021  -0.873  1.00  1.00      DIA  H
     131             : ATOM     18  HA  ALA     2       0.099  -0.774  -2.218  1.00  1.00      DIA  H
     132             : ATOM     19  CB  ALA     2       2.063  -1.223  -1.276  1.00  1.00      DIA  C
     133             : ATOM     20  HB1 ALA     2       2.670  -0.716  -2.057  1.00  1.00      DIA  H
     134             : ATOM     21  HB2 ALA     2       2.556  -1.051  -0.295  1.00  1.00      DIA  H
     135             : ATOM     22  HB3 ALA     2       2.070  -2.314  -1.490  1.00  1.00      DIA  H
     136             : END
     137             : \endauxfile
     138             : 
     139             : 
     140             : */
     141             : //+ENDPLUMEDOC
     142             : 
     143       10425 : PLUMED_REGISTER_ACTION(MultiRMSD,"MULTI_RMSD")
     144             : 
     145           4 : void MultiRMSD::registerKeywords(Keywords& keys) {
     146           4 :   Colvar::registerKeywords(keys);
     147           8 :   keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
     148           8 :   keys.add("compulsory","TYPE","MULTI-SIMPLE","the manner in which RMSD alignment is performed.  Should be MULTI-OPTIMAL, MULTI-OPTIMAL-FAST,  MULTI-SIMPLE or MULTI-DRMSD.");
     149           8 :   keys.addFlag("SQUARED",false," This should be set if you want the mean squared displacement instead of the root mean squared displacement");
     150           4 : }
     151             : 
     152           3 : MultiRMSD::MultiRMSD(const ActionOptions&ao):
     153           3 :   PLUMED_COLVAR_INIT(ao),squared(false),myvals(1,0), mypack(0,0,myvals),nopbc(false)
     154             : {
     155             :   std::string reference;
     156           6 :   parse("REFERENCE",reference);
     157             :   std::string type;
     158           3 :   type.assign("SIMPLE");
     159           3 :   parse("TYPE",type);
     160           3 :   parseFlag("SQUARED",squared);
     161           3 :   parseFlag("NOPBC",nopbc);
     162           3 :   checkRead();
     163             : 
     164           3 :   addValueWithDerivatives(); setNotPeriodic();
     165           3 :   PDB pdb;
     166             : 
     167             :   // read everything in ang and transform to nm if we are not in natural units
     168           6 :   if( !pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength()) )
     169           0 :     error("missing input file " + reference );
     170             : 
     171           6 :   rmsd=metricRegister().create<MultiDomainRMSD>(type,pdb);
     172             :   // Do not align molecule if we are doing DRMSD for domains and NOPBC has been specified in input
     173           6 :   if( pdb.hasFlag("NOPBC") ) nopbc=true;
     174             : 
     175             :   std::vector<AtomNumber> atoms;
     176           3 :   rmsd->getAtomRequests( atoms );
     177           3 :   requestAtoms( atoms );
     178             : 
     179           3 :   myvals.resize( 1, 3*atoms.size()+9 ); mypack.resize( 0, atoms.size() );
     180          48 :   for(unsigned i=0; i<atoms.size(); ++i) mypack.setAtomIndex( i, i );
     181             : 
     182           3 :   log.printf("  reference from file %s\n",reference.c_str());
     183           3 :   log.printf("  which contains %d atoms\n",getNumberOfAtoms());
     184           3 :   log.printf("  with indices : ");
     185          48 :   for(unsigned i=0; i<atoms.size(); ++i) {
     186          45 :     if(i%25==0) log<<"\n";
     187          45 :     log.printf("%d ",atoms[i].serial());
     188             :   }
     189           3 :   log.printf("\n");
     190           3 :   log.printf("  method for alignment : %s \n",type.c_str() );
     191           3 :   if(squared)log.printf("  chosen to use SQUARED option for MSD instead of RMSD\n");
     192           6 : }
     193             : 
     194             : // calculator
     195          15 : void MultiRMSD::calculate() {
     196          15 :   if(!nopbc) makeWhole();
     197          15 :   double r=rmsd->calculate( getPositions(), getPbc(), mypack, squared );
     198             : 
     199          15 :   setValue(r);
     200         240 :   for(unsigned i=0; i<getNumberOfAtoms(); i++) setAtomsDerivatives( i, mypack.getAtomDerivative(i) );
     201             : 
     202          15 :   if( !mypack.virialWasSet() ) setBoxDerivativesNoPbc();
     203           5 :   else setBoxDerivatives( mypack.getBoxDerivatives() );
     204          15 : }
     205             : 
     206             : }
     207             : }
     208             : 
     209             : 
     210             :