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 depracated.  In previous versions of PLUMED a more complex version of this method was implemented.
      33             : We felt, however, that the input syntax for the method was not very transparant.  We have thus provided this minimal action
      34             : that creates the input for calculating the MultiDomain RMSD for simple cases.  This action is a shortcut.  If you look at the log you can see how we
      35             : use the various actions that are in PLUMED to calculate the final quantity.  If you would like to implement some of the more
      36             : complicated CVs things that this could do with MULTI_RMSD looking at how this shortcut works will help you start.
      37             : 
      38             : 
      39             : \par Examples
      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           3 : void MultiRMSD::registerKeywords(Keywords& keys) {
      56           3 :   ActionShortcut::registerKeywords( keys );
      57           6 :   keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
      58           6 :   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           6 :   keys.addFlag("SQUARED",false," This should be set if you want the mean squared displacement instead of the root mean squared displacement");
      60           6 :   keys.addFlag("NOPBC",false,"don't use periodic boundary conditions");
      61           3 :   keys.setValueDescription("the sum of the multiple RMSD distances");
      62           9 :   keys.needsAction("CONSTANT"); keys.needsAction("WHOLEMOLECULES"); keys.needsAction("POSITION");
      63          12 :   keys.needsAction("CONCATENATE"); keys.needsAction("RMSD_VECTOR"); keys.needsAction("COMBINE"); keys.needsAction("CUSTOM");
      64           3 : }
      65             : 
      66           1 : MultiRMSD::MultiRMSD(const ActionOptions& ao):
      67             :   Action(ao),
      68           1 :   ActionShortcut(ao)
      69             : {
      70           2 :   warning("this action is depracated.  look at the log to see how it is implemented using the new syntax");
      71           2 :   std::string type; parse("TYPE",type); bool nopbc; parseFlag("NOPBC",nopbc);
      72           1 :   std::size_t dash=type.find_first_of("-");
      73           1 :   if( dash!=std::string::npos ) {
      74           2 :     if( type.substr(0,dash)=="MULTI" ) warning("MULTI is deprecated.  You can just use OPTIMAL/SIMPLE");
      75           0 :     else error("cannot understand type " + type );
      76           2 :     type = type.substr(dash+1);
      77             :   }
      78           2 :   std::string reference; parse("REFERENCE",reference); PDB pdb;
      79           1 :   if( !pdb.read(reference,usingNaturalUnits(),0.1/getUnits().getLength()) ) error("missing input file " + reference );
      80             : 
      81           1 :   unsigned nblocks =  pdb.getNumberOfAtomBlocks();
      82           1 :   if( nblocks<2 ) error("multidomain RMSD only has one block of atoms");
      83           1 :   std::string num; std::vector<unsigned> blocks( nblocks+1 ); blocks[0]=0;
      84           3 :   for(unsigned i=0; i<nblocks; ++i) blocks[i+1]=pdb.getAtomBlockEnds()[i];
      85             : 
      86           3 :   for(unsigned i=1; i<=nblocks; ++i) {
      87             :     // Setup a constant
      88           2 :     double asum=0; std::string bnum; Tools::convert( i, bnum );
      89          17 :     for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) asum += pdb.getOccupancy()[j];
      90           2 :     Vector center; center.zero();
      91          17 :     for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) center += ( pdb.getOccupancy()[j] / asum )*pdb.getPositions()[j];
      92             :     std::vector<double> vals;
      93           8 :     for(unsigned k=0; k<3; ++k) {
      94          51 :       for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) vals.push_back( pdb.getPositions()[j][k] - center[k] );
      95             :     }
      96           2 :     std::string valstr; Tools::convert( vals[0], valstr );
      97          45 :     for(unsigned i=1; i<vals.size(); ++i) { std::string rnum; Tools::convert( vals[i], rnum ); valstr += "," + rnum; }
      98             :     // Create the reference value
      99           4 :     readInputLine( getShortcutLabel() + "_ref" + bnum + ": CONSTANT VALUES=" + valstr );
     100             :     // Do whole molecules
     101           2 :     if( !nopbc ) {
     102           0 :       std::string num; Tools::convert( pdb.getAtomNumbers()[blocks[i-1]].serial(), num ); std::string wm_line = "WHOLEMOLECULES ENTITY0=" + num;
     103           0 :       for(unsigned j=blocks[i-1]+1; j<blocks[i]; ++j) { Tools::convert( pdb.getAtomNumbers()[j].serial(), num ); wm_line += "," + num; }
     104           0 :       readInputLine( wm_line );
     105             :     }
     106             :     // Get the positions of the atoms in this block
     107           4 :     std::string num; Tools::convert( pdb.getAtomNumbers()[blocks[i-1]].serial(), num ); std::string pos_line = getShortcutLabel() + "_cpos" + bnum + ": POSITION NOPBC ATOMS=" + num;
     108          15 :     for(unsigned j=blocks[i-1]+1; j<blocks[i]; ++j) { Tools::convert( pdb.getAtomNumbers()[j].serial(), num ); pos_line += "," + num; }
     109           2 :     readInputLine( pos_line );
     110             :     // Concatenate the positiosn together
     111           4 :     readInputLine( getShortcutLabel() + "_pos" + bnum + ": CONCATENATE ARG=" + getShortcutLabel() + "_cpos" + bnum + ".x," + getShortcutLabel() + "_cpos" + bnum + ".y," + getShortcutLabel() + "_cpos" + bnum + ".z");
     112             :     // Computer the RMSD for this block
     113           4 :     std::string rmsd_line = getShortcutLabel() + "_rmsd" + bnum + ": RMSD_VECTOR SQUARED ARG=" + getShortcutLabel() + "_pos" + bnum + "," + getShortcutLabel() + "_ref" + bnum;
     114             :     // Now align
     115           4 :     Tools::convert( pdb.getOccupancy()[blocks[i-1]], num ); rmsd_line += " ALIGN=" + num;
     116          15 :     for(unsigned j=blocks[i-1]+1; j<blocks[i]; ++j) { Tools::convert( pdb.getOccupancy()[j], num ); rmsd_line += "," + num; }
     117             :     // And displace
     118           4 :     Tools::convert( pdb.getBeta()[blocks[i-1]], num ); rmsd_line += " DISPLACE=" + num;
     119          15 :     for(unsigned j=blocks[i-1]+1; j<blocks[i]; ++j) { Tools::convert( pdb.getBeta()[j], num ); rmsd_line += "," + num; }
     120           4 :     readInputLine( rmsd_line + " TYPE=" + type );
     121             :   }
     122           3 :   std::string argstr = getShortcutLabel() + "_rmsd1"; for(unsigned i=1; i<nblocks; ++i) { std::string bnum; Tools::convert( i+1, bnum); argstr += "," + getShortcutLabel() + "_rmsd" + bnum; }
     123           1 :   bool squared; parseFlag("SQUARED",squared);
     124           1 :   if( !squared ) {
     125           2 :     readInputLine( getShortcutLabel() + "_2: COMBINE ARG=" + argstr + " PERIODIC=NO");
     126           2 :     readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_2 FUNC=sqrt(x) PERIODIC=NO");
     127           0 :   } else readInputLine( getShortcutLabel() + ": COMBINE ARG=" + argstr + " PERIODIC=NO");
     128           2 : }
     129             : 
     130             : }
     131             : }