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             : 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 of the more
      39             : complicated CVs things that this could do with MULTI_RMSD 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           3 : void MultiRMSD::registerKeywords(Keywords& keys) {
      56           3 :   ActionShortcut::registerKeywords( keys );
      57           3 :   keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
      58           3 :   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           3 :   keys.addFlag("SQUARED",false," This should be set if you want the mean squared displacement instead of the root mean squared displacement");
      60           3 :   keys.addFlag("NOPBC",false,"don't use periodic boundary conditions");
      61           6 :   keys.setValueDescription("scalar","the sum of the multiple RMSD distances");
      62           3 :   keys.needsAction("CONSTANT");
      63           3 :   keys.needsAction("WHOLEMOLECULES");
      64           3 :   keys.needsAction("POSITION");
      65           3 :   keys.needsAction("CONCATENATE");
      66           3 :   keys.needsAction("RMSD_VECTOR");
      67           3 :   keys.needsAction("COMBINE");
      68           3 :   keys.needsAction("CUSTOM");
      69           3 : }
      70             : 
      71           1 : MultiRMSD::MultiRMSD(const ActionOptions& ao):
      72             :   Action(ao),
      73           1 :   ActionShortcut(ao) {
      74           2 :   warning("this action is depracated.  look at the log to see how it is implemented using the new syntax");
      75             :   std::string type;
      76           1 :   parse("TYPE",type);
      77             :   bool nopbc;
      78           1 :   parseFlag("NOPBC",nopbc);
      79           1 :   std::size_t dash=type.find_first_of("-");
      80           1 :   if( dash!=std::string::npos ) {
      81           2 :     if( type.substr(0,dash)=="MULTI" ) {
      82           1 :       warning("MULTI is deprecated.  You can just use OPTIMAL/SIMPLE");
      83             :     } else {
      84           0 :       error("cannot understand type " + type );
      85             :     }
      86           2 :     type = type.substr(dash+1);
      87             :   }
      88             :   std::string reference;
      89           1 :   parse("REFERENCE",reference);
      90           1 :   PDB pdb;
      91           1 :   if( !pdb.read(reference,usingNaturalUnits(),0.1/getUnits().getLength()) ) {
      92           0 :     error("missing input file " + reference );
      93             :   }
      94             : 
      95           1 :   unsigned nblocks =  pdb.getNumberOfAtomBlocks();
      96           1 :   if( nblocks<2 ) {
      97           0 :     error("multidomain RMSD only has one block of atoms");
      98             :   }
      99             :   std::string num;
     100           1 :   std::vector<unsigned> blocks( nblocks+1 );
     101           1 :   blocks[0]=0;
     102           3 :   for(unsigned i=0; i<nblocks; ++i) {
     103           2 :     blocks[i+1]=pdb.getAtomBlockEnds()[i];
     104             :   }
     105             : 
     106           3 :   for(unsigned i=1; i<=nblocks; ++i) {
     107             :     // Setup a constant
     108             :     double asum=0;
     109             :     std::string bnum;
     110           2 :     Tools::convert( i, bnum );
     111          17 :     for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) {
     112          15 :       asum += pdb.getOccupancy()[j];
     113             :     }
     114           2 :     Vector center;
     115           2 :     center.zero();
     116          17 :     for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) {
     117          15 :       center += ( pdb.getOccupancy()[j] / asum )*pdb.getPositions()[j];
     118             :     }
     119             :     std::vector<double> vals;
     120           8 :     for(unsigned k=0; k<3; ++k) {
     121          51 :       for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) {
     122          45 :         vals.push_back( pdb.getPositions()[j][k] - center[k] );
     123             :       }
     124             :     }
     125             :     std::string valstr;
     126           2 :     Tools::convert( vals[0], valstr );
     127          45 :     for(unsigned i=1; i<vals.size(); ++i) {
     128             :       std::string rnum;
     129          43 :       Tools::convert( vals[i], rnum );
     130          86 :       valstr += "," + rnum;
     131             :     }
     132             :     // Create the reference value
     133           4 :     readInputLine( getShortcutLabel() + "_ref" + bnum + ": CONSTANT VALUES=" + valstr );
     134             :     // Do whole molecules
     135           2 :     if( !nopbc ) {
     136             :       std::string num;
     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             :     std::string num;
     147           2 :     Tools::convert( pdb.getAtomNumbers()[blocks[i-1]].serial(), num );
     148           4 :     std::string pos_line = getShortcutLabel() + "_cpos" + bnum + ": POSITION NOPBC ATOMS=" + num;
     149          15 :     for(unsigned j=blocks[i-1]+1; j<blocks[i]; ++j) {
     150          13 :       Tools::convert( pdb.getAtomNumbers()[j].serial(), num );
     151          26 :       pos_line += "," + num;
     152             :     }
     153           2 :     readInputLine( pos_line );
     154             :     // Concatenate the positiosn together
     155           4 :     readInputLine( getShortcutLabel() + "_pos" + bnum + ": CONCATENATE ARG=" + getShortcutLabel() + "_cpos" + bnum + ".x," + getShortcutLabel() + "_cpos" + bnum + ".y," + getShortcutLabel() + "_cpos" + bnum + ".z");
     156             :     // Computer the RMSD for this block
     157           4 :     std::string rmsd_line = getShortcutLabel() + "_rmsd" + bnum + ": RMSD_VECTOR SQUARED ARG=" + getShortcutLabel() + "_pos" + bnum + "," + getShortcutLabel() + "_ref" + bnum;
     158             :     // Now align
     159           2 :     Tools::convert( pdb.getOccupancy()[blocks[i-1]], num );
     160           4 :     rmsd_line += " ALIGN=" + num;
     161          15 :     for(unsigned j=blocks[i-1]+1; j<blocks[i]; ++j) {
     162          13 :       Tools::convert( pdb.getOccupancy()[j], num );
     163          26 :       rmsd_line += "," + num;
     164             :     }
     165             :     // And displace
     166           2 :     Tools::convert( pdb.getBeta()[blocks[i-1]], num );
     167           4 :     rmsd_line += " DISPLACE=" + num;
     168          15 :     for(unsigned j=blocks[i-1]+1; j<blocks[i]; ++j) {
     169          13 :       Tools::convert( pdb.getBeta()[j], num );
     170          26 :       rmsd_line += "," + num;
     171             :     }
     172           4 :     readInputLine( rmsd_line + " TYPE=" + type );
     173             :   }
     174           1 :   std::string argstr = getShortcutLabel() + "_rmsd1";
     175           2 :   for(unsigned i=1; i<nblocks; ++i) {
     176             :     std::string bnum;
     177           1 :     Tools::convert( i+1, bnum);
     178           2 :     argstr += "," + getShortcutLabel() + "_rmsd" + bnum;
     179             :   }
     180             :   bool squared;
     181           1 :   parseFlag("SQUARED",squared);
     182           1 :   if( !squared ) {
     183           2 :     readInputLine( getShortcutLabel() + "_2: COMBINE ARG=" + argstr + " PERIODIC=NO");
     184           2 :     readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_2 FUNC=sqrt(x) PERIODIC=NO");
     185             :   } else {
     186           0 :     readInputLine( getShortcutLabel() + ": COMBINE ARG=" + argstr + " PERIODIC=NO");
     187             :   }
     188           2 : }
     189             : 
     190             : }
     191             : }