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        8049 : void CenterShortcut::registerKeywords( Keywords& keys ) {
      36        8049 :   ActionShortcut::registerKeywords( keys );
      37       16098 :   keys.add("atoms","ATOMS","the group of atoms that you are calculating the Gyration Tensor for");
      38       16098 :   keys.add("compulsory","TYPE","RADIUS","The type of calculation relative to the Gyration Tensor you want to perform");
      39       16098 :   keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
      40       16098 :   keys.add("optional","WEIGHTS","what weights should be used when calculating the center.  If this keyword is not present the geometric center is computed. "
      41             :            "If WEIGHTS=@Masses is used the center of mass is computed.  If WEIGHTS=@charges the center of charge is computed.  If "
      42             :            "the label of an action is provided PLUMED assumes that that action calculates a list of symmetry functions that can be used "
      43             :            "as weights. Lastly, an explicit list of numbers to use as weights can be provided");
      44       16098 :   keys.addFlag("PHASES",false,"use trigonometric phases when computing position of center");
      45       16098 :   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. "
      46             :                "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");
      47       16098 :   keys.addFlag("MASS",false,"calculate the center of mass"); keys.addActionNameSuffix("_FAST");
      48        8049 :   keys.setValueDescription("the position of the center of mass");
      49       32196 :   keys.needsAction("MASS"); keys.needsAction("SUM"); keys.needsAction("CHARGE"); keys.needsAction("CONSTANT");
      50       24147 :   keys.needsAction("CUSTOM"); keys.needsAction("POSITION"); keys.needsAction("ARGS2VATOM");
      51        8049 : }
      52             : 
      53        7748 : CenterShortcut::CenterShortcut(const ActionOptions& ao):
      54             :   Action(ao),
      55        7748 :   ActionShortcut(ao)
      56             : {
      57             :   // Read in what we are doing with the weights
      58       15496 :   bool usemass; parseFlag("MASS",usemass);
      59        7748 :   std::vector<std::string> str_weights; parseVector("WEIGHTS",str_weights);
      60        7757 :   if( usemass || str_weights.size()==0 || str_weights.size()>1 || (str_weights.size()==1 && str_weights[0]=="@Masses") ) {
      61        7739 :     if( usemass && str_weights.size()!=0 ) error("WEIGHTS and MASS keywords cannot not be used simultaneously");
      62             :     std::string wt_str;
      63        7738 :     if( str_weights.size()>0 ) {
      64       35278 :       wt_str="WEIGHTS=" + str_weights[0]; for(unsigned i=1; i<str_weights.size(); ++i) wt_str += "," + str_weights[i];
      65             :     }
      66        7738 :     if( usemass || (str_weights.size()==1 && str_weights[0]=="@Masses") ) wt_str = "MASS";
      67       15480 :     readInputLine( getShortcutLabel() + ": CENTER_FAST " + wt_str + " " + convertInputLineToString() );
      68             :     return;
      69             :   }
      70             :   // Read in the atoms
      71           9 :   std::string atlist; parse("ATOMS",atlist);
      72             :   // Calculate the mass of the vatom
      73          18 :   readInputLine( getShortcutLabel() + "_m: MASS ATOMS=" + atlist );
      74          18 :   readInputLine( getShortcutLabel() + "_mass: SUM PERIODIC=NO ARG=" + getShortcutLabel() + "_m" );
      75             :   // Calculate the charge of the vatom
      76          18 :   readInputLine( getShortcutLabel() + "_q: CHARGE ATOMS=" + atlist );
      77          18 :   readInputLine( getShortcutLabel() + "_charge: SUM PERIODIC=NO ARG=" + getShortcutLabel() + "_q" );
      78             :   // Retrieve the number of atoms
      79           9 :   ActionWithValue* am = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_m" );
      80           9 :   unsigned nat=am->copyOutput(0)->getNumberOfValues();
      81             :   // Get the weights to use for each atom
      82           9 :   std::string wlab = getShortcutLabel() + "_w";
      83           9 :   if( str_weights.size()>0 ) {
      84           9 :     if( str_weights.size()==1 ) {
      85           9 :       if( str_weights[0]=="@Charges" ) wlab = getShortcutLabel() + "_q";
      86             :       else wlab=str_weights[0];
      87           0 :     } else if( str_weights.size()==nat ) {
      88           0 :       std::string vals=str_weights[0]; for(unsigned i=1; i<str_weights.size(); ++i) vals += "," + str_weights[i];
      89           0 :       readInputLine( getShortcutLabel() + "_w: CONSTANT VALUES=" + vals );
      90           0 :     } else error("invalid input for WEIGHTS keyword " + str_weights[0] );
      91             :   } else {
      92           0 :     std::string ones="1"; for(unsigned i=1; i<nat; ++i) ones += ",1";
      93           0 :     readInputLine( getShortcutLabel() + "_w: CONSTANT VALUES=" + ones );
      94             :   }
      95             :   // Read in the instructions on how to compute the center of mass
      96          18 :   bool safe_phases, phases, nopbc; parseFlag("SAFE_PHASES",safe_phases); parseFlag("NOPBC",nopbc);
      97          18 :   if( safe_phases ) phases=true; else parseFlag("PHASES",phases);
      98             :   // This computes a center in the conventional way
      99           9 :   if( !phases || safe_phases ) {
     100             :     // Calculate the sum of the weights
     101           4 :     readInputLine( getShortcutLabel() + "_wnorm: SUM PERIODIC=NO ARG=" + wlab );
     102             :     // Compute the normalised weights
     103           4 :     readInputLine( getShortcutLabel() + "_weights: CUSTOM ARG=" + getShortcutLabel() + "_wnorm," + wlab + " FUNC=y/x PERIODIC=NO");
     104             :     // Get the positions into a multicolvar
     105           3 :     if( phases || nopbc ) readInputLine( getShortcutLabel() + "_pos: POSITION NOPBC ATOMS=" + atlist );
     106           2 :     else readInputLine( getShortcutLabel() + "_pos: POSITION WHOLEMOLECULES ATOMS=" + atlist );
     107             :     // Multiply each vector of positions by the weight
     108           4 :     readInputLine( getShortcutLabel() + "_xwvec: CUSTOM ARG=" + getShortcutLabel() + "_weights," + getShortcutLabel() + "_pos.x FUNC=x*y PERIODIC=NO");
     109           4 :     readInputLine( getShortcutLabel() + "_ywvec: CUSTOM ARG=" + getShortcutLabel() + "_weights," + getShortcutLabel() + "_pos.y FUNC=x*y PERIODIC=NO");
     110           4 :     readInputLine( getShortcutLabel() + "_zwvec: CUSTOM ARG=" + getShortcutLabel() + "_weights," + getShortcutLabel() + "_pos.z FUNC=x*y PERIODIC=NO");
     111             :     // And sum the weighted vectors
     112           4 :     readInputLine( getShortcutLabel() + "_x: SUM ARG=" + getShortcutLabel() + "_xwvec PERIODIC=NO");
     113           4 :     readInputLine( getShortcutLabel() + "_y: SUM ARG=" + getShortcutLabel() + "_ywvec PERIODIC=NO");
     114           4 :     readInputLine( getShortcutLabel() + "_z: SUM ARG=" + getShortcutLabel() + "_zwvec PERIODIC=NO");
     115             :   }
     116             :   // This computes a center using the trigonometric phases
     117           9 :   if( phases ) {
     118             :     // Get the positions into a multicolvar
     119          14 :     readInputLine( getShortcutLabel() + "_fpos: POSITION SCALED_COMPONENTS ATOMS=" + atlist );
     120             :     // Calculate the sines and cosines of the positions and multiply by the weights
     121          14 :     readInputLine( getShortcutLabel() + "_sina: CUSTOM ARG=" + wlab + "," + getShortcutLabel() + "_fpos.a FUNC=x*sin(2*pi*y) PERIODIC=NO");
     122          14 :     readInputLine( getShortcutLabel() + "_cosa: CUSTOM ARG=" + wlab + "," + getShortcutLabel() + "_fpos.a FUNC=x*cos(2*pi*y) PERIODIC=NO");
     123          14 :     readInputLine( getShortcutLabel() + "_sinb: CUSTOM ARG=" + wlab + "," + getShortcutLabel() + "_fpos.b FUNC=x*sin(2*pi*y) PERIODIC=NO");
     124          14 :     readInputLine( getShortcutLabel() + "_cosb: CUSTOM ARG=" + wlab + "," + getShortcutLabel() + "_fpos.b FUNC=x*cos(2*pi*y) PERIODIC=NO");
     125          14 :     readInputLine( getShortcutLabel() + "_sinc: CUSTOM ARG=" + wlab + "," + getShortcutLabel() + "_fpos.c FUNC=x*sin(2*pi*y) PERIODIC=NO");
     126          14 :     readInputLine( getShortcutLabel() + "_cosc: CUSTOM ARG=" + wlab + "," + getShortcutLabel() + "_fpos.c FUNC=x*cos(2*pi*y) PERIODIC=NO");
     127             :     // Sum the sines and cosines
     128          14 :     readInputLine( getShortcutLabel() + "_sinsuma: SUM ARG=" + getShortcutLabel() + "_sina PERIODIC=NO");
     129          14 :     readInputLine( getShortcutLabel() + "_cossuma: SUM ARG=" + getShortcutLabel() + "_cosa PERIODIC=NO");
     130          14 :     readInputLine( getShortcutLabel() + "_sinsumb: SUM ARG=" + getShortcutLabel() + "_sinb PERIODIC=NO");
     131          14 :     readInputLine( getShortcutLabel() + "_cossumb: SUM ARG=" + getShortcutLabel() + "_cosb PERIODIC=NO");
     132          14 :     readInputLine( getShortcutLabel() + "_sinsumc: SUM ARG=" + getShortcutLabel() + "_sinc PERIODIC=NO");
     133          14 :     readInputLine( getShortcutLabel() + "_cossumc: SUM ARG=" + getShortcutLabel() + "_cosc PERIODIC=NO");
     134             :     // And get the final position in fractional coordinates
     135          14 :     readInputLine( getShortcutLabel() + "_a: CUSTOM ARG=" + getShortcutLabel() + "_sinsuma," + getShortcutLabel() + "_cossuma FUNC=atan2(x,y)/(2*pi) PERIODIC=NO");
     136          14 :     readInputLine( getShortcutLabel() + "_b: CUSTOM ARG=" + getShortcutLabel() + "_sinsumb," + getShortcutLabel() + "_cossumb FUNC=atan2(x,y)/(2*pi) PERIODIC=NO");
     137          14 :     readInputLine( getShortcutLabel() + "_c: CUSTOM ARG=" + getShortcutLabel() + "_sinsumc," + getShortcutLabel() + "_cossumc FUNC=atan2(x,y)/(2*pi) PERIODIC=NO");
     138             :     // And create the virtual atom
     139           7 :     if( safe_phases ) {
     140           0 :       readInputLine( getShortcutLabel() + ": ARGS2VATOM XPOS=" + getShortcutLabel() + "_a YPOS=" + getShortcutLabel() + "_b ZPOS=" + getShortcutLabel() + "_c "
     141           0 :                      + " XBKP=" + getShortcutLabel() + "_x YBKP=" + getShortcutLabel() + "_y ZBKP=" + getShortcutLabel() + "_z "
     142           0 :                      + " MASS=" + getShortcutLabel() + "_mass CHARGE=" + getShortcutLabel() + "_charge FRACTIONAL");
     143             :     } else {
     144          21 :       readInputLine( getShortcutLabel() + ": ARGS2VATOM XPOS=" + getShortcutLabel() + "_a YPOS=" + getShortcutLabel() + "_b ZPOS=" + getShortcutLabel() + "_c "
     145          21 :                      + " MASS=" + getShortcutLabel() + "_mass CHARGE=" + getShortcutLabel() + "_charge FRACTIONAL");
     146             :     }
     147             :   } else {
     148             :     // And create the virtual atom
     149           6 :     readInputLine( getShortcutLabel() + ": ARGS2VATOM XPOS=" + getShortcutLabel() + "_x YPOS=" + getShortcutLabel() + "_y ZPOS=" + getShortcutLabel() + "_z "
     150           6 :                    + " MASS=" + getShortcutLabel() + "_mass CHARGE=" + getShortcutLabel() + "_charge ");
     151             :   }
     152        7752 : }
     153             : 
     154             : }
     155             : }