Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 : namespace PLMD {
26 : namespace gridtools {
27 :
28 : //+PLUMEDOC GRIDCALC MULTICOLVARDENS
29 : /*
30 : Evaluate the average value of a multicolvar on a grid.
31 :
32 : This keyword allows one to construct a phase field representation for a symmetry function from
33 : an atomistic description. If each atom has an associated order parameter, \f$\phi_i\f$ then a
34 : smooth phase field function \f$\phi(r)\f$ can be computed using:
35 :
36 : \f[
37 : \phi(\mathbf{r}) = \frac{\sum_i K(\mathbf{r}-\mathbf{r}_i) \phi_i }{ \sum_i K(\mathbf{r} - \mathbf{r}_i )}
38 : \f]
39 :
40 : where \f$\mathbf{r}_i\f$ is the position of atom \f$i\f$, the sums run over all the atoms input
41 : and \f$K(\mathbf{r} - \mathbf{r}_i)\f$ is one of the \ref kernelfunctions implemented in plumed.
42 : This action calculates the above function on a grid, which can then be used in the input to further
43 : actions.
44 :
45 : \par Examples
46 :
47 : The following example shows perhaps the simplest way in which this action can be used. The following
48 : input computes the density of atoms at each point on the grid and outputs this quantity to a file. In
49 : other words this input instructs plumed to calculate \f$\rho(\mathbf{r}) = \sum_i K(\mathbf{r} - \mathbf{r}_i )\f$
50 :
51 : \plumedfile
52 : dens: DENSITY SPECIES=1-100
53 : grid: MULTICOLVARDENS DATA=dens ORIGIN=1 DIR=xyz NBINS=100,100,100 BANDWIDTH=0.05,0.05,0.05 STRIDE=1
54 : DUMPGRID GRID=grid STRIDE=500 FILE=density
55 : \endplumedfile
56 :
57 : In the above example density is added to the grid on every step. The PRINT_GRID instruction thus tells PLUMED to
58 : output the average density at each point on the grid every 500 steps of simulation. Notice that the that grid output
59 : on step 1000 is an average over all 1000 frames of the trajectory. If you would like to analyze these two blocks
60 : of data separately you must use the CLEAR flag.
61 :
62 : This second example computes an order parameter (in this case \ref FCCUBIC) and constructs a phase field model
63 : for this order parameter using the equation above.
64 :
65 : \plumedfile
66 : fcc: FCCUBIC SPECIES=1-5184 SWITCH={CUBIC D_0=1.2 D_MAX=1.5} ALPHA=27
67 : dens: MULTICOLVARDENS DATA=fcc ORIGIN=1 DIR=xyz NBINS=14,14,28 BANDWIDTH=1.0,1.0,1.0 STRIDE=1 CLEAR=1
68 : DUMPCUBE GRID=dens STRIDE=1 FILE=dens.cube
69 : \endplumedfile
70 :
71 : In this example the phase field model is computed and output to a file on every step of the simulation. Furthermore,
72 : because the CLEAR=1 keyword is set on the MULTICOLVARDENS line each Gaussian cube file output is a phase field
73 : model for a particular trajectory frame. The average value accumulated thus far is cleared at the start of every single
74 : timestep and there is no averaging over trajectory frames in this case.
75 :
76 : */
77 : //+ENDPLUMEDOC
78 :
79 : class MultiColvarDensity : public ActionShortcut {
80 : public:
81 : explicit MultiColvarDensity(const ActionOptions&);
82 : static void registerKeywords( Keywords& keys );
83 : };
84 :
85 : PLUMED_REGISTER_ACTION(MultiColvarDensity,"MULTICOLVARDENS")
86 :
87 10 : void MultiColvarDensity::registerKeywords( Keywords& keys ) {
88 10 : ActionShortcut::registerKeywords( keys );
89 10 : keys.add("compulsory","STRIDE","1","the frequency with which to accumulate the densities");
90 10 : keys.add("compulsory","CLEAR","0","the frequency with which to clear the density");
91 10 : keys.add("compulsory","ORIGIN","we will use the position of this atom as the origin");
92 10 : keys.add("compulsory","DIR","the direction in which to calculate the density profile");
93 10 : keys.add("optional","BANDWIDTH","the bandwidths for kernel density esimtation");
94 10 : keys.add("compulsory","KERNEL","GAUSSIAN","the kernel function you are using. More details on the kernels available "
95 : "in plumed plumed can be found in \\ref kernelfunctions.");
96 10 : keys.add("optional","NBINS","the number of bins to use in each direction (alternative to GRID_NBIN)");
97 10 : keys.add("optional","GRID_MIN","the lower bounds for the grid (default boxlengths)");
98 10 : keys.add("optional","GRID_MAX","the upper bounds for the grid (default boxlengths)");
99 10 : keys.add("optional","DATA","the multicolvar which you would like to calculate the density profile for");
100 10 : keys.add("optional","ATOMS","if you are calculating a atomic density you use this keyword to specify the atoms that are involved");
101 10 : keys.addFlag("UNORMALIZED",false,"do not divide by the density");
102 10 : keys.add("optional","NORMALIZATION","set true/false to determine how to the data is normalised");
103 20 : keys.setValueDescription("grid","the average value of the order parameters at each point on the grid");
104 10 : keys.needsAction("DISTANCES");
105 10 : keys.needsAction("KDE");
106 10 : keys.needsAction("ACCUMULATE");
107 10 : keys.needsAction("CUSTOM");
108 10 : keys.needsAction("ONES");
109 10 : keys.needsAction("CUSTOM");
110 10 : }
111 :
112 8 : MultiColvarDensity::MultiColvarDensity(const ActionOptions&ao):
113 : Action(ao),
114 8 : ActionShortcut(ao) {
115 : // Read in the position of the origin
116 : std::string origin_str;
117 16 : parse("ORIGIN",origin_str);
118 : // Read in the quantity we are calculating the density for
119 : std::string atoms_str, data_str;
120 8 : parse("ATOMS",atoms_str);
121 16 : parse("DATA",data_str);
122 8 : if( atoms_str.length()==0 && data_str.length()==0 ) {
123 0 : error("quantity to calculate the density for was not specified used DATA/ATOMS");
124 : }
125 : // Get the information on the direction for the density
126 : std::string dir, direction_string;
127 8 : parse("DIR",dir);
128 8 : std::string nbins="";
129 16 : parse("NBINS",nbins);
130 8 : if(nbins.length()>0) {
131 0 : nbins=" GRID_BIN=" + nbins;
132 : }
133 8 : if( dir=="x" ) {
134 16 : direction_string = "ARG=" + getShortcutLabel() + "_dist.x " + nbins;
135 0 : } else if( dir=="y" ) {
136 0 : direction_string = "ARG=" + getShortcutLabel() + "_dist.y " + nbins;
137 0 : } else if( dir=="z" ) {
138 0 : direction_string = "ARG=" + getShortcutLabel() + "_dist.z " + nbins;
139 0 : } else if( dir=="xy" ) {
140 0 : direction_string = "ARG=" + getShortcutLabel() + "_dist.x," + getShortcutLabel() + "_dist.y " + nbins;
141 0 : } else if( dir=="xz" ) {
142 0 : direction_string = "ARG=" + getShortcutLabel() + "_dist.x," + getShortcutLabel() + "_dist.z " + nbins;
143 0 : } else if( dir=="yz" ) {
144 0 : direction_string = "ARG=" + getShortcutLabel() + "_dist.y," + getShortcutLabel() + "_dist.z " + nbins;
145 0 : } else if( dir=="xyz" ) {
146 0 : direction_string = "ARG=" + getShortcutLabel() + "_dist.x," + getShortcutLabel() + "_dist.y," + getShortcutLabel() + "_dist.z " + nbins;
147 : } else {
148 0 : error( dir + " is invalid dir specification use x/y/z/xy/xz/yz/xyz");
149 : }
150 :
151 : // Parse the keymap for this averaging stuff
152 : std::string stride, clear;
153 8 : parse("STRIDE",stride);
154 8 : parse("CLEAR",clear);
155 : bool unorm;
156 8 : parseFlag("UNORMALIZED",unorm);
157 8 : if( !unorm ) {
158 : std::string normstr;
159 12 : parse("NORMALIZATION",normstr);
160 6 : if( normstr=="false" ) {
161 0 : unorm=true;
162 : }
163 : }
164 : // Create distance action
165 : bool hasheights;
166 16 : std::string dist_words = getShortcutLabel() + "_dist: DISTANCES COMPONENTS ORIGIN=" + origin_str;
167 8 : if( atoms_str.length()>0 ) {
168 : hasheights=false;
169 6 : dist_words += " ATOMS=" + atoms_str;
170 : } else {
171 : hasheights=true;
172 10 : dist_words += " ATOMS=" + data_str;
173 : }
174 : // plumed_massert( keys.count("ORIGIN"), "you must specify the position of the origin" );
175 8 : readInputLine( dist_words );
176 :
177 8 : std::string inputLine = convertInputLineToString();
178 : // Make the kde object for the numerator if needed
179 8 : if( hasheights ) {
180 10 : readInputLine( getShortcutLabel() + "_inumer: KDE VOLUMES=" + data_str + " " + direction_string + " " + inputLine );
181 5 : if( unorm ) {
182 4 : readInputLine( getShortcutLabel() + ": ACCUMULATE ARG=" + getShortcutLabel() + "_inumer STRIDE=" + stride + " CLEAR=" + clear );
183 : return;
184 : } else {
185 6 : readInputLine( getShortcutLabel() + "_numer: ACCUMULATE ARG=" + getShortcutLabel() + "_inumer STRIDE=" + stride + " CLEAR=" + clear );
186 : }
187 : }
188 : // Make the kde object
189 12 : readInputLine( getShortcutLabel() + "_kde: KDE " + inputLine + " " + direction_string );
190 : // Make the division object if it is required
191 6 : if( hasheights && !unorm ) {
192 6 : readInputLine( getShortcutLabel() + "_denom: ACCUMULATE ARG=" + getShortcutLabel() + "_kde STRIDE=" + stride + " CLEAR=" + clear );
193 6 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_numer," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
194 3 : } else if( !hasheights ) {
195 3 : readInputLine( getShortcutLabel() + "_weight: ONES SIZE=1" );
196 6 : readInputLine( getShortcutLabel() + "_numer: ACCUMULATE ARG=" + getShortcutLabel() + "_kde STRIDE=" + stride + " CLEAR=" + clear );
197 6 : readInputLine( getShortcutLabel() + "_denom: ACCUMULATE ARG=" + getShortcutLabel() + "_weight STRIDE=" + stride + " CLEAR=" + clear );
198 6 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_numer," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
199 : }
200 0 : }
201 :
202 : }
203 : }
|
