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 { return false; }
101 : void clearAverage() override;
102 : void prepareForAveraging() override;
103 : void compute( const unsigned&, MultiValue& ) const override;
104 34 : void apply() override {}
105 : };
106 :
107 10441 : PLUMED_REGISTER_ACTION(MultiColvarDensity,"MULTICOLVARDENS")
108 :
109 12 : void MultiColvarDensity::registerKeywords( Keywords& keys ) {
110 12 : gridtools::ActionWithGrid::registerKeywords( keys );
111 24 : keys.add("atoms","ORIGIN","we will use the position of this atom as the origin");
112 24 : keys.add("compulsory","DATA","the multicolvar which you would like to calculate the density profile for");
113 24 : keys.add("compulsory","DIR","the direction in which to calculate the density profile");
114 24 : keys.add("optional","NBINS","the number of bins to use to represent the density profile");
115 24 : keys.add("optional","SPACING","the approximate grid spacing (to be used as an alternative or together with NBINS)");
116 24 : keys.addFlag("FRACTIONAL",false,"use fractional coordinates for the various axes");
117 24 : keys.addFlag("XREDUCED",false,"limit the calculation of the density/average to a portion of the z-axis only");
118 24 : 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 "
119 : "which you are calculating the density/average");
120 24 : 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 "
121 : "which you are calculating the density/average");
122 24 : keys.addFlag("YREDUCED",false,"limit the calculation of the density/average to a portion of the y-axis only");
123 24 : 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 "
124 : "which you are calculating the density/average");
125 24 : 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 "
126 : "which you are calculating the density/average");
127 24 : keys.addFlag("ZREDUCED",false,"limit the calculation of the density/average to a portion of the z-axis only");
128 24 : 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 "
129 : "which you are calculating the density/average");
130 24 : 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 "
131 : "which you are calculating the density/average");
132 12 : }
133 :
134 11 : MultiColvarDensity::MultiColvarDensity(const ActionOptions&ao):
135 : Action(ao),
136 11 : ActionWithGrid(ao)
137 : {
138 : std::vector<AtomNumber> atom;
139 22 : parseAtomList("ORIGIN",atom);
140 11 : if( atom.size()!=1 ) error("should only be one atom specified");
141 11 : log.printf(" origin is at position of atom : %d\n",atom[0].serial() );
142 :
143 11 : std::string mlab; parse("DATA",mlab);
144 11 : mycolv = plumed.getActionSet().selectWithLabel<MultiColvarBase*>(mlab);
145 11 : if(!mycolv) error("action labelled " + mlab + " does not exist or is not a MultiColvar");
146 11 : stash = mycolv->buildDataStashes( NULL );
147 :
148 22 : parseFlag("FRACTIONAL",fractional);
149 11 : std::string direction; parse("DIR",direction);
150 11 : log.printf(" calculating for %s density profile along ", mycolv->getLabel().c_str() );
151 11 : if( direction=="x" ) {
152 8 : log.printf("x axis");
153 8 : directions.resize(1); directions[0]=0;
154 3 : } else if( direction=="y" ) {
155 0 : log.printf("y axis");
156 0 : directions.resize(1); directions[0]=1;
157 3 : } else if( direction=="z" ) {
158 0 : log.printf("z axis");
159 0 : directions.resize(1); directions[0]=2;
160 3 : } else if( direction=="xy" ) {
161 0 : log.printf("x and y axes");
162 0 : directions.resize(2); directions[0]=0; directions[1]=1;
163 3 : } else if( direction=="xz" ) {
164 0 : log.printf("x and z axes");
165 0 : directions.resize(2); directions[0]=0; directions[1]=2;
166 3 : } else if( direction=="yz" ) {
167 0 : log.printf("y and z axis");
168 0 : directions.resize(2); directions[0]=1; directions[1]=2;
169 3 : } else if( direction=="xyz" ) {
170 3 : log.printf("x, y and z axes");
171 3 : directions.resize(3); directions[0]=0; directions[1]=1; directions[2]=2;
172 : } else {
173 0 : error( direction + " is not valid gradient direction");
174 : }
175 11 : log.printf(" for colvars calculated by action %s \n",mycolv->getLabel().c_str() );
176 33 : parseVector("NBINS",nbins); parseVector("SPACING",gspacing);
177 11 : if( nbins.size()!=directions.size() && gspacing.size()!=directions.size() ) error("NBINS or SPACING must be set");
178 :
179 11 : confined.resize( directions.size() ); cmin.resize( directions.size(), 0 ); cmax.resize( directions.size(), 0 );
180 28 : for(unsigned i=0; i<directions.size(); ++i) {
181 17 : if( directions[i]==0 ) {
182 22 : bool tflag; parseFlag("XREDUCED",tflag); confined[i]=tflag;
183 11 : if( confined[i] ) {
184 0 : cmin[i]=cmax[i]=0.0; parse("XLOWER",cmin[i]); parse("XUPPER",cmax[i]);
185 0 : if( fractional ) error("XREDUCED is incompatible with FRACTIONAL");
186 0 : if( std::abs(cmin[i]-cmax[i])<epsilon ) error("range set for x axis makes no sense");
187 0 : log.printf(" confining calculation in x direction to between %f and %f \n",cmin[i],cmax[i]);
188 : }
189 6 : } else if( directions[i]==1 ) {
190 6 : bool tflag; parseFlag("YREDUCED",tflag); confined[i]=tflag;
191 3 : if( confined[i] ) {
192 0 : cmin[i]=cmax[i]=0.0; parse("YLOWER",cmin[i]); parse("YUPPER",cmax[i]);
193 0 : if( fractional ) error("YREDUCED is incompatible with FRACTIONAL");
194 0 : if( std::abs(cmin[i]-cmax[i])<epsilon ) error("range set for y axis makes no sense");
195 0 : log.printf(" confining calculation in y direction to between %f and %f \n",cmin[i],cmax[i]);
196 : }
197 3 : } else if( directions[i]==2 ) {
198 6 : bool tflag; parseFlag("ZREDUCED",tflag); confined[i]=tflag;
199 3 : if( confined[i] ) {
200 2 : cmin[i]=cmax[i]=0.0; parse("ZLOWER",cmin[i]); parse("ZUPPER",cmax[i]);
201 1 : if( fractional ) error("ZREDUCED is incompatible with FRACTIONAL");
202 1 : if( std::abs(cmin[i]-cmax[i])<epsilon ) error("range set for z axis search makes no sense");
203 1 : log.printf(" confining calculation in z direction to between %f and %f \n",cmin[i],cmax[i]);
204 : }
205 : }
206 : }
207 :
208 : std::string vstring;
209 11 : if( confined[0] ) vstring +="PBC=F";
210 : else vstring += " PBC=T";
211 17 : for(unsigned i=1; i<directions.size(); ++i) {
212 6 : if( confined[i] ) vstring += ",F";
213 : else vstring += ",T";
214 : }
215 22 : vstring +=" COMPONENTS=" + mycolv->getLabel() + ".dens";
216 : vstring +=" COORDINATES=";
217 11 : if( directions[0]==0 ) vstring+="x";
218 0 : else if( directions[0]==1 ) vstring+="y";
219 0 : else if( directions[0]==2 ) vstring+="z";
220 17 : for(unsigned i=1; i<directions.size(); ++i) {
221 6 : if( directions[i]==0 ) vstring+=",x";
222 6 : else if( directions[i]==1 ) vstring+=",y";
223 3 : else if( directions[i]==2 ) vstring+=",z";
224 : }
225 : // Create a task list
226 10899 : for(unsigned i=0; i<mycolv->getFullNumberOfTasks(); ++i) addTaskToList(i);
227 : // And create the grid
228 11 : std::unique_ptr<gridtools::GridVessel> grid;
229 14 : if( mycolv->isDensity() ) grid=createGrid( "histogram", vstring );
230 16 : else grid=createGrid( "average", vstring );
231 11 : mygrid=grid.get();
232 : // And finish the grid setup
233 11 : setAveragingAction( std::move(grid), true );
234 :
235 : // Enusre units for cube files are set correctly
236 11 : if( !fractional ) {
237 11 : if( plumed.getAtoms().usingNaturalUnits() ) mygrid->setCubeUnits( 1.0/0.5292 );
238 1 : else mygrid->setCubeUnits( plumed.getAtoms().getUnits().getLength()/.05929 );
239 : }
240 :
241 11 : checkRead(); requestAtoms(atom);
242 : // Stupid dependencies cleared by requestAtoms - why GBussi why? That's got me so many times
243 11 : addDependency( mycolv );
244 22 : }
245 :
246 62 : unsigned MultiColvarDensity::getNumberOfQuantities() const {
247 62 : return directions.size() + 2;
248 : }
249 :
250 18 : void MultiColvarDensity::clearAverage() {
251 18 : std::vector<double> min(directions.size()), max(directions.size());
252 18 : std::vector<std::string> gmin(directions.size()), gmax(directions.size());;
253 46 : for(unsigned i=0; i<directions.size(); ++i) { min[i]=-0.5; max[i]=0.5; }
254 18 : if( !fractional ) {
255 18 : if( !mycolv->getPbc().isOrthorombic() ) {
256 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");
257 : }
258 :
259 46 : for(unsigned i=0; i<directions.size(); ++i) {
260 28 : if( !confined[i] ) {
261 25 : min[i]*=mycolv->getBox()(directions[i],directions[i]);
262 25 : max[i]*=mycolv->getBox()(directions[i],directions[i]);
263 : } else {
264 3 : min[i]=cmin[i]; max[i]=cmax[i];
265 : }
266 : }
267 : }
268 46 : for(unsigned i=0; i<directions.size(); ++i) { Tools::convert(min[i],gmin[i]); Tools::convert(max[i],gmax[i]); }
269 18 : ActionWithAveraging::clearAverage();
270 18 : mygrid->setBounds( gmin, gmax, nbins, gspacing ); resizeFunctions();
271 36 : }
272 :
273 26 : void MultiColvarDensity::prepareForAveraging() {
274 62 : for(unsigned i=0; i<directions.size(); ++i) {
275 36 : if( confined[i] ) continue;
276 33 : std::string max; Tools::convert( 0.5*mycolv->getBox()(directions[i],directions[i]), max );
277 33 : if( max!=mygrid->getMax()[i] ) error("box size should be fixed. Use FRACTIONAL");
278 : }
279 : // Ensure we only work with active multicolvars
280 26 : deactivateAllTasks();
281 21178 : for(unsigned i=0; i<stash->getNumberOfStoredValues(); ++i) taskFlags[i]=1;
282 26 : lockContributors();
283 : // Retrieve the origin
284 26 : origin = getPosition(0);
285 26 : }
286 :
287 21152 : void MultiColvarDensity::compute( const unsigned& current, MultiValue& myvals ) const {
288 21152 : std::vector<double> cvals( mycolv->getNumberOfQuantities() ); stash->retrieveSequentialValue( current, false, cvals );
289 21152 : Vector fpos, apos = pbcDistance( origin, mycolv->getCentralAtomPos( mycolv->getPositionInFullTaskList(current) ) );
290 21152 : if( fractional ) { fpos = getPbc().realToScaled( apos ); } else { fpos=apos; }
291 :
292 84566 : myvals.setValue( 0, cweight*cvals[0] ); for(unsigned j=0; j<directions.size(); ++j) myvals.setValue( 1+j, fpos[ directions[j] ] );
293 21152 : myvals.setValue( 1+directions.size(), cvals[1] );
294 21152 : }
295 :
296 : }
297 : }
|
