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 "tools/Pbc.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/Atoms.h"
26 : #include "tools/Units.h"
27 : #include <cstdio>
28 : #include "core/ActionSet.h"
29 : #include "MultiColvarBase.h"
30 : #include "gridtools/ActionWithGrid.h"
31 :
32 : namespace PLMD {
33 : namespace multicolvar {
34 :
35 : //+PLUMEDOC GRIDCALC MULTICOLVARDENS
36 : /*
37 : Evaluate the average value of a multicolvar on a grid.
38 :
39 : This keyword allows one to construct a phase field representation for a symmetry function from
40 : an atomistic description. If each atom has an associated order parameter, \f$\phi_i\f$ then a
41 : smooth phase field function \f$\phi(r)\f$ can be computed using:
42 :
43 : \f[
44 : \phi(\mathbf{r}) = \frac{\sum_i K(\mathbf{r}-\mathbf{r}_i) \phi_i }{ \sum_i K(\mathbf{r} - \mathbf{r}_i )}
45 : \f]
46 :
47 : where \f$\mathbf{r}_i\f$ is the position of atom \f$i\f$, the sums run over all the atoms input
48 : and \f$K(\mathbf{r} - \mathbf{r}_i)\f$ is one of the \ref kernelfunctions implemented in plumed.
49 : This action calculates the above function on a grid, which can then be used in the input to further
50 : actions.
51 :
52 : \par Examples
53 :
54 : The following example shows perhaps the simplest way in which this action can be used. The following
55 : input computes the density of atoms at each point on the grid and outputs this quantity to a file. In
56 : other words this input instructs plumed to calculate \f$\rho(\mathbf{r}) = \sum_i K(\mathbf{r} - \mathbf{r}_i )\f$
57 :
58 : \plumedfile
59 : dens: DENSITY SPECIES=1-100
60 : grid: MULTICOLVARDENS DATA=dens ORIGIN=1 DIR=xyz NBINS=100,100,100 BANDWIDTH=0.05,0.05,0.05 STRIDE=1
61 : DUMPGRID GRID=grid STRIDE=500 FILE=density
62 : \endplumedfile
63 :
64 : In the above example density is added to the grid on every step. The PRINT_GRID instruction thus tells PLUMED to
65 : output the average density at each point on the grid every 500 steps of simulation. Notice that the that grid output
66 : on step 1000 is an average over all 1000 frames of the trajectory. If you would like to analyze these two blocks
67 : of data separately you must use the CLEAR flag.
68 :
69 : This second example computes an order parameter (in this case \ref FCCUBIC) and constructs a phase field model
70 : for this order parameter using the equation above.
71 :
72 : \plumedfile
73 : fcc: FCCUBIC SPECIES=1-5184 SWITCH={CUBIC D_0=1.2 D_MAX=1.5} ALPHA=27
74 : dens: MULTICOLVARDENS DATA=fcc ORIGIN=1 DIR=xyz NBINS=14,14,28 BANDWIDTH=1.0,1.0,1.0 STRIDE=1 CLEAR=1
75 : DUMPCUBE GRID=dens STRIDE=1 FILE=dens.cube
76 : \endplumedfile
77 :
78 : In this example the phase field model is computed and output to a file on every step of the simulation. Furthermore,
79 : because the CLEAR=1 keyword is set on the MULTICOLVARDENS line each Gaussian cube file output is a phase field
80 : model for a particular trajectory frame. The average value accumulated thus far is cleared at the start of every single
81 : timestep and there is no averaging over trajectory frames in this case.
82 :
83 : */
84 : //+ENDPLUMEDOC
85 :
86 : class MultiColvarDensity : public gridtools::ActionWithGrid {
87 : bool fractional;
88 : MultiColvarBase* mycolv;
89 : std::vector<unsigned> nbins;
90 : std::vector<double> gspacing;
91 : std::vector<bool> confined;
92 : std::vector<double> cmin, cmax;
93 : vesselbase::StoreDataVessel* stash;
94 : Vector origin;
95 : std::vector<unsigned> directions;
96 : public:
97 : explicit MultiColvarDensity(const ActionOptions&);
98 : static void registerKeywords( Keywords& keys );
99 : unsigned getNumberOfQuantities() const override;
100 0 : bool isPeriodic() override {
101 0 : return false;
102 : }
103 : void clearAverage() override;
104 : void prepareForAveraging() override;
105 : void compute( const unsigned&, MultiValue& ) const override;
106 34 : void apply() override {}
107 : };
108 :
109 13807 : PLUMED_REGISTER_ACTION(MultiColvarDensity,"MULTICOLVARDENS")
110 :
111 15 : void MultiColvarDensity::registerKeywords( Keywords& keys ) {
112 15 : gridtools::ActionWithGrid::registerKeywords( keys );
113 30 : keys.add("atoms","ORIGIN","we will use the position of this atom as the origin");
114 30 : keys.add("compulsory","DATA","the multicolvar which you would like to calculate the density profile for");
115 30 : keys.add("compulsory","DIR","the direction in which to calculate the density profile");
116 30 : keys.add("optional","NBINS","the number of bins to use to represent the density profile");
117 30 : keys.add("optional","SPACING","the approximate grid spacing (to be used as an alternative or together with NBINS)");
118 30 : keys.addFlag("FRACTIONAL",false,"use fractional coordinates for the various axes");
119 30 : keys.addFlag("XREDUCED",false,"limit the calculation of the density/average to a portion of the z-axis only");
120 30 : keys.add("optional","XLOWER","this is required if you are using XREDUCED. It specifies the lower bound for the region of the x-axis that for "
121 : "which you are calculating the density/average");
122 30 : keys.add("optional","XUPPER","this is required if you are using XREDUCED. It specifies the upper bound for the region of the x-axis that for "
123 : "which you are calculating the density/average");
124 30 : keys.addFlag("YREDUCED",false,"limit the calculation of the density/average to a portion of the y-axis only");
125 30 : keys.add("optional","YLOWER","this is required if you are using YREDUCED. It specifies the lower bound for the region of the y-axis that for "
126 : "which you are calculating the density/average");
127 30 : keys.add("optional","YUPPER","this is required if you are using YREDUCED. It specifies the upper bound for the region of the y-axis that for "
128 : "which you are calculating the density/average");
129 30 : keys.addFlag("ZREDUCED",false,"limit the calculation of the density/average to a portion of the z-axis only");
130 30 : keys.add("optional","ZLOWER","this is required if you are using ZREDUCED. It specifies the lower bound for the region of the z-axis that for "
131 : "which you are calculating the density/average");
132 30 : keys.add("optional","ZUPPER","this is required if you are using ZREDUCED. It specifies the upper bound for the region of the z-axis that for "
133 : "which you are calculating the density/average");
134 15 : }
135 :
136 11 : MultiColvarDensity::MultiColvarDensity(const ActionOptions&ao):
137 : Action(ao),
138 11 : ActionWithGrid(ao) {
139 : std::vector<AtomNumber> atom;
140 22 : parseAtomList("ORIGIN",atom);
141 11 : if( atom.size()!=1 ) {
142 0 : error("should only be one atom specified");
143 : }
144 11 : log.printf(" origin is at position of atom : %d\n",atom[0].serial() );
145 :
146 : std::string mlab;
147 11 : parse("DATA",mlab);
148 11 : mycolv = plumed.getActionSet().selectWithLabel<MultiColvarBase*>(mlab);
149 11 : if(!mycolv) {
150 0 : error("action labelled " + mlab + " does not exist or is not a MultiColvar");
151 : }
152 11 : stash = mycolv->buildDataStashes( NULL );
153 :
154 22 : parseFlag("FRACTIONAL",fractional);
155 : std::string direction;
156 11 : parse("DIR",direction);
157 11 : log.printf(" calculating for %s density profile along ", mycolv->getLabel().c_str() );
158 11 : if( direction=="x" ) {
159 8 : log.printf("x axis");
160 8 : directions.resize(1);
161 8 : directions[0]=0;
162 3 : } else if( direction=="y" ) {
163 0 : log.printf("y axis");
164 0 : directions.resize(1);
165 0 : directions[0]=1;
166 3 : } else if( direction=="z" ) {
167 0 : log.printf("z axis");
168 0 : directions.resize(1);
169 0 : directions[0]=2;
170 3 : } else if( direction=="xy" ) {
171 0 : log.printf("x and y axes");
172 0 : directions.resize(2);
173 0 : directions[0]=0;
174 0 : directions[1]=1;
175 3 : } else if( direction=="xz" ) {
176 0 : log.printf("x and z axes");
177 0 : directions.resize(2);
178 0 : directions[0]=0;
179 0 : directions[1]=2;
180 3 : } else if( direction=="yz" ) {
181 0 : log.printf("y and z axis");
182 0 : directions.resize(2);
183 0 : directions[0]=1;
184 0 : directions[1]=2;
185 3 : } else if( direction=="xyz" ) {
186 3 : log.printf("x, y and z axes");
187 3 : directions.resize(3);
188 3 : directions[0]=0;
189 3 : directions[1]=1;
190 3 : directions[2]=2;
191 : } else {
192 0 : error( direction + " is not valid gradient direction");
193 : }
194 11 : log.printf(" for colvars calculated by action %s \n",mycolv->getLabel().c_str() );
195 11 : parseVector("NBINS",nbins);
196 22 : parseVector("SPACING",gspacing);
197 11 : if( nbins.size()!=directions.size() && gspacing.size()!=directions.size() ) {
198 0 : error("NBINS or SPACING must be set");
199 : }
200 :
201 11 : confined.resize( directions.size() );
202 11 : cmin.resize( directions.size(), 0 );
203 11 : cmax.resize( directions.size(), 0 );
204 28 : for(unsigned i=0; i<directions.size(); ++i) {
205 17 : if( directions[i]==0 ) {
206 : bool tflag;
207 11 : parseFlag("XREDUCED",tflag);
208 11 : confined[i]=tflag;
209 11 : if( confined[i] ) {
210 0 : cmin[i]=cmax[i]=0.0;
211 0 : parse("XLOWER",cmin[i]);
212 0 : parse("XUPPER",cmax[i]);
213 0 : if( fractional ) {
214 0 : error("XREDUCED is incompatible with FRACTIONAL");
215 : }
216 0 : if( std::abs(cmin[i]-cmax[i])<epsilon ) {
217 0 : error("range set for x axis makes no sense");
218 : }
219 0 : log.printf(" confining calculation in x direction to between %f and %f \n",cmin[i],cmax[i]);
220 : }
221 6 : } else if( directions[i]==1 ) {
222 : bool tflag;
223 3 : parseFlag("YREDUCED",tflag);
224 3 : confined[i]=tflag;
225 3 : if( confined[i] ) {
226 0 : cmin[i]=cmax[i]=0.0;
227 0 : parse("YLOWER",cmin[i]);
228 0 : parse("YUPPER",cmax[i]);
229 0 : if( fractional ) {
230 0 : error("YREDUCED is incompatible with FRACTIONAL");
231 : }
232 0 : if( std::abs(cmin[i]-cmax[i])<epsilon ) {
233 0 : error("range set for y axis makes no sense");
234 : }
235 0 : log.printf(" confining calculation in y direction to between %f and %f \n",cmin[i],cmax[i]);
236 : }
237 3 : } else if( directions[i]==2 ) {
238 : bool tflag;
239 3 : parseFlag("ZREDUCED",tflag);
240 3 : confined[i]=tflag;
241 3 : if( confined[i] ) {
242 1 : cmin[i]=cmax[i]=0.0;
243 2 : parse("ZLOWER",cmin[i]);
244 1 : parse("ZUPPER",cmax[i]);
245 1 : if( fractional ) {
246 0 : error("ZREDUCED is incompatible with FRACTIONAL");
247 : }
248 1 : if( std::abs(cmin[i]-cmax[i])<epsilon ) {
249 0 : error("range set for z axis search makes no sense");
250 : }
251 1 : log.printf(" confining calculation in z direction to between %f and %f \n",cmin[i],cmax[i]);
252 : }
253 : }
254 : }
255 :
256 : std::string vstring;
257 11 : if( confined[0] ) {
258 : vstring +="PBC=F";
259 : } else {
260 : vstring += " PBC=T";
261 : }
262 17 : for(unsigned i=1; i<directions.size(); ++i) {
263 6 : if( confined[i] ) {
264 : vstring += ",F";
265 : } else {
266 : vstring += ",T";
267 : }
268 : }
269 22 : vstring +=" COMPONENTS=" + mycolv->getLabel() + ".dens";
270 : vstring +=" COORDINATES=";
271 11 : if( directions[0]==0 ) {
272 : vstring+="x";
273 0 : } else if( directions[0]==1 ) {
274 : vstring+="y";
275 0 : } else if( directions[0]==2 ) {
276 : vstring+="z";
277 : }
278 17 : for(unsigned i=1; i<directions.size(); ++i) {
279 6 : if( directions[i]==0 ) {
280 : vstring+=",x";
281 6 : } else if( directions[i]==1 ) {
282 : vstring+=",y";
283 3 : } else if( directions[i]==2 ) {
284 : vstring+=",z";
285 : }
286 : }
287 : // Create a task list
288 10899 : for(unsigned i=0; i<mycolv->getFullNumberOfTasks(); ++i) {
289 10888 : addTaskToList(i);
290 : }
291 : // And create the grid
292 11 : std::unique_ptr<gridtools::GridVessel> grid;
293 11 : if( mycolv->isDensity() ) {
294 6 : grid=createGrid( "histogram", vstring );
295 : } else {
296 16 : grid=createGrid( "average", vstring );
297 : }
298 : // cppcheck-suppress danglingLifetime
299 11 : mygrid=grid.get();
300 : // And finish the grid setup
301 11 : setAveragingAction( std::move(grid), true );
302 :
303 : // Enusre units for cube files are set correctly
304 11 : if( !fractional ) {
305 11 : if( plumed.getAtoms().usingNaturalUnits() ) {
306 10 : mygrid->setCubeUnits( 1.0/0.5292 );
307 : } else {
308 1 : mygrid->setCubeUnits( plumed.getAtoms().getUnits().getLength()/.05929 );
309 : }
310 : }
311 :
312 11 : checkRead();
313 11 : requestAtoms(atom);
314 : // Stupid dependencies cleared by requestAtoms - why GBussi why? That's got me so many times
315 11 : addDependency( mycolv );
316 22 : }
317 :
318 62 : unsigned MultiColvarDensity::getNumberOfQuantities() const {
319 62 : return directions.size() + 2;
320 : }
321 :
322 18 : void MultiColvarDensity::clearAverage() {
323 18 : std::vector<double> min(directions.size()), max(directions.size());
324 18 : std::vector<std::string> gmin(directions.size()), gmax(directions.size());;
325 46 : for(unsigned i=0; i<directions.size(); ++i) {
326 28 : min[i]=-0.5;
327 28 : max[i]=0.5;
328 : }
329 18 : if( !fractional ) {
330 18 : if( !mycolv->getPbc().isOrthorombic() ) {
331 0 : error("I think that density profiles with non-orthorhombic cells don't work. If you want it have a look and see if you can work it out");
332 : }
333 :
334 46 : for(unsigned i=0; i<directions.size(); ++i) {
335 28 : if( !confined[i] ) {
336 25 : min[i]*=mycolv->getBox()(directions[i],directions[i]);
337 25 : max[i]*=mycolv->getBox()(directions[i],directions[i]);
338 : } else {
339 3 : min[i]=cmin[i];
340 3 : max[i]=cmax[i];
341 : }
342 : }
343 : }
344 46 : for(unsigned i=0; i<directions.size(); ++i) {
345 28 : Tools::convert(min[i],gmin[i]);
346 28 : Tools::convert(max[i],gmax[i]);
347 : }
348 18 : ActionWithAveraging::clearAverage();
349 18 : mygrid->setBounds( gmin, gmax, nbins, gspacing );
350 18 : resizeFunctions();
351 36 : }
352 :
353 26 : void MultiColvarDensity::prepareForAveraging() {
354 62 : for(unsigned i=0; i<directions.size(); ++i) {
355 36 : if( confined[i] ) {
356 3 : continue;
357 : }
358 : std::string max;
359 33 : Tools::convert( 0.5*mycolv->getBox()(directions[i],directions[i]), max );
360 33 : if( max!=mygrid->getMax()[i] ) {
361 0 : error("box size should be fixed. Use FRACTIONAL");
362 : }
363 : }
364 : // Ensure we only work with active multicolvars
365 26 : deactivateAllTasks();
366 21178 : for(unsigned i=0; i<stash->getNumberOfStoredValues(); ++i) {
367 21152 : taskFlags[i]=1;
368 : }
369 26 : lockContributors();
370 : // Retrieve the origin
371 26 : origin = getPosition(0);
372 26 : }
373 :
374 21152 : void MultiColvarDensity::compute( const unsigned& current, MultiValue& myvals ) const {
375 21152 : std::vector<double> cvals( mycolv->getNumberOfQuantities() );
376 21152 : stash->retrieveSequentialValue( current, false, cvals );
377 21152 : Vector fpos, apos = pbcDistance( origin, mycolv->getCentralAtomPos( mycolv->getPositionInFullTaskList(current) ) );
378 21152 : if( fractional ) {
379 0 : fpos = getPbc().realToScaled( apos );
380 : } else {
381 21152 : fpos=apos;
382 : }
383 :
384 21152 : myvals.setValue( 0, cweight*cvals[0] );
385 84566 : for(unsigned j=0; j<directions.size(); ++j) {
386 63414 : myvals.setValue( 1+j, fpos[ directions[j] ] );
387 : }
388 21152 : myvals.setValue( 1+directions.size(), cvals[1] );
389 21152 : }
390 :
391 : }
392 : }
|
