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