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/PlumedMain.h"
24 : #include "core/ActionSet.h"
25 : #include "core/ActionWithArguments.h"
26 : #include "core/ActionRegister.h"
27 :
28 : //+PLUMEDOC ANALYSIS RDF
29 : /*
30 : Calculate the radial distribution function
31 :
32 : This shortcut can be used to evaluate the [radial distribution function](https://en.wikipedia.org/wiki/Radial_distribution_function) for
33 : a collection of atoms. The following example illustrates how it can be used to compute the average radial distribution from a
34 : trajectory:
35 :
36 : ```plumed
37 : rdf: RDF GROUP=1-108 MAXR=2.5 GRID_BIN=25 KERNEL=DISCRETE STRIDE=1 CLEAR=0
38 : DUMPGRID ARG=rdf FILE=rdf.dat
39 : ```
40 :
41 : If you expand the shortcut in this input you can see how the radial distribution function is computed using [KDE](KDE.md) actions and that the
42 : average is obtained by using [ACCUMULATE](ACCUMULATE.md) actions. By using expanded versions of this shortcut you can thus calculate and print
43 : the instantaneous value of the radial distribution function. Alternatively, if you use the CLEAR option in the [ACCUMULATE](ACCUMULATE.md) commands
44 : you can calculate radial distribution functions from various parts of the trajectory.
45 :
46 : The following example shows how you can use this function to calculate the radial distribution for the atoms in GROUPB around the atoms in GROUPA.
47 : Notice that the DISCRETE keyword has been removed here so we are using [kernel density estimation](https://en.wikipedia.org/wiki/Kernel_density_estimation)
48 : to compute the radial distribution function.
49 :
50 : ```plumed
51 : rdf: RDF GROUPA=1-108 GROUPB=109-300 MAXR=2.5 GRID_BIN=25 BANDWIDTH=0.01 DENSITY=1 STRIDE=1 CLEAR=0
52 : DUMPGRID ARG=rdf FILE=rdf.dat
53 : ```
54 :
55 : 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}$.
56 : 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
57 : first input above.
58 :
59 : Notice that you can also use the `NO_AVERAGE` keyword as shown below to calculate an independent RDF around each of the atoms:
60 :
61 : ```plumed
62 : rdf: RDF GROUPA=1-3 GROUPB=4-108 MAXR=2.5 GRID_BIN=25 BANDWIDTH=0.01 DENSITY=1 STRIDE=1 NO_AVERAGE CLEAR=0
63 : DUMPGRID ARG=rdf FILE=rdf.dat
64 : ```
65 :
66 : This command will output a function on a $25 \times 3$ grid. The first coordinate is the distance between atoms. The second coordinate is then then an index that tells
67 : you, which central atom. In other words, the input above will output three radial distribution functions. The first of these is the radial distribution function around atom 1,
68 : the second is the radial distribution function around atom 2 and the last is the radial distribution function around atom 3.
69 :
70 : */
71 : //+ENDPLUMEDOC
72 :
73 : namespace PLMD {
74 : namespace gridtools {
75 :
76 : PLUMED_REGISTER_ACTION(RDF,"RDF")
77 :
78 13 : void RDF::getDistanceMatrixShape( const std::string& lab, ActionShortcut* action, std::vector<std::string>& shape_str ) {
79 13 : std::vector<Value*> args;
80 13 : std::vector<std::string> arglabs(1);
81 : arglabs[0] = lab;
82 13 : ActionWithArguments::interpretArgumentList( arglabs, action->plumed.getActionSet(), action, args );
83 13 : Tools::convert( args[0]->getShape()[0], shape_str[0] );
84 13 : Tools::convert( args[0]->getShape()[1], shape_str[1] );
85 26 : }
86 :
87 0 : void RDF::createX2ReferenceObject( const std::string& lab, const std::string& grid_setup, const bool& calc_dens, const bool& no_average, ActionShortcut* action ) {
88 : // Create grid with normalizing function
89 0 : if( no_average ) {
90 0 : action->readInputLine( lab + "_x2: REFERENCE_GRID PERIODIC=NO,NO FUNC=x*x+0*y " + grid_setup );
91 : } else {
92 0 : action->readInputLine( lab + "_x2: REFERENCE_GRID PERIODIC=NO FUNC=x*x " + grid_setup );
93 : }
94 : // Compute density if required
95 0 : if( calc_dens ) {
96 0 : action->readInputLine( lab + "_vol: VOLUME" );
97 : }
98 0 : }
99 :
100 32 : void RDF::registerKeywords( Keywords& keys ) {
101 32 : ActionShortcut::registerKeywords( keys );
102 32 : keys.add("atoms","GROUP","the atoms that are being used to calculate the RDF");
103 32 : keys.addDeprecatedKeyword("ATOMS","GROUP");
104 32 : keys.add("atoms-2","GROUPA","the atoms that you would like to compute the RDF about. Must be used with GROUPB.");
105 32 : keys.add("atoms-2","GROUPB","the atoms that you would like to use when computing the RDF around the atoms that were specified with GROUPA");
106 32 : keys.add("compulsory","GRID_BIN","the number of bins to use when computing the RDF");
107 : std::string strcut;
108 32 : Tools::convert( sqrt(2*dp2cutoff), strcut );
109 32 : keys.add("compulsory","KERNEL","GAUSSIAN D_MAX=" + strcut,"the type of kernel to use for computing the histograms for the RDF");
110 32 : keys.add("compulsory","MAXR","the maximum distance to use for the rdf");
111 32 : keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density esimtation");
112 32 : if( keys.getDisplayName()=="RDF") {
113 19 : keys.add("compulsory","CLEAR","1","the frequency with which to clear the estimate of the rdf. Set equal to 0 if you want to compute an rdf over the whole trajectory");
114 19 : keys.add("compulsory","STRIDE","1","the frequency with which to compute the rdf and accumulate averages");
115 : }
116 32 : keys.add("optional","DENSITY","the reference density to use when normalizing the RDF");
117 32 : keys.addFlag("NO_AVERAGE",false,"output a two-dimensional grid with separate RDF functions for each of the central atoms");
118 32 : keys.add("hidden","REFERENCE","this is the label of the reference objects");
119 64 : keys.setValueDescription("grid","the radial distribution function");
120 32 : keys.needsAction("REFERENCE_GRID");
121 32 : keys.needsAction("VOLUME");
122 32 : keys.needsAction("DISTANCE_MATRIX");
123 32 : keys.needsAction("CUSTOM");
124 32 : keys.needsAction("KDE");
125 32 : keys.needsAction("ACCUMULATE");
126 32 : keys.needsAction("CONSTANT");
127 32 : keys.needsAction("ONES");
128 64 : keys.needsAction("OUTER_PRODUCT");
129 64 : }
130 :
131 10 : RDF::RDF(const ActionOptions&ao):
132 : Action(ao),
133 10 : ActionShortcut(ao) {
134 : bool noaverage;
135 20 : parseFlag("NO_AVERAGE",noaverage);
136 : // Read in grid extent and number of bins
137 : std::string maxr, nbins, dens;
138 10 : parse("MAXR",maxr);
139 10 : parse("GRID_BIN",nbins);
140 20 : parse("DENSITY",dens);
141 : double dval;
142 10 : if( dens.length()>0 && Tools::convertNoexcept( dens, dval ) && dval<=0 ) {
143 2 : error("DENSITY should be positive");
144 : }
145 : // Read input to histogram
146 : std::string cutoff, kernel, bandwidth, kernel_data;
147 9 : parse("KERNEL",kernel);
148 9 : if( kernel=="DISCRETE" ) {
149 : cutoff = maxr;
150 : kernel_data="KERNEL=DISCRETE";
151 4 : warning("rdf is normalised by dividing by the surface area at the grid value and not by the volume of the bin as it should be with discrete kernels");
152 : } else {
153 7 : parse("BANDWIDTH",bandwidth);
154 : double bw;
155 7 : Tools::convert( bandwidth, bw );
156 7 : if( noaverage ) {
157 : bandwidth += ",0";
158 : }
159 14 : kernel_data="KERNEL={" + kernel + "} BANDWIDTH=" + bandwidth;
160 : double fcut;
161 7 : Tools::convert( maxr, fcut );
162 7 : double rcut = 0;
163 7 : std::vector<std::string> kwords( Tools::getWords(kernel) );
164 7 : Tools::parse( kwords, "D_MAX", rcut );
165 7 : if( rcut==0 ) {
166 1 : rcut = sqrt(2.0*dp2cutoff);
167 : }
168 7 : Tools::convert( fcut + rcut*bw, cutoff );
169 7 : }
170 :
171 : // Create contact matrix
172 : std::string atom_str, group_str, groupa_str, groupb_str;
173 18 : parse("GROUP",group_str);
174 9 : if( group_str.length()>0 ) {
175 16 : atom_str="GROUP=" + group_str;
176 : } else {
177 1 : parse("GROUPA",groupa_str);
178 1 : parse("GROUPB",groupb_str);
179 3 : atom_str="GROUPA=" + groupa_str + " GROUPB=" + groupb_str;
180 : }
181 :
182 : // Calculate all the distances
183 18 : readInputLine( getShortcutLabel() + "_mat: DISTANCE_MATRIX LINKCELL_CUTOFF=" + cutoff + " " + atom_str);
184 : // Get the shape of the matrix
185 9 : std::vector<std::string> shape_str(2);
186 9 : getDistanceMatrixShape( getShortcutLabel() + "_mat", this, shape_str );
187 :
188 : // Setup the grid
189 18 : std::string grid_setup = "GRID_MIN=0 GRID_MAX=" + maxr + " GRID_BIN=" + nbins;
190 9 : if( noaverage ) {
191 : int iatoms;
192 : std::string num, str_values;
193 2 : Tools::convert( 0, str_values );
194 2 : Tools::convert( shape_str[0], iatoms );
195 67 : for(int i=1; i<iatoms; ++i) {
196 65 : Tools::convert( i, num );
197 130 : str_values += "," + num;
198 : }
199 2 : Tools::convert( iatoms-1, num );
200 2 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + shape_str[1] );
201 2 : readInputLine( getShortcutLabel() + "_index: CONSTANT VALUES=" + str_values );
202 6 : readInputLine( getShortcutLabel() + "_atno: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_index," + getShortcutLabel() + "_ones" );
203 10 : grid_setup = "GRID_MIN=0,0 GRID_MAX=" + maxr + "," + num + " GRID_BIN=" + nbins + "," + num;
204 : shape_str[0] = "1";
205 : }
206 : // Create grid with normalizing function on it
207 : std::string refstr;
208 18 : parse("REFERENCE",refstr);
209 9 : if( refstr.length()==0 ) {
210 9 : if( noaverage ) {
211 4 : readInputLine( getShortcutLabel() + "_x2: REFERENCE_GRID PERIODIC=NO,NO FUNC=x*x+0*y " + grid_setup );
212 : } else {
213 14 : readInputLine( getShortcutLabel() + "_x2: REFERENCE_GRID PERIODIC=NO FUNC=x*x " + grid_setup );
214 : }
215 9 : refstr = getShortcutLabel();
216 : }
217 9 : if( dens.length()==0 ) {
218 18 : readInputLine( getShortcutLabel() + "_vol: VOLUME" );
219 : }
220 : // Calculate weights of distances
221 27 : readInputLine( getShortcutLabel() + "_wmat: CUSTOM ARG=" + getShortcutLabel() + "_mat FUNC=step(" + cutoff + "-x) PERIODIC=NO");
222 : // Now create a histogram from the contact matrix
223 : unsigned clear, stride;
224 9 : parse("CLEAR",clear);
225 9 : parse("STRIDE",stride);
226 9 : if( clear==1 ) {
227 6 : if( noaverage ) {
228 6 : readInputLine( getShortcutLabel() + "_kde: KDE ARG=" + getShortcutLabel() + "_mat," + getShortcutLabel() + "_atno VOLUMES=" + getShortcutLabel() + "_wmat " + grid_setup + " " + kernel_data);
229 : } else {
230 25 : readInputLine( getShortcutLabel() + "_kde: KDE ARG=" + getShortcutLabel() + "_mat VOLUMES=" + getShortcutLabel() + "_wmat " + grid_setup + " " + kernel_data);
231 : }
232 : } else {
233 : std::string stridestr, clearstr;
234 3 : Tools::convert( stride, stridestr );
235 3 : Tools::convert( clear, clearstr );
236 3 : if( noaverage ) {
237 6 : readInputLine( getShortcutLabel() + "_okde: KDE ARG=" + getShortcutLabel() + "_mat," + getShortcutLabel() + "_atno HEIGHTS=" + getShortcutLabel() + "_wmat " + grid_setup + " " + kernel_data);
238 : } else {
239 10 : readInputLine( getShortcutLabel() + "_okde: KDE ARG=" + getShortcutLabel() + "_mat HEIGHTS=" + getShortcutLabel() + "_wmat " + grid_setup + " " + kernel_data);
240 : }
241 9 : readInputLine( getShortcutLabel() + "_kde: ACCUMULATE ARG=" + getShortcutLabel() + "_okde STRIDE=" + stridestr + " CLEAR=" + clearstr );
242 3 : readInputLine( getShortcutLabel() + "_one: CONSTANT VALUE=1");
243 12 : readInputLine( getShortcutLabel() + "_norm: ACCUMULATE ARG=" + getShortcutLabel() + "_one STRIDE=" + stridestr + " CLEAR=" + clearstr );
244 : }
245 : // Transform the histogram by normalizing factor for rdf
246 27 : readInputLine( getShortcutLabel() + "_vrdf: CUSTOM ARG=" + getShortcutLabel() + "_kde," + refstr + "_x2 FUNC=x/(4*pi*y) PERIODIC=NO");
247 : // And normalize by density and number of atoms (separated from above to avoid nans)
248 18 : std::string func_str = "PERIODIC=NO ARG=" + getShortcutLabel() + "_vrdf";
249 9 : if( dens.length()>0 ) {
250 0 : if( clear==1 ) {
251 0 : func_str += " FUNC=x/(" + dens + "*" + shape_str[0] + ")";
252 : } else {
253 0 : func_str += "," + getShortcutLabel() + "_norm FUNC=x/(y*" + dens + "*" + shape_str[0] + ")";
254 : }
255 : } else {
256 9 : if( clear==1 ) {
257 24 : func_str += "," + refstr + "_vol FUNC=x*y/(" + shape_str[1] + "*" + shape_str[0] + ")";
258 : } else {
259 15 : func_str += "," + refstr + "_vol," + getShortcutLabel() + "_norm FUNC=x*y/(z*" + shape_str[1] + "*" + shape_str[0] + ")";
260 : }
261 : }
262 18 : readInputLine( getShortcutLabel() + ": CUSTOM " + func_str);
263 20 : }
264 :
265 : }
266 : }
|
