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