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 "core/ActionShortcut.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "core/ActionWithArguments.h"
      25             : #include "tools/PDB.h"
      26             : #include "core/PlumedMain.h"
      27             : 
      28             : //+PLUMEDOC DCOLVAR MULTI_RMSD
      29             : /*
      30             : Calculate RMSD distances for different domains and combine them.
      31             : 
      32             : This action is largely deprecated.  In previous versions of PLUMED a more complex version of this method was implemented.
      33             : You can see an example in [example from the plumed nest](https://www.plumed-nest.org/eggs/20/026/).
      34             : 
      35             : We felt that the input syntax for the method was not very transparant.  We have thus provided this minimal action
      36             : that creates the input for calculating the MultiDomain RMSD for simple cases.  This action is a shortcut.  If you look at the inputs in the
      37             : egg that is linked above you can see how we
      38             : use the various actions that are in PLUMED to calculate the final quantity.  If you would like to implement some
      39             : complicated CVs that are linear combinations of multiple RMSD calculations looking at how this shortcut works will help you start.
      40             : 
      41             : */
      42             : //+ENDPLUMEDOC
      43             : 
      44             : namespace PLMD {
      45             : namespace colvar {
      46             : 
      47             : class MultiRMSD : public ActionShortcut {
      48             : public:
      49             :   static void registerKeywords(Keywords& keys);
      50             :   explicit MultiRMSD(const ActionOptions&);
      51             : };
      52             : 
      53             : PLUMED_REGISTER_ACTION(MultiRMSD,"MULTI_RMSD")
      54             : 
      55           5 : void MultiRMSD::registerKeywords(Keywords& keys) {
      56           5 :   ActionShortcut::registerKeywords( keys );
      57           5 :   keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
      58           5 :   keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed.  Should be MULTI-OPTIMAL, MULTI-OPTIMAL-FAST,  MULTI-SIMPLE or MULTI-DRMSD.");
      59           5 :   keys.addFlag("SQUARED",false," This should be set if you want the mean squared displacement instead of the root mean squared displacement");
      60           5 :   keys.addFlag("NOPBC",false,"don't use periodic boundary conditions");
      61          10 :   keys.setValueDescription("scalar","the sum of the multiple RMSD distances");
      62           5 :   keys.setDeprecated("RMSD");
      63           5 :   keys.needsAction("CONSTANT");
      64           5 :   keys.needsAction("WHOLEMOLECULES");
      65           5 :   keys.needsAction("POSITION");
      66           5 :   keys.needsAction("CONCATENATE");
      67           5 :   keys.needsAction("RMSD_VECTOR");
      68           5 :   keys.needsAction("COMBINE");
      69           5 :   keys.needsAction("CUSTOM");
      70           5 : }
      71             : 
      72           1 : MultiRMSD::MultiRMSD(const ActionOptions& ao):
      73             :   Action(ao),
      74           1 :   ActionShortcut(ao) {
      75           2 :   warning("this action is deprecated.  look at the log to see how it is implemented using the new syntax");
      76             :   std::string type;
      77           1 :   parse("TYPE",type);
      78             :   bool nopbc;
      79           1 :   parseFlag("NOPBC",nopbc);
      80           1 :   std::size_t dash=type.find_first_of("-");
      81           1 :   if( dash!=std::string::npos ) {
      82           1 :     if( type.substr(0,dash)=="MULTI" ) {
      83           1 :       warning("MULTI is deprecated.  You can just use OPTIMAL/SIMPLE");
      84             :     } else {
      85           0 :       error("cannot understand type " + type );
      86             :     }
      87           2 :     type = type.substr(dash+1);
      88             :   }
      89             :   std::string reference;
      90           1 :   parse("REFERENCE",reference);
      91           1 :   PDB pdb;
      92           1 :   if( !pdb.read(reference,usingNaturalUnits(),0.1/getUnits().getLength()) ) {
      93           0 :     error("missing input file " + reference );
      94             :   }
      95             : 
      96           1 :   unsigned nblocks =  pdb.getNumberOfAtomBlocks();
      97           1 :   if( nblocks<2 ) {
      98           0 :     error("multidomain RMSD only has one block of atoms");
      99             :   }
     100             :   std::string num;
     101           1 :   std::vector<unsigned> blocks( nblocks+1 );
     102           1 :   blocks[0]=0;
     103           3 :   for(unsigned i=0; i<nblocks; ++i) {
     104           2 :     blocks[i+1]=pdb.getAtomBlockEnds()[i];
     105             :   }
     106             : 
     107           3 :   for(unsigned i=1; i<=nblocks; ++i) {
     108             :     // Setup a constant
     109             :     double asum=0;
     110             :     std::string bnum;
     111           2 :     Tools::convert( i, bnum );
     112          17 :     for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) {
     113          15 :       asum += pdb.getOccupancy()[j];
     114             :     }
     115           2 :     Vector center;
     116             :     center.zero();
     117          17 :     for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) {
     118          15 :       center += ( pdb.getOccupancy()[j] / asum )*pdb.getPositions()[j];
     119             :     }
     120           2 :     std::vector<double> vals;
     121           8 :     for(unsigned k=0; k<3; ++k) {
     122          51 :       for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) {
     123          45 :         vals.push_back( pdb.getPositions()[j][k] - center[k] );
     124             :       }
     125             :     }
     126             :     std::string valstr;
     127           2 :     Tools::convert( vals[0], valstr );
     128          45 :     for(unsigned ii=1; ii<vals.size(); ++ii) {
     129             :       std::string rnum;
     130          43 :       Tools::convert( vals[ii], rnum );
     131          43 :       valstr += "," + rnum;
     132             :     }
     133             :     // Create the reference value
     134           4 :     readInputLine( getShortcutLabel() + "_ref" + bnum + ": CONSTANT VALUES=" + valstr );
     135             :     // Do whole molecules
     136           2 :     if( !nopbc ) {
     137           0 :       Tools::convert( pdb.getAtomNumbers()[blocks[i-1]].serial(), num );
     138           0 :       std::string wm_line = "WHOLEMOLECULES ENTITY0=" + num;
     139           0 :       for(unsigned j=blocks[i-1]+1; j<blocks[i]; ++j) {
     140           0 :         Tools::convert( pdb.getAtomNumbers()[j].serial(), num );
     141           0 :         wm_line += "," + num;
     142             :       }
     143           0 :       readInputLine( wm_line );
     144             :     }
     145             :     // Get the positions of the atoms in this block
     146           2 :     Tools::convert( pdb.getAtomNumbers()[blocks[i-1]].serial(), num );
     147           6 :     std::string pos_line = getShortcutLabel() + "_cpos" + bnum + ": POSITION NOPBC ATOMS=" + num;
     148          15 :     for(unsigned j=blocks[i-1]+1; j<blocks[i]; ++j) {
     149          13 :       Tools::convert( pdb.getAtomNumbers()[j].serial(), num );
     150          26 :       pos_line += "," + num;
     151             :     }
     152           2 :     readInputLine( pos_line );
     153             :     // Concatenate the positiosn together
     154          16 :     readInputLine( getShortcutLabel() + "_pos" + bnum + ": CONCATENATE ARG=" + getShortcutLabel() + "_cpos" + bnum + ".x," + getShortcutLabel() + "_cpos" + bnum + ".y," + getShortcutLabel() + "_cpos" + bnum + ".z");
     155             :     // Computer the RMSD for this block
     156          10 :     std::string rmsd_line = getShortcutLabel() + "_rmsd" + bnum + ": RMSD_VECTOR SQUARED ARG=" + getShortcutLabel() + "_pos" + bnum + "," + getShortcutLabel() + "_ref" + bnum;
     157             :     // Now align
     158           2 :     Tools::convert( pdb.getOccupancy()[blocks[i-1]], num );
     159           4 :     rmsd_line += " ALIGN=" + num;
     160          15 :     for(unsigned j=blocks[i-1]+1; j<blocks[i]; ++j) {
     161          13 :       Tools::convert( pdb.getOccupancy()[j], num );
     162          26 :       rmsd_line += "," + num;
     163             :     }
     164             :     // And displace
     165           2 :     Tools::convert( pdb.getBeta()[blocks[i-1]], num );
     166           4 :     rmsd_line += " DISPLACE=" + num;
     167          15 :     for(unsigned j=blocks[i-1]+1; j<blocks[i]; ++j) {
     168          13 :       Tools::convert( pdb.getBeta()[j], num );
     169          26 :       rmsd_line += "," + num;
     170             :     }
     171           4 :     readInputLine( rmsd_line + " TYPE=" + type );
     172           2 :   }
     173           1 :   std::string argstr = getShortcutLabel() + "_rmsd1";
     174           2 :   for(unsigned i=1; i<nblocks; ++i) {
     175             :     std::string bnum;
     176           1 :     Tools::convert( i+1, bnum);
     177           3 :     argstr += "," + getShortcutLabel() + "_rmsd" + bnum;
     178             :   }
     179             :   bool squared;
     180           1 :   parseFlag("SQUARED",squared);
     181           1 :   if( !squared ) {
     182           2 :     readInputLine( getShortcutLabel() + "_2: COMBINE ARG=" + argstr + " PERIODIC=NO");
     183           3 :     readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_2 FUNC=sqrt(x) PERIODIC=NO");
     184             :   } else {
     185           0 :     readInputLine( getShortcutLabel() + ": COMBINE ARG=" + argstr + " PERIODIC=NO");
     186             :   }
     187          46 : }
     188             : 
     189             : }
     190             : }