Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2020 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 "ActionVolume.h"
25 : #include "VolumeShortcut.h"
26 :
27 : //+PLUMEDOC VOLUMES AROUND
28 : /*
29 : This quantity can be used to calculate functions of the distribution of collective variables for the atoms that lie in a particular, user-specified part of of the cell.
30 :
31 : This action can be used to calculate whether each of the atoms are within a particular part of the simulation box or not as illustrated by the following example:
32 :
33 : ```plumed
34 : f: FIXEDATOM AT=0,0,0
35 : a: AROUND ATOMS=1-100 ORIGIN=f SIGMA=0.2 XLOWER=-1.0 XUPPER=1.0 YLOWER=-1.0 YUPPER=1.0 ZLOWER=-1.0 ZUPPER=1.0
36 : PRINT ARG=a FILE=colvar
37 : ```
38 :
39 : The 100 elements of the vector `a` that is returned from the AROUND action in the above input are calculated using:
40 :
41 : $$
42 : w(x_i,y_i,z_i) = \int_{xl}^{xu} \int_{yl}^{yu} \int_{zl}^{zu} \textrm{d}x\textrm{d}y\textrm{d}z K\left( \frac{x - x_i}{\sigma} \right)K\left( \frac{y - y_i}{\sigma} \right)K\left( \frac{z - z_i}{\sigma} \right)
43 : $$
44 :
45 : where $K$ is one of the kernel functions described in the documentation for the function [BETWEEN](BETWEEN.md), $\sigma$ is a bandwidth parameter and the limits
46 : for the integrals are the values specified using the keywords XLOWER, XUPPER, YLOWER, YUPPER, YUPPER, ZLOWER and ZUPPER. $x_i$, $y_i$ and $z_i$ are then the components
47 : of the vector that connects the $i$th atom that was specified using the ATOMS keyword to the atom that was specified using the ORIGIN keyword. In other words,
48 : $w(x_i,y_i,z_i)$ is 1 if the atom is within a rectangular box that is centered on the atom that is specified as the origin and zero otherwise.
49 :
50 : ## Calculating the number of atoms in a particular part of the box
51 :
52 : Lets suppose that you want to calculate how many atoms are in have an $x$ coordinate that is between -1.0 and 1.0. You can do this using the following PLUMED input:
53 :
54 : ```plumed
55 : f: FIXEDATOM AT=0,0,0
56 : a: AROUND ATOMS=1-100 ORIGIN=f SIGMA=0.2 XLOWER=-1.0 XUPPER=1.0
57 : s: SUM ARG=a PERIODIC=NO
58 : PRINT ARG=s FILE=colvar
59 : ```
60 :
61 : In this example the components of `a` are calculated as:
62 :
63 : $$
64 : w(x_i,y_i,z_i) = \int_{xl}^{xu} \textrm{d}x K\left( \frac{x - x_i}{\sigma} \right)
65 : $$
66 :
67 : as the YLOWER, YUPPER, YUPPER, ZLOWER and ZUPPER flags have not been included. The [SUM](SUM.md) command then adds together all the elements of the vector `a` to calculate the total number of atoms in the region
68 : of the box that is of interest.
69 :
70 : ## Calculating the average value for an order parameter in a particular part of the box
71 :
72 : Suppose that you have calculated a vector of order parameters that can be assigned to a particular point in three dimensional space.
73 : The symmetry functions in the [symfunc](module_symfunc.md) module are examples of order parameters that satisfy this criteria. You can use
74 : the AROUND command to calculate the average value of the symmetry function in a particular part of the box as follows:
75 :
76 : ```plumed
77 : c: COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=1.0}
78 : f: FIXEDATOM AT=0,0,0
79 : a: AROUND ATOMS=1-100 ORIGIN=f SIGMA=0.2 XLOWER=-1.0 XUPPER=1.0 YLOWER=-1.0 YUPPER=1.0
80 : p: CUSTOM ARG=c,a FUNC=x*y PERIODIC=NO
81 : n: SUM ARG=p PERIODIC=NO
82 : d: SUM ARG=a PERIODIC=NO
83 : av: CUSTOM ARG=n,d FUNC=x/y PERIODIC=NO
84 : PRINT ARG=av FILE=colvar
85 : ```
86 :
87 : The final quantity `av` here is:
88 :
89 : $$
90 : \overline{s}_{\tau} = \frac{ \sum_i c_i w(x_i,y_i,z_i) }{ \sum_i w(x_i,y_i,z_i) }
91 : $$
92 :
93 : where $c_i$ are the coordination numbers and $w_i$ is:
94 :
95 : $$
96 : w(x_i,y_i,z_i) = \int_{xl}^{xu} \int_{yl}^{yu} \textrm{d}x \textrm{d}y K\left( \frac{x - x_i}{\sigma} \right) K\left( \frac{y - y_i}{\sigma} \right)
97 : $$
98 :
99 : ## Old syntax
100 :
101 : In earlier versions of PLUMED the syntax for the calculation in the previous section is as follows:
102 :
103 : ```plumed
104 : c: COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=1.0}
105 : f: FIXEDATOM AT=0,0,0
106 : a: AROUND DATA=c ORIGIN=f SIGMA=0.2 XLOWER=-1.0 XUPPER=1.0 YLOWER=-1.0 YUPPER=1.0 MEAN
107 : PRINT ARG=a.mean FILE=colvar
108 : ```
109 :
110 : This old syntax still works but we highly recommend you use the newer syntax as it is easlier to understand
111 : and more flexible. You will also notice that AROUND in the input above is a shortcut that expands to the longer input
112 : that was given in the previous section.
113 :
114 : */
115 : //+ENDPLUMEDOC
116 :
117 : //+PLUMEDOC MCOLVAR AROUND_CALC
118 : /*
119 : Calculate a vector from the input positions with elements equal to one when the positions are in a particular part of the cell and elements equal to zero otherwise
120 :
121 : \par Examples
122 :
123 : */
124 : //+ENDPLUMEDOC
125 :
126 : namespace PLMD {
127 : namespace volumes {
128 :
129 : class VolumeAround : public ActionVolume {
130 : private:
131 : Vector origin;
132 : bool dox, doy, doz;
133 : double xlow, xhigh;
134 : double ylow, yhigh;
135 : double zlow, zhigh;
136 : public:
137 : static void registerKeywords( Keywords& keys );
138 : explicit VolumeAround(const ActionOptions& ao);
139 : void setupRegions() override;
140 : double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const override;
141 : };
142 :
143 : PLUMED_REGISTER_ACTION(VolumeAround,"AROUND_CALC")
144 : char glob_around[] = "AROUND";
145 : typedef VolumeShortcut<glob_around> VolumeAroundShortcut;
146 : PLUMED_REGISTER_ACTION(VolumeAroundShortcut,"AROUND")
147 :
148 156 : void VolumeAround::registerKeywords( Keywords& keys ) {
149 156 : ActionVolume::registerKeywords( keys );
150 312 : keys.setDisplayName("AROUND");
151 156 : keys.add("atoms","ORIGIN","the atom whose vicinity we are interested in examining");
152 156 : keys.add("atoms-2","ATOM","an alternative to ORIGIN");
153 156 : keys.add("compulsory","XLOWER","0.0","the lower boundary in x relative to the x coordinate of the atom (0 indicates use full extent of box).");
154 156 : keys.add("compulsory","XUPPER","0.0","the upper boundary in x relative to the x coordinate of the atom (0 indicates use full extent of box).");
155 156 : keys.add("compulsory","YLOWER","0.0","the lower boundary in y relative to the y coordinate of the atom (0 indicates use full extent of box).");
156 156 : keys.add("compulsory","YUPPER","0.0","the upper boundary in y relative to the y coordinate of the atom (0 indicates use full extent of box).");
157 156 : keys.add("compulsory","ZLOWER","0.0","the lower boundary in z relative to the z coordinate of the atom (0 indicates use full extent of box).");
158 156 : keys.add("compulsory","ZUPPER","0.0","the upper boundary in z relative to the z coordinate of the atom (0 indicates use full extent of box).");
159 156 : }
160 :
161 46 : VolumeAround::VolumeAround(const ActionOptions& ao):
162 : Action(ao),
163 46 : ActionVolume(ao) {
164 : std::vector<AtomNumber> atom;
165 92 : parseAtomList("ORIGIN",atom);
166 46 : if( atom.size()==0 ) {
167 0 : parseAtomList("ATOM",atom);
168 : }
169 46 : if( atom.size()!=1 ) {
170 0 : error("should only be one atom specified");
171 : }
172 46 : log.printf(" boundaries for region are calculated based on positions of atom : %d\n",atom[0].serial() );
173 :
174 46 : dox=true;
175 46 : parse("XLOWER",xlow);
176 46 : parse("XUPPER",xhigh);
177 46 : doy=true;
178 46 : parse("YLOWER",ylow);
179 46 : parse("YUPPER",yhigh);
180 46 : doz=true;
181 46 : parse("ZLOWER",zlow);
182 46 : parse("ZUPPER",zhigh);
183 46 : if( xlow==0.0 && xhigh==0.0 ) {
184 6 : dox=false;
185 : }
186 46 : if( ylow==0.0 && yhigh==0.0 ) {
187 14 : doy=false;
188 : }
189 46 : if( zlow==0.0 && zhigh==0.0 ) {
190 14 : doz=false;
191 : }
192 46 : if( !dox && !doy && !doz ) {
193 0 : error("no subregion defined use XLOWER, XUPPER, YLOWER, YUPPER, ZLOWER, ZUPPER");
194 : }
195 46 : log.printf(" boundaries for region (region of interest about atom) : x %f %f, y %f %f, z %f %f \n",xlow,xhigh,ylow,yhigh,zlow,zhigh);
196 46 : checkRead();
197 46 : requestAtoms(atom);
198 46 : }
199 :
200 388 : void VolumeAround::setupRegions() { }
201 :
202 98837 : double VolumeAround::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const {
203 : // Setup the histogram bead
204 98837 : HistogramBead bead;
205 : bead.isNotPeriodic();
206 98837 : bead.setKernelType( getKernelType() );
207 :
208 : // Calculate position of atom wrt to origin
209 98837 : Vector fpos=pbcDistance( getPosition(0), cpos );
210 : double xcontr, ycontr, zcontr, xder, yder, zder;
211 98837 : if( dox ) {
212 66341 : bead.set( xlow, xhigh, getSigma() );
213 66341 : xcontr=bead.calculate( fpos[0], xder );
214 : } else {
215 : xcontr=1.;
216 32496 : xder=0.;
217 : }
218 98837 : if( doy ) {
219 63681 : bead.set( ylow, yhigh, getSigma() );
220 63681 : ycontr=bead.calculate( fpos[1], yder );
221 : } else {
222 : ycontr=1.;
223 35156 : yder=0.;
224 : }
225 98837 : if( doz ) {
226 62589 : bead.set( zlow, zhigh, getSigma() );
227 62589 : zcontr=bead.calculate( fpos[2], zder );
228 : } else {
229 : zcontr=1.;
230 36248 : zder=0.;
231 : }
232 98837 : derivatives[0]=xder*ycontr*zcontr;
233 98837 : derivatives[1]=xcontr*yder*zcontr;
234 98837 : derivatives[2]=xcontr*ycontr*zder;
235 : // Add derivatives wrt to position of origin atom
236 98837 : refders[0] = -derivatives;
237 : // Add virial contribution
238 98837 : vir -= Tensor(fpos,derivatives);
239 98837 : return xcontr*ycontr*zcontr;
240 : }
241 :
242 : }
243 : }
|
