Line data    Source code

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