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 :
25 : //+PLUMEDOC ANALYSIS RDF
26 : /*
27 : Calculate the radial distribution function
28 :
29 : \par Examples
30 :
31 : */
32 : //+ENDPLUMEDOC
33 :
34 : namespace PLMD {
35 : namespace gridtools {
36 :
37 : PLUMED_REGISTER_ACTION(RDF,"RDF")
38 :
39 3 : void RDF::createX2ReferenceObject( const std::string& lab, const std::string& grid_setup, const bool& calc_dens, ActionShortcut* action ) {
40 : // Create grid with normalizing function
41 6 : action->readInputLine( lab + "_x2: REFERENCE_GRID PERIODIC=NO FUNC=x*x " + grid_setup );
42 : // Compute density if required
43 3 : if( calc_dens ) {
44 6 : action->readInputLine( lab + "_vol: VOLUME" );
45 : }
46 3 : }
47 :
48 267 : void RDF::registerKeywords( Keywords& keys ) {
49 267 : ActionShortcut::registerKeywords( keys );
50 534 : keys.add("atoms","GROUP","");
51 534 : keys.add("atoms-2","GROUPA","");
52 534 : keys.add("atoms-2","GROUPB","");
53 534 : keys.add("compulsory","GRID_BIN","the number of bins to use when computing the RDF");
54 534 : keys.add("compulsory","KERNEL","GAUSSIAN","the type of kernel to use for computing the histograms for the RDF");
55 534 : keys.add("compulsory","CUTOFF","6.25","the cutoff at which to stop evaluating the kernel functions is set equal to sqrt(2*x)*bandwidth in each direction where x is this number");
56 534 : keys.add("compulsory","MAXR","the maximum distance to use for the rdf");
57 534 : keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density esimtation");
58 534 : 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");
59 534 : keys.add("compulsory","STRIDE","1","the frequency with which to compute the rdf and accumulate averages");
60 534 : keys.add("optional","DENSITY","the reference density to use when normalizing the RDF");
61 534 : keys.add("hidden","REFERENCE","this is the label of the reference objects");
62 267 : keys.setValueDescription("the radial distribution function");
63 267 : keys.needsAction("REFERENCE_GRID");
64 267 : keys.needsAction("VOLUME");
65 267 : keys.needsAction("DISTANCE_MATRIX");
66 267 : keys.needsAction("CUSTOM");
67 267 : keys.needsAction("KDE");
68 267 : keys.needsAction("ACCUMULATE");
69 267 : keys.needsAction("CONSTANT");
70 267 : }
71 :
72 66 : RDF::RDF(const ActionOptions&ao):
73 : Action(ao),
74 66 : ActionShortcut(ao) {
75 : // Read in grid extent and number of bins
76 : std::string maxr, nbins, dens;
77 66 : parse("MAXR",maxr);
78 66 : parse("GRID_BIN",nbins);
79 66 : parse("DENSITY",dens);
80 132 : std::string grid_setup = "GRID_MIN=0 GRID_MAX=" + maxr + " GRID_BIN=" + nbins;
81 : // Create grid with normalizing function on it
82 : std::string refstr;
83 132 : parse("REFERENCE",refstr);
84 66 : if( refstr.length()==0 ) {
85 2 : createX2ReferenceObject( getShortcutLabel(), grid_setup, dens.length()==0, this );
86 2 : refstr = getShortcutLabel();
87 : }
88 : // Read input to histogram
89 : std::string cutoff, kernel, bandwidth, kernel_data;
90 132 : parse("KERNEL",kernel);
91 66 : if( kernel=="DISCRETE" ) {
92 : cutoff = maxr;
93 : kernel_data="KERNEL=DISCRETE";
94 2 : 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");
95 : } else {
96 65 : parse("BANDWIDTH",bandwidth);
97 : double rcut;
98 65 : parse("CUTOFF",rcut);
99 130 : kernel_data="KERNEL=" + kernel + " IGNORE_IF_OUT_OF_RANGE BANDWIDTH=" + bandwidth;
100 : double bw;
101 65 : Tools::convert( bandwidth, bw );
102 : double fcut;
103 65 : Tools::convert( maxr, fcut );
104 65 : Tools::convert( fcut + sqrt(2.0*rcut)*bw, cutoff );
105 : }
106 :
107 : // Create contact matrix
108 : std::string natoms, str_norm_atoms, atom_str, group_str, groupa_str, groupb_str;
109 132 : parse("GROUP",group_str);
110 66 : if( group_str.length()>0 ) {
111 2 : atom_str="GROUP=" + group_str;
112 2 : std::vector<std::string> awords=Tools::getWords(group_str,"\t\n ,");
113 2 : Tools::interpretRanges( awords );
114 2 : Tools::convert( awords.size(), natoms );
115 : str_norm_atoms = natoms;
116 2 : } else {
117 64 : parse("GROUPA",groupa_str);
118 64 : parse("GROUPB",groupb_str);
119 64 : std::vector<std::string> awords=Tools::getWords(groupb_str,"\t\n ,");
120 64 : Tools::interpretRanges( awords );
121 64 : Tools::convert( awords.size(), natoms );
122 128 : atom_str="GROUPA=" + groupa_str + " GROUPB=" + groupb_str;
123 64 : std::vector<std::string> bwords=Tools::getWords(groupa_str,"\t\n ,");
124 64 : Tools::interpretRanges( bwords );
125 64 : Tools::convert( bwords.size()+1, str_norm_atoms );
126 64 : }
127 : // Retrieve the number of atoms
128 132 : readInputLine( getShortcutLabel() + "_mat: DISTANCE_MATRIX CUTOFF=" + cutoff + " " + atom_str);
129 :
130 : // Calculate weights of distances
131 132 : readInputLine( getShortcutLabel() + "_wmat: CUSTOM ARG=" + getShortcutLabel() + "_mat FUNC=step(" + cutoff + "-x) PERIODIC=NO");
132 : // Now create a histogram from the contact matrix
133 : unsigned clear, stride;
134 66 : parse("CLEAR",clear);
135 66 : parse("STRIDE",stride);
136 66 : if( clear==1 ) {
137 130 : readInputLine( getShortcutLabel() + "_kde: KDE ARG=" + getShortcutLabel() + "_mat VOLUMES=" + getShortcutLabel() + "_wmat " + grid_setup + " " + kernel_data);
138 : } else {
139 : std::string stridestr, clearstr;
140 1 : Tools::convert( stride, stridestr );
141 1 : Tools::convert( clear, clearstr );
142 2 : readInputLine( getShortcutLabel() + "_okde: KDE ARG=" + getShortcutLabel() + "_mat HEIGHTS=" + getShortcutLabel() + "_wmat " + grid_setup + " " + kernel_data);
143 2 : readInputLine( getShortcutLabel() + "_kde: ACCUMULATE ARG=" + getShortcutLabel() + "_okde STRIDE=" + stridestr + " CLEAR=" + clearstr );
144 1 : readInputLine( getShortcutLabel() + "_one: CONSTANT VALUE=1");
145 2 : readInputLine( getShortcutLabel() + "_norm: ACCUMULATE ARG=" + getShortcutLabel() + "_one STRIDE=" + stridestr + " CLEAR=" + clearstr );
146 : }
147 : // Transform the histogram by normalizing factor for rdf
148 132 : readInputLine( getShortcutLabel() + "_vrdf: CUSTOM ARG=" + getShortcutLabel() + "_kde," + refstr + "_x2 FUNC=x/(4*pi*y) PERIODIC=NO");
149 : // And normalize by density and number of atoms (separated from above to avoid nans)
150 132 : std::string func_str = "PERIODIC=NO ARG=" + getShortcutLabel() + "_vrdf";
151 66 : if( dens.length()>0 ) {
152 0 : if( clear==1 ) {
153 0 : func_str += " FUNC=x/(" + dens + "*" + str_norm_atoms + ")";
154 : } else {
155 0 : func_str += "," + getShortcutLabel() + "_norm FUNC=x/(y*" + dens + "*" + str_norm_atoms + ")";
156 : }
157 : } else {
158 66 : if( clear==1 ) {
159 130 : func_str += "," + refstr + "_vol FUNC=x*y/(" + natoms + "*" + str_norm_atoms + ")";
160 : } else {
161 2 : func_str += "," + refstr + "_vol," + getShortcutLabel() + "_norm FUNC=x*y/(z*" + natoms + "*" + str_norm_atoms + ")";
162 : }
163 : }
164 132 : readInputLine( getShortcutLabel() + ": CUSTOM " + func_str);
165 66 : }
166 :
167 : }
168 : }
|
