Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2012-2020 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             :    See http://www.plumed.org for more information.
       5             :    This file is part of plumed, version 2.
       6             :    plumed is free software: you can redistribute it and/or modify
       7             :    it under the terms of the GNU Lesser General Public License as published by
       8             :    the Free Software Foundation, either version 3 of the License, or
       9             :    (at your option) any later version.
      10             :    plumed is distributed in the hope that it will be useful,
      11             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      12             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      13             :    GNU Lesser General Public License for more details.
      14             :    You should have received a copy of the GNU Lesser General Public License
      15             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      16             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      17             : #include "core/ActionRegister.h"
      18             : #include "core/PlumedMain.h"
      19             : #include "core/ActionSet.h"
      20             : #include "core/Group.h"
      21             : #include "core/ActionWithValue.h"
      22             : #include "core/ActionShortcut.h"
      23             : 
      24             : namespace PLMD {
      25             : namespace vatom {
      26             : 
      27             : class CenterShortcut : public ActionShortcut {
      28             : public:
      29             :   static void registerKeywords( Keywords& keys );
      30             :   explicit CenterShortcut(const ActionOptions&);
      31             : };
      32             : 
      33             : PLUMED_REGISTER_ACTION(CenterShortcut,"CENTER")
      34             : 
      35        8034 : void CenterShortcut::registerKeywords( Keywords& keys ) {
      36        8034 :   ActionShortcut::registerKeywords( keys );
      37        8034 :   keys.add("atoms","ATOMS","the group of atoms that appear in the definition of this center");
      38        8034 :   keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
      39        8034 :   keys.add("optional","WEIGHTS","what weights should be used when calculating the center.  If this keyword is not present the geometric center is computed. "
      40             :            "If WEIGHTS=@Masses is used the center of mass is computed.  If WEIGHTS=@charges the center of charge is computed.  If "
      41             :            "the label of an action is provided PLUMED assumes that that action calculates a list of symmetry functions that can be used "
      42             :            "as weights. Lastly, an explicit list of numbers to use as weights can be provided");
      43        8034 :   keys.addFlag("PHASES",false,"use trigonometric phases when computing position of center");
      44        8034 :   keys.addFlag("SAFE_PHASES",false,"use trignomentric phases when computing position of center but also compute the center in ths usual way and use this when the pbc are not set. "
      45             :                "There are two reasons for using this option (1) you are doing something that you know is really weird or (2) you are an idiot");
      46        8034 :   keys.addFlag("MASS",false,"calculate the center of mass");
      47        8034 :   keys.addActionNameSuffix("_FAST");
      48       16068 :   keys.setValueDescription("atom","the position of the center of mass");
      49        8034 :   keys.needsAction("MASS");
      50        8034 :   keys.needsAction("SUM");
      51        8034 :   keys.needsAction("CHARGE");
      52        8034 :   keys.needsAction("CONSTANT");
      53        8034 :   keys.needsAction("CUSTOM");
      54        8034 :   keys.needsAction("POSITION");
      55        8034 :   keys.needsAction("ARGS2VATOM");
      56        8034 : }
      57             : 
      58        7707 : CenterShortcut::CenterShortcut(const ActionOptions& ao):
      59             :   Action(ao),
      60        7707 :   ActionShortcut(ao) {
      61             :   // Read in what we are doing with the weights
      62             :   bool usemass;
      63        7707 :   parseFlag("MASS",usemass);
      64             :   bool safe_phases;
      65        7707 :   parseFlag("SAFE_PHASES",safe_phases);
      66        7707 :   std::vector<std::string> str_weights;
      67        7707 :   parseVector("WEIGHTS",str_weights);
      68        7707 :   if( usemass || str_weights.size()==0 || str_weights.size()>1 || (str_weights.size()==1 && str_weights[0]=="@Masses") ) {
      69        7698 :     std::string useph = "";
      70        7698 :     if( safe_phases ) {
      71             :       useph = " PHASES";
      72             :     }
      73        7698 :     if( usemass && str_weights.size()!=0 ) {
      74           3 :       error("WEIGHTS and MASS keywords cannot not be used simultaneously");
      75             :     }
      76             :     std::string wt_str;
      77        7697 :     if( str_weights.size()>0 ) {
      78        7044 :       wt_str="WEIGHTS=" + str_weights[0];
      79       28238 :       for(unsigned i=1; i<str_weights.size(); ++i) {
      80       42388 :         wt_str += "," + str_weights[i];
      81             :       }
      82             :     }
      83        7697 :     if( usemass || (str_weights.size()==1 && str_weights[0]=="@Masses") ) {
      84             :       wt_str = "MASS";
      85             :     }
      86       23092 :     readInputLine( getShortcutLabel() + ": CENTER_FAST " + wt_str + " " + convertInputLineToString() + useph );
      87             :     return;
      88             :   }
      89             :   // Read in the atoms
      90             :   std::string atlist;
      91           9 :   parse("ATOMS",atlist);
      92             :   // Calculate the mass of the vatom
      93           9 :   readInputLine( getShortcutLabel() + "_m: MASS ATOMS=" + atlist );
      94          18 :   readInputLine( getShortcutLabel() + "_mass: SUM PERIODIC=NO ARG=" + getShortcutLabel() + "_m" );
      95             :   // Calculate the charge of the vatom
      96           9 :   readInputLine( getShortcutLabel() + "_q: CHARGE ATOMS=" + atlist );
      97          18 :   readInputLine( getShortcutLabel() + "_charge: SUM PERIODIC=NO ARG=" + getShortcutLabel() + "_q" );
      98             :   // Retrieve the number of atoms
      99           9 :   ActionWithValue* am = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_m" );
     100           9 :   unsigned nat=am->copyOutput(0)->getNumberOfValues();
     101             :   // Get the weights to use for each atom
     102           9 :   std::string wlab = getShortcutLabel() + "_w";
     103           9 :   if( str_weights.size()>0 ) {
     104           9 :     if( str_weights.size()==1 ) {
     105           9 :       if( str_weights[0]=="@Charges" ) {
     106           0 :         wlab = getShortcutLabel() + "_q";
     107             :       } else {
     108             :         wlab=str_weights[0];
     109             :       }
     110           0 :     } else if( str_weights.size()==nat ) {
     111           0 :       std::string vals=str_weights[0];
     112           0 :       for(unsigned i=1; i<str_weights.size(); ++i) {
     113           0 :         vals += "," + str_weights[i];
     114             :       }
     115           0 :       readInputLine( getShortcutLabel() + "_w: CONSTANT VALUES=" + vals );
     116             :     } else {
     117           0 :       error("invalid input for WEIGHTS keyword " + str_weights[0] );
     118             :     }
     119             :   } else {
     120           0 :     std::string ones="1";
     121           0 :     for(unsigned i=1; i<nat; ++i) {
     122             :       ones += ",1";
     123             :     }
     124           0 :     readInputLine( getShortcutLabel() + "_w: CONSTANT VALUES=" + ones );
     125             :   }
     126             :   // Read in the instructions on how to compute the center of mass
     127             :   bool phases, nopbc;
     128           9 :   parseFlag("NOPBC",nopbc);
     129           9 :   if( safe_phases ) {
     130           0 :     phases=true;
     131             :   } else {
     132          18 :     parseFlag("PHASES",phases);
     133             :   }
     134             :   // This computes a center in the conventional way
     135           9 :   if( !phases || safe_phases ) {
     136             :     // Calculate the sum of the weights
     137           2 :     readInputLine( getShortcutLabel() + "_wnorm: SUM PERIODIC=NO ARG=" + wlab );
     138             :     // Compute the normalised weights
     139           6 :     readInputLine( getShortcutLabel() + "_weights: CUSTOM ARG=" + wlab + "," + getShortcutLabel() + "_wnorm FUNC=x/y PERIODIC=NO");
     140             :     // Get the positions into a multicolvar
     141           2 :     if( phases || nopbc ) {
     142           2 :       readInputLine( getShortcutLabel() + "_pos: POSITION NOPBC ATOMS=" + atlist );
     143             :     } else {
     144           2 :       readInputLine( getShortcutLabel() + "_pos: POSITION WHOLEMOLECULES ATOMS=" + atlist );
     145             :     }
     146             :     // Multiply each vector of positions by the weight
     147           6 :     readInputLine( getShortcutLabel() + "_xwvec: CUSTOM ARG=" + getShortcutLabel() + "_weights," + getShortcutLabel() + "_pos.x FUNC=x*y PERIODIC=NO");
     148           6 :     readInputLine( getShortcutLabel() + "_ywvec: CUSTOM ARG=" + getShortcutLabel() + "_weights," + getShortcutLabel() + "_pos.y FUNC=x*y PERIODIC=NO");
     149           6 :     readInputLine( getShortcutLabel() + "_zwvec: CUSTOM ARG=" + getShortcutLabel() + "_weights," + getShortcutLabel() + "_pos.z FUNC=x*y PERIODIC=NO");
     150             :     // And sum the weighted vectors
     151           4 :     readInputLine( getShortcutLabel() + "_x: SUM ARG=" + getShortcutLabel() + "_xwvec PERIODIC=NO");
     152           4 :     readInputLine( getShortcutLabel() + "_y: SUM ARG=" + getShortcutLabel() + "_ywvec PERIODIC=NO");
     153           6 :     readInputLine( getShortcutLabel() + "_z: SUM ARG=" + getShortcutLabel() + "_zwvec PERIODIC=NO");
     154             :   }
     155             :   // This computes a center using the trigonometric phases
     156           9 :   if( phases ) {
     157             :     // Get the positions into a multicolvar
     158           7 :     readInputLine( getShortcutLabel() + "_fpos: POSITION SCALED_COMPONENTS ATOMS=" + atlist );
     159             :     // Calculate the sines and cosines of the positions and multiply by the weights
     160          21 :     readInputLine( getShortcutLabel() + "_sina: CUSTOM ARG=" + wlab + "," + getShortcutLabel() + "_fpos.a FUNC=x*sin(2*pi*y) PERIODIC=NO");
     161          21 :     readInputLine( getShortcutLabel() + "_cosa: CUSTOM ARG=" + wlab + "," + getShortcutLabel() + "_fpos.a FUNC=x*cos(2*pi*y) PERIODIC=NO");
     162          21 :     readInputLine( getShortcutLabel() + "_sinb: CUSTOM ARG=" + wlab + "," + getShortcutLabel() + "_fpos.b FUNC=x*sin(2*pi*y) PERIODIC=NO");
     163          21 :     readInputLine( getShortcutLabel() + "_cosb: CUSTOM ARG=" + wlab + "," + getShortcutLabel() + "_fpos.b FUNC=x*cos(2*pi*y) PERIODIC=NO");
     164          21 :     readInputLine( getShortcutLabel() + "_sinc: CUSTOM ARG=" + wlab + "," + getShortcutLabel() + "_fpos.c FUNC=x*sin(2*pi*y) PERIODIC=NO");
     165          21 :     readInputLine( getShortcutLabel() + "_cosc: CUSTOM ARG=" + wlab + "," + getShortcutLabel() + "_fpos.c FUNC=x*cos(2*pi*y) PERIODIC=NO");
     166             :     // Sum the sines and cosines
     167          14 :     readInputLine( getShortcutLabel() + "_sinsuma: SUM ARG=" + getShortcutLabel() + "_sina PERIODIC=NO");
     168          14 :     readInputLine( getShortcutLabel() + "_cossuma: SUM ARG=" + getShortcutLabel() + "_cosa PERIODIC=NO");
     169          14 :     readInputLine( getShortcutLabel() + "_sinsumb: SUM ARG=" + getShortcutLabel() + "_sinb PERIODIC=NO");
     170          14 :     readInputLine( getShortcutLabel() + "_cossumb: SUM ARG=" + getShortcutLabel() + "_cosb PERIODIC=NO");
     171          14 :     readInputLine( getShortcutLabel() + "_sinsumc: SUM ARG=" + getShortcutLabel() + "_sinc PERIODIC=NO");
     172          14 :     readInputLine( getShortcutLabel() + "_cossumc: SUM ARG=" + getShortcutLabel() + "_cosc PERIODIC=NO");
     173             :     // And get the final position in fractional coordinates
     174          21 :     readInputLine( getShortcutLabel() + "_a: CUSTOM ARG=" + getShortcutLabel() + "_sinsuma," + getShortcutLabel() + "_cossuma FUNC=atan2(x,y)/(2*pi) PERIODIC=NO");
     175          21 :     readInputLine( getShortcutLabel() + "_b: CUSTOM ARG=" + getShortcutLabel() + "_sinsumb," + getShortcutLabel() + "_cossumb FUNC=atan2(x,y)/(2*pi) PERIODIC=NO");
     176          21 :     readInputLine( getShortcutLabel() + "_c: CUSTOM ARG=" + getShortcutLabel() + "_sinsumc," + getShortcutLabel() + "_cossumc FUNC=atan2(x,y)/(2*pi) PERIODIC=NO");
     177             :     // And create the virtual atom
     178           7 :     if( safe_phases ) {
     179           0 :       readInputLine( getShortcutLabel() + ": ARGS2VATOM XPOS=" + getShortcutLabel() + "_a YPOS=" + getShortcutLabel() + "_b ZPOS=" + getShortcutLabel() + "_c "
     180           0 :                      + " XBKP=" + getShortcutLabel() + "_x YBKP=" + getShortcutLabel() + "_y ZBKP=" + getShortcutLabel() + "_z "
     181           0 :                      + " MASS=" + getShortcutLabel() + "_mass CHARGE=" + getShortcutLabel() + "_charge FRACTIONAL");
     182             :     } else {
     183          28 :       readInputLine( getShortcutLabel() + ": ARGS2VATOM XPOS=" + getShortcutLabel() + "_a YPOS=" + getShortcutLabel() + "_b ZPOS=" + getShortcutLabel() + "_c "
     184          28 :                      + " MASS=" + getShortcutLabel() + "_mass CHARGE=" + getShortcutLabel() + "_charge FRACTIONAL");
     185             :     }
     186             :   } else {
     187             :     // And create the virtual atom
     188           8 :     readInputLine( getShortcutLabel() + ": ARGS2VATOM XPOS=" + getShortcutLabel() + "_x YPOS=" + getShortcutLabel() + "_y ZPOS=" + getShortcutLabel() + "_z "
     189           8 :                    + " MASS=" + getShortcutLabel() + "_mass CHARGE=" + getShortcutLabel() + "_charge ");
     190             :   }
     191        7711 : }
     192             : 
     193             : }
     194             : }