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           4 :   keys.needsAction("DISTANCE_MATRIX");
      65           4 :   keys.needsAction("CUSTOM");
      66           4 :   keys.needsAction("ONES");
      67           4 :   keys.needsAction("MATRIX_VECTOR_PRODUCT");
      68           4 : }
      69             : 
      70           2 : DopsShortcut::DopsShortcut(const ActionOptions&ao):
      71             :   Action(ao),
      72           2 :   ActionShortcut(ao) {
      73             :   // Open a file and read in the kernels
      74           2 :   double cutoff=0, h;
      75             :   std::string kfunc,fname;
      76             :   double dp2cutoff;
      77           2 :   parse("CUTOFF",dp2cutoff);
      78           2 :   parse("KERNELFILE",fname);
      79           2 :   IFile ifile;
      80           2 :   ifile.open(fname);
      81          20 :   for(unsigned k=0;; ++k) {
      82          44 :     if( !ifile.scanField("height",h) ) {
      83             :       break;
      84             :     }
      85             :     std::string ktype;
      86          40 :     ifile.scanField("kerneltype",ktype);
      87          20 :     if( ktype!="gaussian" ) {
      88           0 :       error("cannot process kernels of type " + ktype );
      89             :     }
      90             :     double mu, sigma;
      91          20 :     ifile.scanField("mu",mu);
      92          20 :     ifile.scanField("sigma",sigma);
      93          20 :     ifile.scanField();
      94             :     std::string hstr, mustr, sigmastr;
      95          20 :     Tools::convert( h, hstr );
      96          20 :     Tools::convert( 2*sigma*sigma, sigmastr );
      97          20 :     Tools::convert( mu, mustr );
      98             :     // Get a sensible value for the cutoff
      99          20 :     double support = sqrt(2.0*dp2cutoff)*(1.0/sigma);
     100          20 :     if( mu+support>cutoff ) {
     101           6 :       cutoff= mu + support;
     102             :     }
     103             :     // And make the kernel
     104          20 :     if( k==0 ) {
     105             :       kfunc = hstr;
     106             :     } else {
     107          36 :       kfunc += "+" + hstr;
     108             :     }
     109          40 :     kfunc += "*exp(-(x-" + mustr +")^2/" + sigmastr + ")";
     110          20 :   }
     111             :   std::string sp_str, specA, specB, grpinfo;
     112           2 :   parse("SPECIES",sp_str);
     113           2 :   parse("SPECIESA",specA);
     114           4 :   parse("SPECIESB",specB);
     115           2 :   if( sp_str.length()>0 ) {
     116           4 :     grpinfo="GROUP=" + sp_str;
     117             :   } else {
     118           0 :     if( specA.length()==0 || specB.length()==0 ) {
     119           0 :       error("no atoms were specified in input use either SPECIES or SPECIESA + SPECIESB");
     120             :     }
     121           0 :     grpinfo="GROUPA=" + specA + " GROUPB=" + specB;
     122             :   }
     123             :   std::string cutstr;
     124           2 :   Tools::convert( cutoff, cutstr );
     125             :   // Setup the contact matrix
     126           4 :   readInputLine( getShortcutLabel() + "_cmat: DISTANCE_MATRIX  " + grpinfo + " CUTOFF=" + cutstr);
     127             :   // And the kernels
     128           4 :   readInputLine( getShortcutLabel() + "_kval: CUSTOM ARG=" + getShortcutLabel() + "_cmat PERIODIC=NO FUNC=" + kfunc );
     129             :   // Find the number of ones we need to multiply by
     130           2 :   ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmat");
     131           2 :   plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
     132             :   std::string size;
     133           2 :   Tools::convert( (av->copyOutput(0))->getShape()[1], size );
     134           4 :   readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
     135             :   // And the final order parameters
     136           4 :   readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_kval," + getShortcutLabel() + "_ones");
     137           4 : }
     138             : 
     139             : }
     140             : }
     141             : 
     142             : 
     143             :