Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2012-2017 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "RDF.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "core/ActionShortcut.h"
      25             : 
      26             : //+PLUMEDOC MCOLVAR PAIRENTROPY
      27             : /*
      28             : Calculate the KL Entropy from the radial distribution function
      29             : 
      30             : This shortcut provides an implementation of the CV that is described in the paper that is cited below. Although, the implementation
      31             : that is offered here is not the one that was used to produce the results described in that paper the values produced by this
      32             : implementation have been tested against the implementation that was used in the paper, which can be found [here](https://sites.google.com/site/pablompiaggi/scripts/pair-entropy/pair-entropy-cv?authuser=0).
      33             : __If you want to use this implementation for a production application we would recommend further testing.__
      34             : 
      35             : An example input that calculates and prints the PAIRENTROPY CV is shown below:
      36             : 
      37             : ```plumed
      38             : pp: PAIRENTROPY ...
      39             :    GROUP=1-108 MAXR=2.0 GRID_BIN=20
      40             :    BANDWIDTH=0.13 KERNEL=gaussian
      41             : ...
      42             : PRINT ARG=pp FILE=colvar
      43             : ```
      44             : 
      45             : By expanding the shortct in the input above you can see how features that are already avaialble within PLUMED can be reused to calculate this CV.  The resulting implementation is likely slower than the
      46             : direct implementation that is available [here](https://sites.google.com/site/pablompiaggi/scripts/pair-entropy/pair-entropy-cv?authuser=0).  We hope, however, that this implementation helps others to understand
      47             : how this CV is constructed.
      48             : 
      49             : If you would like to calculate the pair entropy based on the radial distribution of the atoms in GROUPB around GROUPA you can use an input similar to the one shown below:
      50             : 
      51             : ```plumed
      52             : pp: PAIRENTROPY ...
      53             :    GROUPA=1-108 GROUPB=109-300
      54             :    GRID_BIN=20 MAXR=2.0
      55             :    BANDWIDTH=0.13 KERNEL=gaussian
      56             :    DENSITY=1.0
      57             : ...
      58             : PRINT ARG=pp FILE=colvar
      59             : ```
      60             : 
      61             : Notice that we have also used the `DENSITY` keyword to set the background density that is used when normalizing the radial distribution function explicity to 1 atom$/nm^{3}$.
      62             : When this keyword is not used, this density is calculated by dividing the number of atoms by the volume of the box as you can see if you expand the shortcut in the
      63             : first input above.
      64             : 
      65             : */
      66             : //+ENDPLUMEDOC
      67             : 
      68             : //+PLUMEDOC MCOLVAR PAIRENTROPIES
      69             : /*
      70             : Calculate the KL entropy from the RDF around each of the atoms
      71             : 
      72             : This shortcut can be used to calculate the descriptors of local structure that are described in the paper that is cited below.  This is not the implementation
      73             : that was used when that paper was written.  However, the values that this implementation generates were tested against the values that the implementation that
      74             : was used in the paper, which can be found [here](https://sites.google.com/site/pablompiaggi/scripts/pair-entropy/pair-entropy-fingerprint?authuser=0).
      75             : __If you want to use this implementation for a production application we would recommend further testing.__
      76             : 
      77             : The following example illustrates how this shortcut can be used to calculate and print the descriptor for a set of 64 atoms:
      78             : 
      79             : ```plumed
      80             : pp: PAIRENTROPIES ...
      81             :    GROUP=1-64 MAXR=2
      82             :    GRID_BIN=20 BANDWIDTH=0.13
      83             :    KERNEL=gaussian
      84             :    DENSITY=1.0
      85             : ...
      86             : DUMPATOMS ATOMS=1-64 ARG=pp FILE=p_entropies.xyz
      87             : ```
      88             : 
      89             : If you expand the shortcut you will notice that this shortcut creates __a lot__ of actions.  It is thus likely to be considerably slower than the implementation
      90             : [here](https://sites.google.com/site/pablompiaggi/scripts/pair-entropy/pair-entropy-fingerprint?authuser=0). However, we hope that the implementation here is a
      91             : useful reference implementation that anyone interested in using this method can use to test their implementations of it. In addition, we hope that the shortcut
      92             : illustrates how the calculation of this complex descriptor can be broken down into a set of simpler calculations.
      93             : 
      94             : Notice also that you do not need to use the DENSITY keyword to set the density as has been done above.  You can instead use an input like the one below:
      95             : 
      96             : ```plumed
      97             : pp: PAIRENTROPIES ...
      98             :    GROUP=1-64 MAXR=2
      99             :    GRID_BIN=20 BANDWIDTH=0.13
     100             :    KERNEL=gaussian
     101             : ...
     102             : DUMPATOMS ATOMS=1-64 ARG=pp FILE=p_entropies.xyz
     103             : ```
     104             : 
     105             : As you can see, if you expand the shortcut in the above input, the density in this case is calculated by dividing the number of atoms by the volume of the box.
     106             : 
     107             : */
     108             : //+ENDPLUMEDOC
     109             : 
     110             : namespace PLMD {
     111             : namespace gridtools {
     112             : 
     113             : class PairEntropy : public ActionShortcut {
     114             : public:
     115             :   static void registerKeywords( Keywords& keys );
     116             :   explicit PairEntropy(const ActionOptions&ao);
     117             : };
     118             : 
     119             : PLUMED_REGISTER_ACTION(PairEntropy,"PAIRENTROPY")
     120             : PLUMED_REGISTER_ACTION(PairEntropy,"PAIRENTROPIES")
     121             : 
     122          13 : void PairEntropy::registerKeywords( Keywords& keys ) {
     123          13 :   RDF::registerKeywords( keys );
     124          13 :   keys.needsAction("RDF");
     125          13 :   keys.addDeprecatedKeyword("SIGMA","BANDWIDTH");
     126          13 :   if( keys.getDisplayName()=="PAIRENTROPIES" ) {
     127          10 :     keys.setValueDescription("grid","the KL-entropy that is computed from the radial distribution function of each atom (a function on a grid with one point per atom)");
     128             :   } else {
     129          16 :     keys.setValueDescription("scalar","the KL-entropy that is computed from the radial distribution function");
     130             :   }
     131          13 :   keys.addDOI("10.1103/PhysRevLett.119.015701");
     132          13 :   keys.needsAction("INTERPOLATE_GRID");
     133          13 :   keys.needsAction("INTEGRATE_GRID");
     134          13 :   keys.needsAction("MARGINAL");
     135          13 : }
     136             : 
     137           5 : PairEntropy::PairEntropy(const ActionOptions&ao):
     138             :   Action(ao),
     139           5 :   ActionShortcut(ao) {
     140             :   std::string ref_str, ref_name;
     141          10 :   parse("REFERENCE",ref_name);
     142           5 :   if( ref_name.length()>0 ) {
     143           0 :     ref_str = "REFERENCE=" + ref_name;
     144             :   } else {
     145          11 :     ref_name = getShortcutLabel() + "_rdf";
     146             :   }
     147             :   // Read in the atoms and get the number of atoms that we are using
     148             :   std::string atom_str, group_str;
     149          10 :   parse("GROUP",group_str);
     150           5 :   if( group_str.length()>0 ) {
     151          11 :     atom_str="GROUP=" + group_str;
     152             :   } else {
     153             :     std::string groupa_str, groupb_str;
     154           0 :     parse("GROUPA",groupa_str);
     155           0 :     parse("GROUPB",groupb_str);
     156           0 :     atom_str="GROUPA=" + groupa_str + " GROUPB=" + groupb_str;
     157             :   }
     158             :   // Read in all other keywords and create the RDF object
     159           6 :   std::string maxr, nbins, dens, bw="", kernel;
     160           5 :   parse("MAXR",maxr);
     161           5 :   parse("GRID_BIN",nbins);
     162           5 :   parse("DENSITY",dens);
     163           5 :   parse("BANDWIDTH",bw);
     164          11 :   parse("KERNEL",kernel);
     165             :   std::string dens_str;
     166           5 :   if( dens.length()>0 ) {
     167           2 :     dens_str = " DENSITY=" + dens;
     168             :   }
     169           6 :   std::string noaverage="";
     170           5 :   if( getName()=="PAIRENTROPIES" ) {
     171             :     noaverage = " NO_AVERAGE";
     172             :   }
     173          30 :   readInputLine( getShortcutLabel() + "_rdf: RDF " + atom_str + " KERNEL={" + kernel + "} GRID_BIN=" + nbins + " MAXR=" + maxr + dens_str + " BANDWIDTH=" + bw + " " + ref_str + noaverage);
     174             :   // Get the shape of the matrix
     175           5 :   std::vector<std::string> shape_str(2);
     176           4 :   RDF::getDistanceMatrixShape( getShortcutLabel() + "_rdf_mat", this, shape_str );
     177             :   // And compute the two functions we are integrating (we use two matheval objects here and sum them in order to avoid nans from taking logarithms of zero)
     178          12 :   readInputLine( getShortcutLabel() + "_conv_t1: CUSTOM ARG=" + getShortcutLabel() + "_rdf," + ref_name + "_x2 FUNC=x*y*log(x) PERIODIC=NO");
     179          12 :   readInputLine( getShortcutLabel() + "_conv_t2: CUSTOM ARG=" + getShortcutLabel() + "_rdf," + ref_name + "_x2 FUNC=(1-x)*y PERIODIC=NO");
     180          12 :   readInputLine( getShortcutLabel() + "_conv: CUSTOM ARG=" + getShortcutLabel() + "_conv_t1," + getShortcutLabel() + "_conv_t2 FUNC=x+y PERIODIC=NO");
     181             :   // Calculate number of midpoints
     182             :   unsigned nb;
     183           4 :   Tools::convert( nbins, nb );
     184           4 :   Tools::convert( nb-1, nbins );
     185             :   // Calculate grid spacing
     186             :   double rmax;
     187           4 :   Tools::convert( maxr, rmax );
     188           4 :   double spacing = rmax / nb;
     189             :   std::string fmid, lmid;
     190           4 :   Tools::convert( 0.5*spacing, fmid );
     191           4 :   Tools::convert( rmax-0.5*spacing, lmid );
     192           4 :   if( getName()=="PAIRENTROPIES" ) {
     193             :     int iatoms;
     194           1 :     Tools::convert( shape_str[0], iatoms );
     195             :     std::string num;
     196           1 :     Tools::convert( iatoms-1, num );
     197             :     // Now integrate using trapezium rule
     198           6 :     readInputLine( getShortcutLabel() + "_midp: INTERPOLATE_GRID ARG=" + getShortcutLabel() + "_conv INTERPOLATION_TYPE=linear GRID_MIN=" + fmid + ",0 GRID_MAX=" + lmid + "," + num + " GRID_BIN=" + nbins + "," + num ); // First interpolate onto midpoints
     199           3 :     readInputLine( getShortcutLabel() + "_int: MARGINAL ARG=" + getShortcutLabel() + "_midp DIR=2");
     200             :   } else {
     201             :     // Now integrate using trapezium rule
     202          12 :     readInputLine( getShortcutLabel() + "_midp: INTERPOLATE_GRID ARG=" + getShortcutLabel() + "_conv INTERPOLATION_TYPE=linear GRID_MIN=" + fmid +" GRID_MAX=" + lmid + " GRID_BIN=" + nbins ); // First interpolate onto midpoints
     203           9 :     readInputLine( getShortcutLabel() + "_int: INTEGRATE_GRID ARG=" + getShortcutLabel() + "_midp PERIODIC=NO"); // And then integrate
     204             :   }
     205             :   // And multiply by final normalizing constant
     206             :   std::string norm_str;
     207           4 :   if( dens.length()>0 ) {
     208           0 :     norm_str = " FUNC=-2*pi*x*" + dens;
     209             :   } else {
     210           8 :     norm_str = "," + ref_name + "_vol FUNC=-(2*pi*x/y)*" + shape_str[1];
     211             :   }
     212          12 :   readInputLine( getShortcutLabel() + ": CUSTOM PERIODIC=NO ARG=" + getShortcutLabel() + "_int" + norm_str );
     213          20 : }
     214             : 
     215             : }
     216             : }