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 6 : if( calc_dens ) action->readInputLine( lab + "_vol: VOLUME" );
44 3 : }
45 :
46 267 : void RDF::registerKeywords( Keywords& keys ) {
47 267 : ActionShortcut::registerKeywords( keys );
48 534 : keys.add("atoms","GROUP","");
49 534 : keys.add("atoms-2","GROUPA","");
50 534 : keys.add("atoms-2","GROUPB","");
51 534 : keys.add("compulsory","GRID_BIN","the number of bins to use when computing the RDF");
52 534 : keys.add("compulsory","KERNEL","GAUSSIAN","the type of kernel to use for computing the histograms for the RDF");
53 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");
54 534 : keys.add("compulsory","MAXR","the maximum distance to use for the rdf");
55 534 : keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density esimtation");
56 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");
57 534 : keys.add("compulsory","STRIDE","1","the frequency with which to compute the rdf and accumulate averages");
58 534 : keys.add("optional","DENSITY","the reference density to use when normalizing the RDF");
59 534 : keys.add("hidden","REFERENCE","this is the label of the reference objects");
60 534 : keys.setValueDescription("grid","the radial distribution function");
61 801 : keys.needsAction("REFERENCE_GRID"); keys.needsAction("VOLUME"); keys.needsAction("DISTANCE_MATRIX");
62 801 : keys.needsAction("CUSTOM"); keys.needsAction("KDE"); keys.needsAction("ACCUMULATE");
63 267 : keys.needsAction("CONSTANT");
64 267 : }
65 :
66 66 : RDF::RDF(const ActionOptions&ao):
67 : Action(ao),
68 66 : ActionShortcut(ao)
69 : {
70 : // Read in grid extent and number of bins
71 198 : std::string maxr, nbins, dens; parse("MAXR",maxr); parse("GRID_BIN",nbins); parse("DENSITY",dens);
72 132 : std::string grid_setup = "GRID_MIN=0 GRID_MAX=" + maxr + " GRID_BIN=" + nbins;
73 : // Create grid with normalizing function on it
74 132 : std::string refstr; parse("REFERENCE",refstr);
75 66 : if( refstr.length()==0 ) {
76 2 : createX2ReferenceObject( getShortcutLabel(), grid_setup, dens.length()==0, this ); refstr = getShortcutLabel();
77 : }
78 : // Read input to histogram
79 132 : std::string cutoff, kernel, bandwidth, kernel_data; parse("KERNEL",kernel);
80 66 : if( kernel=="DISCRETE" ) {
81 : cutoff = maxr; kernel_data="KERNEL=DISCRETE";
82 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");
83 : } else {
84 260 : parse("BANDWIDTH",bandwidth); double rcut; parse("CUTOFF",rcut); kernel_data="KERNEL=" + kernel + " IGNORE_IF_OUT_OF_RANGE BANDWIDTH=" + bandwidth;
85 65 : double bw; Tools::convert( bandwidth, bw ); double fcut; Tools::convert( maxr, fcut ); Tools::convert( fcut + sqrt(2.0*rcut)*bw, cutoff );
86 : }
87 :
88 : // Create contact matrix
89 132 : std::string natoms, str_norm_atoms, atom_str, group_str, groupa_str, groupb_str; parse("GROUP",group_str);
90 66 : if( group_str.length()>0 ) {
91 4 : atom_str="GROUP=" + group_str; std::vector<std::string> awords=Tools::getWords(group_str,"\t\n ,");
92 2 : Tools::interpretRanges( awords ); Tools::convert( awords.size(), natoms ); str_norm_atoms = natoms;
93 2 : } else {
94 128 : parse("GROUPA",groupa_str); parse("GROUPB",groupb_str);
95 64 : std::vector<std::string> awords=Tools::getWords(groupb_str,"\t\n ,");
96 64 : Tools::interpretRanges( awords ); Tools::convert( awords.size(), natoms );
97 128 : atom_str="GROUPA=" + groupa_str + " GROUPB=" + groupb_str;
98 64 : std::vector<std::string> bwords=Tools::getWords(groupa_str,"\t\n ,"); Tools::interpretRanges( bwords );
99 64 : Tools::convert( bwords.size()+1, str_norm_atoms );
100 64 : }
101 : // Retrieve the number of atoms
102 132 : readInputLine( getShortcutLabel() + "_mat: DISTANCE_MATRIX CUTOFF=" + cutoff + " " + atom_str);
103 :
104 : // Calculate weights of distances
105 132 : readInputLine( getShortcutLabel() + "_wmat: CUSTOM ARG=" + getShortcutLabel() + "_mat FUNC=step(" + cutoff + "-x) PERIODIC=NO");
106 : // Now create a histogram from the contact matrix
107 132 : unsigned clear, stride; parse("CLEAR",clear); parse("STRIDE",stride);
108 66 : if( clear==1 ) {
109 130 : readInputLine( getShortcutLabel() + "_kde: KDE ARG=" + getShortcutLabel() + "_mat VOLUMES=" + getShortcutLabel() + "_wmat " + grid_setup + " " + kernel_data);
110 : } else {
111 1 : std::string stridestr, clearstr; Tools::convert( stride, stridestr ); Tools::convert( clear, clearstr );
112 2 : readInputLine( getShortcutLabel() + "_okde: KDE ARG=" + getShortcutLabel() + "_mat HEIGHTS=" + getShortcutLabel() + "_wmat " + grid_setup + " " + kernel_data);
113 2 : readInputLine( getShortcutLabel() + "_kde: ACCUMULATE ARG=" + getShortcutLabel() + "_okde STRIDE=" + stridestr + " CLEAR=" + clearstr );
114 1 : readInputLine( getShortcutLabel() + "_one: CONSTANT VALUE=1");
115 2 : readInputLine( getShortcutLabel() + "_norm: ACCUMULATE ARG=" + getShortcutLabel() + "_one STRIDE=" + stridestr + " CLEAR=" + clearstr );
116 : }
117 : // Transform the histogram by normalizing factor for rdf
118 132 : readInputLine( getShortcutLabel() + "_vrdf: CUSTOM ARG=" + getShortcutLabel() + "_kde," + refstr + "_x2 FUNC=x/(4*pi*y) PERIODIC=NO");
119 : // And normalize by density and number of atoms (separated from above to avoid nans)
120 132 : std::string func_str = "PERIODIC=NO ARG=" + getShortcutLabel() + "_vrdf";
121 66 : if( dens.length()>0 ) {
122 0 : if( clear==1 ) func_str += " FUNC=x/(" + dens + "*" + str_norm_atoms + ")";
123 0 : else func_str += "," + getShortcutLabel() + "_norm FUNC=x/(y*" + dens + "*" + str_norm_atoms + ")";
124 : } else {
125 131 : if( clear==1 ) func_str += "," + refstr + "_vol FUNC=x*y/(" + natoms + "*" + str_norm_atoms + ")";
126 2 : else func_str += "," + refstr + "_vol," + getShortcutLabel() + "_norm FUNC=x*y/(z*" + natoms + "*" + str_norm_atoms + ")";
127 : }
128 132 : readInputLine( getShortcutLabel() + ": CUSTOM " + func_str);
129 66 : }
130 :
131 : }
132 : }
|
