Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) crystdistrib 2023-2023 The code team
       3             :    (see the PEOPLE-crystdistrib file at the root of this folder for a list of names)
       4             : 
       5             :    This file is part of crystdistrib code module.
       6             : 
       7             :    The crystdistrib code module is free software: you can redistribute it and/or modify
       8             :    it under the terms of the GNU Lesser General Public License as published by
       9             :    the Free Software Foundation, either version 3 of the License, or
      10             :    (at your option) any later version.
      11             : 
      12             :    The crystdistrib code module is distributed in the hope that it will be useful,
      13             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      14             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      15             :    GNU Lesser General Public License for more details.
      16             : 
      17             :    You should have received a copy of the GNU Lesser General Public License
      18             :    along with the crystdistrib code module.  If not, see <http://www.gnu.org/licenses/>.
      19             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      20             : #include "core/ActionShortcut.h"
      21             : #include "core/ActionRegister.h"
      22             : #include "core/PlumedMain.h"
      23             : #include "core/ActionSet.h"
      24             : #include "core/ActionWithValue.h"
      25             : #include "tools/IFile.h"
      26             : 
      27             : namespace PLMD {
      28             : namespace crystdistrib {
      29             : 
      30             : //+PLUMEDOC COLVAR DOPS
      31             : /*
      32             : Calculate the DOPS order parameter
      33             : 
      34             : \par Examples
      35             : 
      36             : */
      37             : //+ENDPLUMEDOC
      38             : 
      39             : class DopsShortcut : public ActionShortcut {
      40             : public:
      41             :   static void registerKeywords( Keywords& keys );
      42             :   explicit DopsShortcut(const ActionOptions&);
      43             : };
      44             : 
      45             : PLUMED_REGISTER_ACTION(DopsShortcut,"DOPS")
      46             : 
      47           4 : void DopsShortcut::registerKeywords( Keywords& keys ) {
      48           4 :   ActionShortcut::registerKeywords( keys );
      49           8 :   keys.add("atoms","SPECIES","this keyword is used for colvars such as coordination number. In that context it specifies that plumed should calculate "
      50             :            "one coordination number for each of the atoms specified.  Each of these coordination numbers specifies how many of the "
      51             :            "other specified atoms are within a certain cutoff of the central atom.  You can specify the atoms here as another multicolvar "
      52             :            "action or using a MultiColvarFilter or ActionVolume action.  When you do so the quantity is calculated for those atoms specified "
      53             :            "in the previous multicolvar.  This is useful if you would like to calculate the Steinhardt parameter for those atoms that have a "
      54             :            "coordination number more than four for example");
      55           8 :   keys.add("atoms-2","SPECIESA","this keyword is used for colvars such as the coordination number.  In that context it species that plumed should calculate "
      56             :            "one coordination number for each of the atoms specified in SPECIESA.  Each of these cooordination numbers specifies how many "
      57             :            "of the atoms specifies using SPECIESB is within the specified cutoff.  As with the species keyword the input can also be specified "
      58             :            "using the label of another multicolvar");
      59           8 :   keys.add("atoms-2","SPECIESB","this keyword is used for colvars such as the coordination number.  It must appear with SPECIESA.  For a full explanation see "
      60             :            "the documentation for that keyword");
      61           8 :   keys.add("compulsory","KERNELFILE","the file containing the list of kernel parameters.  We expect h, mu and sigma parameters for a 1D Gaussian kernel of the form h*exp(-(x-mu)^2/2sigma^2)");
      62           8 :   keys.add("compulsory","CUTOFF","6.25","to make the calculation faster we calculate a cutoff value on the distances.  The input to this keyword determines x in this expreession max(mu + sqrt(2*x)/sigma)");
      63           4 :   keys.setValueDescription("the values of the DOPS order parameters");
      64          16 :   keys.needsAction("DISTANCE_MATRIX"); keys.needsAction("CUSTOM"); keys.needsAction("ONES"); keys.needsAction("MATRIX_VECTOR_PRODUCT");
      65           4 : }
      66             : 
      67           2 : DopsShortcut::DopsShortcut(const ActionOptions&ao):
      68             :   Action(ao),
      69           2 :   ActionShortcut(ao)
      70             : {
      71             :   // Open a file and read in the kernels
      72           2 :   double cutoff=0, h; std::string kfunc,fname; double dp2cutoff; parse("CUTOFF",dp2cutoff);
      73           4 :   parse("KERNELFILE",fname); IFile ifile; ifile.open(fname);
      74          20 :   for(unsigned k=0;; ++k) {
      75          44 :     if( !ifile.scanField("height",h) ) break;
      76          40 :     std::string ktype; ifile.scanField("kerneltype",ktype); if( ktype!="gaussian" ) error("cannot process kernels of type " + ktype );
      77          60 :     double mu, sigma; ifile.scanField("mu",mu); ifile.scanField("sigma",sigma); ifile.scanField();
      78          20 :     std::string hstr, mustr, sigmastr; Tools::convert( h, hstr );
      79          20 :     Tools::convert( 2*sigma*sigma, sigmastr ); Tools::convert( mu, mustr );
      80             :     // Get a sensible value for the cutoff
      81          20 :     double support = sqrt(2.0*dp2cutoff)*(1.0/sigma);
      82          20 :     if( mu+support>cutoff ) cutoff= mu + support;
      83             :     // And make the kernel
      84          38 :     if( k==0 ) kfunc = hstr; else kfunc += "+" + hstr;
      85          40 :     kfunc += "*exp(-(x-" + mustr +")^2/" + sigmastr + ")";
      86          20 :   }
      87             :   std::string sp_str, specA, specB, grpinfo;
      88           8 :   parse("SPECIES",sp_str); parse("SPECIESA",specA); parse("SPECIESB",specB);
      89           2 :   if( sp_str.length()>0 ) {
      90           4 :     grpinfo="GROUP=" + sp_str;
      91             :   } else {
      92           0 :     if( specA.length()==0 || specB.length()==0 ) error("no atoms were specified in input use either SPECIES or SPECIESA + SPECIESB");
      93           0 :     grpinfo="GROUPA=" + specA + " GROUPB=" + specB;
      94             :   }
      95           2 :   std::string cutstr; Tools::convert( cutoff, cutstr );
      96             :   // Setup the contact matrix
      97           4 :   readInputLine( getShortcutLabel() + "_cmat: DISTANCE_MATRIX  " + grpinfo + " CUTOFF=" + cutstr);
      98             :   // And the kernels
      99           4 :   readInputLine( getShortcutLabel() + "_kval: CUSTOM ARG=" + getShortcutLabel() + "_cmat PERIODIC=NO FUNC=" + kfunc );
     100             :   // Find the number of ones we need to multiply by
     101           2 :   ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmat");
     102           2 :   plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
     103           2 :   std::string size; Tools::convert( (av->copyOutput(0))->getShape()[1], size );
     104           4 :   readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
     105             :   // And the final order parameters
     106           4 :   readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_kval," + getShortcutLabel() + "_ones");
     107           4 : }
     108             : 
     109             : }
     110             : }
     111             : 
     112             : 
     113             :