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/ActionWithValue.h"
      21             : #include "core/ActionShortcut.h"
      22             : 
      23             : namespace PLMD {
      24             : namespace colvar {
      25             : 
      26             : //+PLUMEDOC COLVAR GYRATION
      27             : /*
      28             : Calculate the radius of gyration, or other properties related to it.
      29             : 
      30             : With this version of the command you can use any property you so choose to define the weights that are used when computing the average.
      31             : If you use the mass or if all the atoms are ascribed weights of one PLUMED defaults to \ref GYRATION_FAST
      32             : 
      33             : \par Examples
      34             : 
      35             : */
      36             : //+ENDPLUMEDOC
      37             : 
      38             : //+PLUMEDOC MCOLVAR GYRATION_TENSOR
      39             : /*
      40             : Calculate the gyration tensor using a user specified vector of weights
      41             : 
      42             : \par Examples
      43             : 
      44             : */
      45             : //+ENDPLUMEDOC
      46             : 
      47             : class GyrationShortcut : public ActionShortcut {
      48             : public:
      49             :   static void registerKeywords( Keywords& keys );
      50             :   explicit GyrationShortcut(const ActionOptions&);
      51             : };
      52             : 
      53             : PLUMED_REGISTER_ACTION(GyrationShortcut,"GYRATION")
      54             : PLUMED_REGISTER_ACTION(GyrationShortcut,"GYRATION_TENSOR")
      55             : 
      56         122 : void GyrationShortcut::registerKeywords( Keywords& keys ) {
      57         122 :   ActionShortcut::registerKeywords( keys );
      58         244 :   keys.add("atoms","ATOMS","the group of atoms that you are calculating the Gyration Tensor for");
      59         244 :   keys.add("compulsory","TYPE","RADIUS","The type of calculation relative to the Gyration Tensor you want to perform");
      60         244 :   keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
      61         244 :   keys.add("optional","WEIGHTS","what weights should be used when calculating the center.  If this keyword is not present the geometric center is computed. "
      62             :            "If WEIGHTS=@Masses is used the center of mass is computed.  If WEIGHTS=@charges the center of charge is computed.  If "
      63             :            "the label of an action is provided PLUMED assumes that that action calculates a list of symmetry functions that can be used "
      64             :            "as weights. Lastly, an explicit list of numbers to use as weights can be provided");
      65         244 :   keys.addFlag("PHASES",false,"use trigonometric phases when computing position of center of mass");
      66         244 :   keys.addFlag("MASS",false,"calculate the center of mass");
      67         244 :   keys.addFlag("MASS_WEIGHTED",false,"set the masses of all the atoms equal to one");
      68         244 :   keys.addFlag("UNORMALIZED",false,"do not divide by the sum of the weights");
      69         122 :   keys.setValueDescription("the radius that was computed from the weights");
      70         122 :   keys.addActionNameSuffix("_FAST");
      71         122 :   keys.needsAction("CENTER");
      72         122 :   keys.needsAction("CONSTANT");
      73         122 :   keys.needsAction("ONES");
      74         122 :   keys.needsAction("MASS");
      75         122 :   keys.needsAction("DISTANCE");
      76         122 :   keys.needsAction("COVARIANCE_MATRIX");
      77         122 :   keys.needsAction("SELECT_COMPONENTS");
      78         122 :   keys.needsAction("SUM");
      79         122 :   keys.needsAction("CUSTOM");
      80         122 :   keys.needsAction("DIAGONALIZE");
      81         122 : }
      82             : 
      83         116 : GyrationShortcut::GyrationShortcut(const ActionOptions& ao):
      84             :   Action(ao),
      85         116 :   ActionShortcut(ao) {
      86             :   bool usemass, phases;
      87         116 :   parseFlag("MASS",usemass);
      88         232 :   parseFlag("PHASES",phases);
      89             :   std::vector<std::string> str_weights;
      90         232 :   parseVector("WEIGHTS",str_weights);
      91             :   std::string wflab;
      92         116 :   if( !phases ) {
      93         115 :     if( usemass || str_weights.size()==0 || (str_weights.size()==1 && str_weights[0]=="@Masses") ) {
      94             :       std::string wt_str;
      95         112 :       if( str_weights.size()>0 ) {
      96           0 :         wt_str="WEIGHTS=" + str_weights[0];
      97           0 :         for(unsigned i=1; i<str_weights.size(); ++i) {
      98           0 :           wt_str += "," + str_weights[i];
      99             :         }
     100             :       }
     101         112 :       if( usemass || (str_weights.size()==1 && str_weights[0]=="@Masses") ) {
     102             :         wt_str = "MASS";
     103             :       }
     104         232 :       readInputLine( getShortcutLabel() + ": GYRATION_FAST " + wt_str + " " + convertInputLineToString() );
     105             :       return;
     106             :     }
     107             :   }
     108           4 :   if( usemass ) {
     109           0 :     str_weights.resize(1);
     110             :     str_weights[0]="@Masses";
     111             :   }
     112          10 :   log<<"  Bibliography "<<plumed.cite("Jirí Vymetal and Jirí Vondrasek, J. Phys. Chem. A 115, 11455 (2011)")<<"\n";
     113             :   // Read in the atoms involved
     114             :   std::vector<std::string> atoms;
     115           4 :   parseVector("ATOMS",atoms);
     116           4 :   Tools::interpretRanges(atoms);
     117           4 :   std::string gtype, atlist=atoms[0];
     118          12 :   for(unsigned i=1; i<atoms.size(); ++i) {
     119          16 :     atlist += "," + atoms[i];
     120             :   }
     121             :   bool nopbc;
     122           8 :   parseFlag("NOPBC",nopbc);
     123             :   std::string pbcstr;
     124           4 :   if(nopbc) {
     125             :     pbcstr = " NOPBC";
     126             :   }
     127             :   std::string phasestr;
     128           4 :   if(phases) {
     129             :     phasestr = " PHASES";
     130             :   }
     131             :   // Create the geometric center of the molecule
     132           4 :   std::string weights_str=" WEIGHTS=" + str_weights[0];
     133           8 :   for(unsigned i=1; i<str_weights.size(); ++i) {
     134           8 :     weights_str += "," + str_weights[i];
     135             :   }
     136           8 :   readInputLine( getShortcutLabel() + "_cent: CENTER ATOMS=" + atlist + pbcstr + phasestr + weights_str );
     137           4 :   if( str_weights.size()==0 ) {
     138           0 :     wflab = getShortcutLabel() + "_w";
     139             :     std::string str_natoms;
     140           0 :     Tools::convert( atoms.size(), str_natoms );
     141           0 :     readInputLine( getShortcutLabel() + "_w: ONES SIZE=" + str_natoms );
     142           6 :   } else if( str_weights.size()==1 && str_weights[0]=="@Masses" ) {
     143           0 :     wflab = getShortcutLabel() + "_m";
     144           0 :     readInputLine( getShortcutLabel() + "_m: MASS ATOMS=" + atlist );
     145           4 :   } else if( str_weights.size()>1 ) {
     146           2 :     std::string vals=str_weights[0];
     147           6 :     for(unsigned i=1; i<str_weights.size(); ++i) {
     148           8 :       vals += "," + str_weights[i];
     149             :     }
     150           4 :     readInputLine( getShortcutLabel() + "_w: CONSTANT VALUES=" + vals );
     151           4 :     wflab=getShortcutLabel() + "_w";
     152             :   } else {
     153           2 :     plumed_assert( str_weights.size()==1 );
     154           2 :     wflab = getShortcutLabel() + "_cent_w";
     155           2 :     ActionWithValue* av=plumed.getActionSet().selectWithLabel<ActionWithValue*>( wflab );
     156           2 :     if( !av ) {
     157             :       wflab = str_weights[0];
     158             :     }
     159             :   }
     160             :   // Check for normalisation
     161             :   bool unorm;
     162           8 :   parseFlag("UNORMALIZED",unorm);
     163             :   // Find out the type
     164           4 :   if( getName()!="GYRATION_TENSOR" ) {
     165           0 :     parse("TYPE",gtype);
     166           0 :     if( gtype!="RADIUS" && gtype!="TRACE" && gtype!="GTPC_1" && gtype!="GTPC_2" && gtype!="GTPC_3" && gtype!="ASPHERICITY" && gtype!="ACYLINDRICITY"
     167           0 :         && gtype!= "KAPPA2" && gtype!="RGYR_1" && gtype!="RGYR_2" && gtype!="RGYR_3" ) {
     168           0 :       error("type " + gtype + " is invalid");
     169             :     }
     170             :     // Check if we need to calculate the unormlised radius
     171           0 :     if( gtype=="TRACE" || gtype=="KAPPA2" ) {
     172           0 :       unorm=true;
     173             :     }
     174             :   }
     175             :   // Compute all the vectors separating all the positions from the center
     176           4 :   std::string distance_act = getShortcutLabel() + "_dists: DISTANCE COMPONENTS" + pbcstr;
     177          16 :   for(unsigned i=0; i<atoms.size(); ++i) {
     178             :     std::string num;
     179          12 :     Tools::convert( i+1, num );
     180          24 :     distance_act += " ATOMS" + num + "=" + getShortcutLabel() + "_cent," + atoms[i];
     181             :   }
     182           4 :   readInputLine( distance_act );
     183             :   // And calculate the covariance
     184             :   std::string norm_str;
     185           4 :   if( unorm ) {
     186             :     norm_str = " UNORMALIZED";
     187             :   }
     188           4 :   if( getName()=="GYRATION_TENSOR" ) {
     189           8 :     readInputLine( getShortcutLabel() + ": COVARIANCE_MATRIX ARG=" + getShortcutLabel() + "_dists.x," + getShortcutLabel() + "_dists.y," + getShortcutLabel() + "_dists.z WEIGHTS=" + wflab + norm_str );
     190             :     return;
     191             :   }
     192           0 :   readInputLine( getShortcutLabel() + "_tensor: COVARIANCE_MATRIX ARG=" + getShortcutLabel() + "_dists.x," + getShortcutLabel() + "_dists.y," + getShortcutLabel() + "_dists.z WEIGHTS=" + wflab + norm_str );
     193             :   // Pick out the diagonal elements
     194           0 :   readInputLine( getShortcutLabel() + "_diag_elements: SELECT_COMPONENTS ARG=" + getShortcutLabel() + "_tensor COMPONENTS=1.1,2.2,3.3");
     195           0 :   if( gtype=="RADIUS") {
     196             :     // And now we need the average trace for the gyration radius
     197           0 :     readInputLine( getShortcutLabel() + "_trace: SUM ARG=" + getShortcutLabel() + "_diag_elements PERIODIC=NO");
     198             :     // Square root the radius
     199           0 :     readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_trace FUNC=sqrt(x) PERIODIC=NO");
     200           0 :   } else if( gtype=="TRACE" ) {
     201             :     // Compte the trace of the gyration tensor
     202           0 :     readInputLine( getShortcutLabel() + "_trace: SUM ARG=" + getShortcutLabel() + "_diag_elements PERIODIC=NO");
     203             :     // And double it
     204           0 :     readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_trace FUNC=2*x PERIODIC=NO");
     205             :   } else {
     206             :     // Diagonalize the gyration tensor
     207           0 :     readInputLine( getShortcutLabel() + "_diag: DIAGONALIZE ARG=" + getShortcutLabel() + "_tensor VECTORS=all" );
     208           0 :     if( gtype.find("GTPC")!=std::string::npos ) {
     209           0 :       std::size_t und=gtype.find_first_of("_");
     210           0 :       if( und==std::string::npos ) {
     211           0 :         error( gtype + " is not a valid type for gyration radius");
     212             :       }
     213           0 :       std::string num = gtype.substr(und+1);
     214           0 :       if( num!="1" && num!="2" && num!="3" ) {
     215           0 :         error( gtype + " is not a valid type for gyration radius");
     216             :       }
     217             :       // Now get the appropriate eigenvalue
     218           0 :       readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-" + num + " FUNC=sqrt(x) PERIODIC=NO");
     219           0 :     } else if( gtype.find("RGYR")!=std::string::npos ) {
     220           0 :       std::size_t und=gtype.find_first_of("_");
     221           0 :       if( und==std::string::npos ) {
     222           0 :         error( gtype + " is not a valid type for gyration radius");
     223             :       }
     224             :       unsigned ind;
     225           0 :       Tools::convert( gtype.substr(und+1), ind );
     226             :       // Now get the appropriate quantity
     227           0 :       if( ind==3 ) {
     228           0 :         readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-2 FUNC=sqrt(x+y) PERIODIC=NO");
     229           0 :       } else if( ind==2 ) {
     230           0 :         readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x+y) PERIODIC=NO");
     231           0 :       } else if( ind==1 ) {
     232           0 :         readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x+y) PERIODIC=NO");
     233             :       } else {
     234           0 :         error( gtype + " is not a valid type for gyration radius");
     235             :       }
     236           0 :     } else if( gtype=="ASPHERICITY" ) {
     237           0 :       readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x-0.5*(y+z)) PERIODIC=NO" );
     238           0 :     } else if( gtype=="ACYLINDRICITY" ) {
     239           0 :       readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x-y) PERIODIC=NO" );
     240           0 :     } else if( gtype=="KAPPA2" ) {
     241           0 :       readInputLine( getShortcutLabel() + "_numer: CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=x*y+x*z+y*z PERIODIC=NO" );
     242           0 :       readInputLine( getShortcutLabel() + "_denom: CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=x+y+z PERIODIC=NO" );
     243           0 :       readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_numer," + getShortcutLabel() + "_denom FUNC=1-3*(x/(y*y)) PERIODIC=NO");
     244             :     } else {
     245           0 :       error( gtype + " is not a valid type for gyration radius");
     246             :     }
     247             :   }
     248         124 : }
     249             : 
     250             : }
     251             : }