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