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 DOPS order parameter for molecules.
      33             : 
      34             : DOPS is a shortcut to calculate the Distance Order Parameters that are described in the paper below.
      35             : As arguments, DOPS takes a list of atoms, a cutoff, and a kernel file detailing the means, variances, and normalization factors of probability distributions. DOPS returns a vector of order parameters.
      36             : 
      37             : The DOPS kernel file has FIELDS height, mu, and sigma corresponding to the normalization factor, mean, and variance of the gaussian distributions used in the order parameters. The SET kerneltype is gaussian.
      38             : 
      39             : DOPS returns one order parameter per atom given, evaluated over each atom's neighbors within the cutoff given. The kernel file defines a radial distribution function. The order parameter is obtained by evaluating the RDF at each neighbor's position within the cutoff, and summing them up.
      40             : 
      41             : This example file calculates the DOPS for a system of 3 molecules.
      42             : ```plumed
      43             : #SETTINGS INPUTFILES=regtest/crystdistrib/rt-dops-forces/kernels.dat
      44             : dops: DOPS SPECIES=1,4,7 CUTOFF=7.4 KERNELFILE=regtest/crystdistrib/rt-dops-forces/kernels.dat
      45             : ```
      46             : 
      47             : */
      48             : //+ENDPLUMEDOC
      49             : 
      50             : 
      51             : class DopsShortcut : public ActionShortcut {
      52             : public:
      53             :   static void registerKeywords( Keywords& keys );
      54             :   explicit DopsShortcut(const ActionOptions&);
      55             : };
      56             : 
      57             : PLUMED_REGISTER_ACTION(DopsShortcut,"DOPS")
      58             : 
      59           4 : void DopsShortcut::registerKeywords( Keywords& keys ) {
      60           4 :   ActionShortcut::registerKeywords( keys );
      61           4 :   keys.add("atoms","SPECIES","this keyword is used for colvars such as coordination number. In that context it specifies that plumed should calculate "
      62             :            "one coordination number for each of the atoms specified.  Each of these coordination numbers specifies how many of the "
      63             :            "other specified atoms are within a certain cutoff of the central atom.  You can specify the atoms here as another multicolvar "
      64             :            "action or using a MultiColvarFilter or ActionVolume action.  When you do so the quantity is calculated for those atoms specified "
      65             :            "in the previous multicolvar.  This is useful if you would like to calculate the Steinhardt parameter for those atoms that have a "
      66             :            "coordination number more than four for example");
      67           4 :   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 "
      68             :            "one coordination number for each of the atoms specified in SPECIESA.  Each of these cooordination numbers specifies how many "
      69             :            "of the atoms specifies using SPECIESB is within the specified cutoff.  As with the species keyword the input can also be specified "
      70             :            "using the label of another multicolvar");
      71           4 :   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 "
      72             :            "the documentation for that keyword");
      73           4 :   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)");
      74           4 :   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)");
      75           8 :   keys.setValueDescription("vector","the values of the DOPS order parameters");
      76           4 :   keys.needsAction("DISTANCE_MATRIX");
      77           4 :   keys.needsAction("CUSTOM");
      78           4 :   keys.needsAction("ONES");
      79           4 :   keys.needsAction("MATRIX_VECTOR_PRODUCT");
      80           4 :   keys.addDOI("10.1063/1.3548889");
      81           4 : }
      82             : 
      83           2 : DopsShortcut::DopsShortcut(const ActionOptions&ao):
      84             :   Action(ao),
      85           2 :   ActionShortcut(ao) {
      86             :   // Open a file and read in the kernels
      87           2 :   double cutoff=0, h;
      88             :   std::string kfunc,fname;
      89             :   double dp2cutoff;
      90           2 :   parse("CUTOFF",dp2cutoff);
      91           2 :   parse("KERNELFILE",fname);
      92           2 :   IFile ifile;
      93           2 :   ifile.open(fname);
      94          20 :   for(unsigned k=0;; ++k) {
      95          44 :     if( !ifile.scanField("height",h) ) {
      96             :       break;
      97             :     }
      98             :     std::string ktype;
      99          40 :     ifile.scanField("kerneltype",ktype);
     100          20 :     if( ktype!="gaussian" ) {
     101           0 :       error("cannot process kernels of type " + ktype );
     102             :     }
     103             :     double mu, sigma;
     104          20 :     ifile.scanField("mu",mu);
     105          20 :     ifile.scanField("sigma",sigma);
     106          20 :     ifile.scanField();
     107             :     std::string hstr, mustr, sigmastr;
     108          20 :     Tools::convert( h, hstr );
     109          20 :     Tools::convert( 2*sigma*sigma, sigmastr );
     110          20 :     Tools::convert( mu, mustr );
     111             :     // Get a sensible value for the cutoff
     112          20 :     double support = sqrt(2.0*dp2cutoff)*(1.0/sigma);
     113          20 :     if( mu+support>cutoff ) {
     114           6 :       cutoff= mu + support;
     115             :     }
     116             :     // And make the kernel
     117          20 :     if( k==0 ) {
     118             :       kfunc = hstr;
     119             :     } else {
     120          36 :       kfunc += "+" + hstr;
     121             :     }
     122          40 :     kfunc += "*exp(-(x-" + mustr +")^2/" + sigmastr + ")";
     123          20 :   }
     124             :   std::string sp_str, specA, specB, grpinfo;
     125           2 :   parse("SPECIES",sp_str);
     126           2 :   parse("SPECIESA",specA);
     127           4 :   parse("SPECIESB",specB);
     128           2 :   if( sp_str.length()>0 ) {
     129           4 :     grpinfo="GROUP=" + sp_str;
     130             :   } else {
     131           0 :     if( specA.length()==0 || specB.length()==0 ) {
     132           0 :       error("no atoms were specified in input use either SPECIES or SPECIESA + SPECIESB");
     133             :     }
     134           0 :     grpinfo="GROUPA=" + specA + " GROUPB=" + specB;
     135             :   }
     136             :   std::string cutstr;
     137           2 :   Tools::convert( cutoff, cutstr );
     138             :   // Setup the contact matrix
     139           4 :   readInputLine( getShortcutLabel() + "_cmat: DISTANCE_MATRIX  " + grpinfo + " CUTOFF=" + cutstr);
     140             :   // And the kernels
     141           4 :   readInputLine( getShortcutLabel() + "_kval: CUSTOM ARG=" + getShortcutLabel() + "_cmat PERIODIC=NO FUNC=" + kfunc );
     142             :   // Find the number of ones we need to multiply by
     143           2 :   ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmat");
     144           2 :   plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
     145             :   std::string size;
     146           2 :   Tools::convert( (av->copyOutput(0))->getShape()[1], size );
     147           4 :   readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
     148             :   // And the final order parameters
     149           4 :   readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_kval," + getShortcutLabel() + "_ones");
     150           4 : }
     151             : 
     152             : }
     153             : }
     154             : 
     155             : 
     156             :