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