Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2023 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 "core/ActionRegister.h"
23 : #include "core/ActionShortcut.h"
24 :
25 : //+PLUMEDOC ANALYSIS KL_ENTROPY
26 : /*
27 : Calculate the KL entropy of a distribution
28 :
29 : This shortcut was implemented in order to make the implementation of CVs like those described in the paper that is cited below straightforward.
30 : An example input that can be used to implement one of the CVs that is introduced in that paper is shown below:
31 :
32 : ```plumed
33 : #SETTINGS INPUTFILES=regtest/gridtools/rt-weights-integral/kde.grid
34 :
35 : # These two commands compute a set of weights that allow us to determine if each distance in our system is a bond
36 : d1: DISTANCE ATOMS1=1,2 ATOMS2=1,3 ATOMS3=1,4 ATOMS4=1,5 ATOMS5=2,3 ATOMS6=2,4 ATOMS7=2,5 ATOMS8=3,4 ATOMS9=3,5 ATOMS10=4,5
37 : d1lt: LESS_THAN ARG=d1 SWITCH={RATIONAL D_0=2.0 R_0=0.5 D_MAX=5.0}
38 :
39 : # These commands caluclate the angle between the z axis and the bond directions
40 : d1c: DISTANCE ATOMS1=2,1 ATOMS2=3,1 ATOMS3=4,1 ATOMS4=5,1 ATOMS5=3,2 ATOMS6=4,2 ATOMS7=5,2 ATOMS8=4,3 ATOMS9=5,3 ATOMS10=5,4 COMPONENTS
41 : d2: COMBINE ARG=d1c.x,d1c.y,d1c.z POWERS=2,2,2 PERIODIC=NO
42 : aa: MATHEVAL ARG=d1c.z,d2 FUNC=acos(x/sqrt(y)) PERIODIC=NO
43 :
44 : # These commands compute the torsion angle between the bond direction and the positive x direction around the z axis
45 : dd0: FIXEDATOM AT=0,0,0
46 : ddx: FIXEDATOM AT=1,0,0
47 : ddz: FIXEDATOM AT=0,0,1
48 :
49 : tt: TORSION ...
50 : VECTORA1=2,1 VECTORB1=ddx,dd0 AXIS1=ddz,dd0
51 : VECTORA2=3,1 VECTORB2=ddx,dd0 AXIS2=ddz,dd0
52 : VECTORA3=4,1 VECTORB3=ddx,dd0 AXIS3=ddz,dd0
53 : VECTORA4=5,1 VECTORB4=ddx,dd0 AXIS4=ddz,dd0
54 : VECTORA5=3,2 VECTORB5=ddx,dd0 AXIS5=ddz,dd0
55 : VECTORA6=4,2 VECTORB6=ddx,dd0 AXIS6=ddz,dd0
56 : VECTORA7=5,2 VECTORB7=ddx,dd0 AXIS7=ddz,dd0
57 : VECTORA8=4,3 VECTORB8=ddx,dd0 AXIS8=ddz,dd0
58 : VECTORA9=5,3 VECTORB9=ddx,dd0 AXIS9=ddz,dd0
59 : VECTORA10=5,4 VECTORB10=ddx,dd0 AXIS10=ddz,dd0
60 : ...
61 :
62 : # We now compute our instaneous normalised histogram of the bond directions. Notice that the weights are used here so that we do not consider
63 : # non-bonded atoms
64 : hu: KDE VOLUMES=d1lt ARG=aa,tt GRID_BIN=20,20 GRID_MIN=0,-pi GRID_MAX=pi,pi BANDWIDTH=0.2,0.2
65 : de: SUM ARG=d1lt PERIODIC=NO
66 : h: CUSTOM ARG=hu,de FUNC=x/y PERIODIC=NO
67 :
68 : # And lastly compute the KL_ENTROPY
69 : kl: KL_ENTROPY ARG=h REFERENCE=regtest/gridtools/rt-weights-integral/kde.grid VALUE=h
70 : PRINT ARG=kl FILE=colvar
71 : ```
72 :
73 : The CV `kl` that is defined in the input above describes the disribution of bond directions that is observed in the structure. To compute it we first define the orientations
74 : of all the bonds in our structure relative to the lab frame by computing the angle between each bond direction and the z axis and the torsional angle around the z axis between
75 : our bond vector and the positive x direction. In other words, we compute two of the three [spherical polar coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system)
76 : for each of the bond directions in our crystal structure.
77 :
78 : A normalised histogram that describes the instaneous distribution for these bond directions is then computed by using the [KDE](KDE.md) action. The final scalar value for the CV is computed by
79 : determining the [Kullbeck-Leibler divergence](https://en.wikipedia.org/wiki/Kullback–Leibler_divergence) between this instaneous distribution and a reference distribution that is provided in input.
80 :
81 : Notice that this action will also works if you use KDEs computed using the [SPHERICAL_KDE](SPHERICAL_KDE.md) action in input as is illustrated in the following input:
82 :
83 : ```plumed
84 : #SETTINGS INPUTFILES=regtest/gridtools/rt-kldiv/averageSp_X240.grid
85 :
86 : # These two commands compute a set of weights that allow us to determine if each distance in our system is a bond
87 : d1: DISTANCE ATOMS1=1,2 ATOMS2=1,3 ATOMS3=1,4 ATOMS4=1,5 ATOMS5=2,3 ATOMS6=2,4 ATOMS7=2,5 ATOMS8=3,4 ATOMS9=3,5 ATOMS10=4,5
88 : d1lt: LESS_THAN ARG=d1 SWITCH={RATIONAL D_0=2.0 R_0=0.5 D_MAX=5.0}
89 :
90 : # These commands compute the directors of the bonds between the atoms
91 : d1c: DISTANCE ATOMS1=2,1 ATOMS2=3,1 ATOMS3=4,1 ATOMS4=5,1 ATOMS5=3,2 ATOMS6=4,2 ATOMS7=5,2 ATOMS8=4,3 ATOMS9=5,3 ATOMS10=5,4 COMPONENTS
92 : d1l: CUSTOM ARG=d1c.x,d1c.y,d1c.z FUNC=sqrt(x*x+y*y+z*z) PERIODIC=NO
93 : d1nx: CUSTOM ARG=d1c.x,d1l FUNC=x/y PERIODIC=NO
94 : d1ny: CUSTOM ARG=d1c.y,d1l FUNC=x/y PERIODIC=NO
95 : d1nz: CUSTOM ARG=d1c.z,d1l FUNC=x/y PERIODIC=NO
96 :
97 : # Now compute our instaneous normalised histogram of the bond directions
98 : hu: SPHERICAL_KDE HEIGHTS=d1lt ARG=d1nx,d1ny,d1nz GRID_BIN=400 CONCENTRATION=100.0
99 : de: SUM ARG=d1lt PERIODIC=NO
100 : h: CUSTOM ARG=hu,de FUNC=x/y PERIODIC=NO
101 :
102 : # And lastly compute the KL_ENTROPY
103 : kl: KL_ENTROPY ARG=h REFERENCE=regtest/gridtools/rt-kldiv/averageSp_X240.grid VALUE=av_Sp_X
104 : PRINT ARG=kl FILE=colvar
105 : ```
106 :
107 : */
108 : //+ENDPLUMEDOC
109 :
110 : namespace PLMD {
111 : namespace gridtools {
112 :
113 : class KLEntropy : public ActionShortcut {
114 : public:
115 : static void registerKeywords( Keywords& keys );
116 : explicit KLEntropy(const ActionOptions&ao);
117 : };
118 :
119 : PLUMED_REGISTER_ACTION(KLEntropy,"KL_ENTROPY")
120 :
121 7 : void KLEntropy::registerKeywords( Keywords& keys ) {
122 7 : ActionShortcut::registerKeywords(keys);
123 7 : keys.add("compulsory","ARG","the grid that you wish to use in the KL entropy calculation");
124 7 : keys.add("compulsory","REFERENCE","a file containing the reference density in grid format");
125 7 : keys.add("compulsory","VALUE","the name of the value that should be read from the grid");
126 14 : keys.setValueDescription("scalar","the value of the KL-Entropy");
127 7 : keys.addDOI("10.1021/acs.jctc.7b01027");
128 7 : keys.needsAction("REFERENCE_GRID");
129 7 : keys.needsAction("CUSTOM");
130 7 : keys.needsAction("INTEGRATE_GRID");
131 7 : }
132 :
133 5 : KLEntropy::KLEntropy( const ActionOptions& ao):
134 : Action(ao),
135 5 : ActionShortcut(ao) {
136 : // Reference grid object
137 : std::string ref_str, val_str, input_g;
138 5 : parse("VALUE",val_str);
139 5 : parse("REFERENCE",ref_str);
140 5 : parse("ARG",input_g);
141 10 : readInputLine( getShortcutLabel() + "_ref: REFERENCE_GRID VALUE=" + val_str + " FILE=" + ref_str );
142 : // Compute KL divergence
143 5 : if( input_g=="") {
144 0 : plumed_merror("could not find ARG keyword in input to KL_ENTROPY");
145 : }
146 10 : readInputLine( getShortcutLabel() + "_kl: CUSTOM ARG=" + input_g + "," + getShortcutLabel() + "_ref FUNC=y*log(y/(0.5*(x+y))) PERIODIC=NO");
147 : // Compute integral
148 10 : readInputLine( getShortcutLabel() + ": INTEGRATE_GRID ARG=" + getShortcutLabel() + "_kl PERIODIC=NO");
149 5 : }
150 :
151 : }
152 : }
|
