Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 : #ifdef __PLUMED_HAS_OPENACC
23 : #define __PLUMED_USE_OPENACC 1
24 : #endif //__PLUMED_HAS_OPENACC
25 : #include "core/ActionRegister.h"
26 : #include "tools/Pbc.h"
27 : #include "tools/SwitchingFunction.h"
28 : #include "ActionVolume.h"
29 : #include "tools/HistogramBead.h"
30 : #include "VolumeShortcut.h"
31 :
32 : //+PLUMEDOC VOLUMES INCYLINDER
33 : /*
34 : 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.
35 :
36 : 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:
37 :
38 : ```plumed
39 : f: FIXEDATOM AT=0,0,0
40 : a: INCYLINDER ...
41 : ATOMS=1-100 CENTER=f DIRECTION=Z
42 : RADIUS={TANH R_0=1.5} SIGMA=0.1
43 : LOWER=-1.0 UPPER=1.0
44 : KERNEL=gaussian
45 : ...
46 : PRINT ARG=a FILE=colvar
47 : ```
48 :
49 : The 100 elements of the vector `a` that is returned from the INCYLINDER action in the above input are calculated using:
50 :
51 : $$
52 : w(x_i,y_i,z_i) = s\left( r_{xy} \right) \int_{zl}^{zu} \textrm{d}z K\left( \frac{z - z_i}{\sigma} \right)
53 : $$
54 :
55 : where
56 :
57 : $$
58 : r_{xy} = \sqrt{ x_i^2 + y_i^2 }
59 : $$
60 :
61 : In these expressions $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
62 : for the integrals are the values specified using the keywords `LOWER` and `UPPER`. $x_i$, $y_i$ and $z_i$ are then the components
63 : of the vector that connects the $i$th atom that was specified using the `ATOMS` keyword to the atom that was specified using the `CENTER` keyword. In other words,
64 : $w(x_i,y_i,z_i)$ is 1 if the atom is within a cylinder that is centered on the $z$ axis that has an extent along the $z$ direction around the position of atom `f` that
65 : is determined by the values specified by `LOWER` and `UPPER` keywords. The radial extent of this cylinder is determined by the parameters of the switching function that is
66 : specified using the `RADIUS` keyword.
67 :
68 : If you want to caluclate and print the number of atoms in the cylinder of interest you can use an input like the one shown below:
69 :
70 : ```plumed
71 : f: FIXEDATOM AT=0,0,0
72 : a: INCYLINDER ...
73 : ATOMS=1-100 CENTER=f DIRECTION=X
74 : RADIUS={TANH R_0=1.5} SIGMA=0.1
75 : LOWER=-1.0 UPPER=1.0
76 : ...
77 : s: SUM ARG=a PERIODIC=NO
78 : PRINT ARG=s FILE=colvar
79 : ```
80 :
81 : Alternatively, if you want to calculate an print the number of atoms that are not in the cylinder of interest you can use the OUTSIDE flag as shown below:
82 :
83 : ```plumed
84 : f: FIXEDATOM AT=0,0,0
85 : a: INCYLINDER ...
86 : ATOMS=1-100 CENTER=f DIRECTION=X
87 : RADIUS={TANH R_0=1.5} SIGMA=0.1
88 : LOWER=-1.0 UPPER=1.0
89 : OUTSIDE
90 : ...
91 : s: SUM ARG=a PERIODIC=NO
92 : PRINT ARG=s FILE=colvar
93 : ```
94 :
95 : Notice that now the cylinder is centered on the `x` axis rather than the `z` axis as we have changed the input for the `DIRECTION` keyword.
96 :
97 : !!! note ""
98 :
99 : You can also calculate the average values of symmetry functions in the cylinder of interest by using inputs similar to those described the documentation for the [AROUND](AROUND.md)
100 : action. In other words, you can swap out AROUND actions for an INCLYLINDER actions. Also as with [AROUND](AROUND.md), earlier versions of PLUMED used a different syntax for doing these
101 : types of calculations, which can still be used with this new version of the command. We strongly recommend using the newer syntax but if you are interested in the
102 : old syntax you can find more information in the old syntax section of the documentation for [AROUND](AROUND.md).
103 :
104 :
105 : */
106 : //+ENDPLUMEDOC
107 :
108 : namespace PLMD {
109 : namespace volumes {
110 :
111 3 : class VolumeInCylinder {
112 : public:
113 : bool docylinder;
114 : double min, max, sigma;
115 : HistogramBead::KernelType kerneltype;
116 : std::array<unsigned,3> dir;
117 : #ifdef __PLUMED_USE_OPENACC
118 : SwitchingFunctionAccelerable switchingFunction;
119 : #else
120 : SwitchingFunction switchingFunction;
121 : #endif //__PLUMED_USE_OPENACC
122 : static void registerKeywords( Keywords& keys );
123 : void parseInput( ActionVolume<VolumeInCylinder>* action );
124 : void setupRegions( ActionVolume<VolumeInCylinder>* action, const Pbc& pbc, const std::vector<Vector>& positions ) {}
125 : static void parseAtoms( ActionVolume<VolumeInCylinder>* action, std::vector<AtomNumber>& atom );
126 : static void calculateNumberInside( const VolumeInput& input, const VolumeInCylinder& actioninput, VolumeOutput& output );
127 : #ifdef __PLUMED_USE_OPENACC
128 : void toACCDevice() const {
129 : #pragma acc enter data copyin(this[0:1], \
130 : docylinder, min, max, sigma, kerneltype, dir[0:3])
131 : switchingFunction.toACCDevice();
132 : }
133 : void removeFromACCDevice() const {
134 : switchingFunction.removeFromACCDevice();
135 : #pragma acc exit data delete(dir[0:3], kerneltype, sigma, max, min, \
136 : docylinder, this[0:1])
137 : }
138 : #endif //__PLUMED_USE_OPENACC
139 : };
140 :
141 : typedef ActionVolume<VolumeInCylinder> Volc;
142 : PLUMED_REGISTER_ACTION(Volc,"INCYLINDER_CALC")
143 : char glob_cylinder[] = "INCYLINDER";
144 : typedef VolumeShortcut<glob_cylinder> VolumeInCylinderShortcut;
145 : PLUMED_REGISTER_ACTION(VolumeInCylinderShortcut,"INCYLINDER")
146 :
147 7 : void VolumeInCylinder::registerKeywords( Keywords& keys ) {
148 14 : keys.setDisplayName("INCYLINDER");
149 7 : keys.add("atoms","CENTER","the atom whose vicinity we are interested in examining");
150 7 : keys.add("optional","SIGMA","the width of the function to be used for kernel density estimation");
151 7 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
152 7 : keys.add("compulsory","DIRECTION","the direction of the long axis of the cylinder. Must be x, y or z");
153 7 : keys.add("compulsory","RADIUS","a switching function that gives the extent of the cylinder in the plane perpendicular to the direction");
154 7 : keys.add("compulsory","LOWER","0.0","the lower boundary on the direction parallel to the long axis of the cylinder");
155 7 : keys.add("compulsory","UPPER","0.0","the upper boundary on the direction parallel to the long axis of the cylinder");
156 14 : keys.linkActionInDocs("RADIUS","LESS_THAN");
157 7 : }
158 :
159 1 : void VolumeInCylinder::parseInput( ActionVolume<VolumeInCylinder>* action ) {
160 2 : action->parse("SIGMA",sigma);
161 : std::string mykerneltype;
162 1 : action->parse("KERNEL",mykerneltype);
163 1 : kerneltype=HistogramBead::getKernelType(mykerneltype);
164 : std::string sdir;
165 2 : action->parse("DIRECTION",sdir);
166 1 : if( sdir=="X") {
167 0 : dir[0]=1;
168 0 : dir[1]=2;
169 0 : dir[2]=0;
170 1 : } else if( sdir=="Y") {
171 0 : dir[0]=0;
172 0 : dir[1]=2;
173 0 : dir[2]=1;
174 1 : } else if( sdir=="Z") {
175 1 : dir[0]=0;
176 1 : dir[1]=1;
177 1 : dir[2]=2;
178 : } else {
179 0 : action->error(sdir + "is not a valid direction. Should be X, Y or Z");
180 : }
181 1 : action->log.printf(" cylinder's long axis is along %s axis\n",sdir.c_str() );
182 :
183 : std::string errors;
184 : std::string swinput;
185 2 : action->parse("RADIUS",swinput);
186 1 : if(swinput.length()==0) {
187 0 : action->error("missing RADIUS keyword");
188 : }
189 1 : switchingFunction.set(swinput,errors);
190 1 : if( errors.length()!=0 ) {
191 0 : action->error("problem reading RADIUS keyword : " + errors );
192 : }
193 1 : action->log.printf(" radius of cylinder is given by %s \n", ( switchingFunction.description() ).c_str() );
194 :
195 1 : docylinder=false;
196 1 : action->parse("LOWER",min);
197 1 : action->parse("UPPER",max);
198 1 : if( min!=0.0 || max!=0.0 ) {
199 1 : if( min>max ) {
200 0 : action->error("minimum of cylinder should be less than maximum");
201 : }
202 1 : docylinder=true;
203 1 : action->log.printf(" cylinder extends from %f to %f along the %s axis\n",min,max,sdir.c_str() );
204 : }
205 1 : }
206 :
207 1 : void VolumeInCylinder::parseAtoms( ActionVolume<VolumeInCylinder>* action, std::vector<AtomNumber>& atom ) {
208 2 : action->parseAtomList("CENTER",atom);
209 1 : if( atom.size()!=1 ) {
210 0 : action->error("should only be one atom specified");
211 : }
212 1 : action->log.printf(" center of cylinder is at position of atom : %d\n",atom[0].serial() );
213 1 : }
214 :
215 8000 : void VolumeInCylinder::calculateNumberInside( const VolumeInput& input, const VolumeInCylinder& actioninput, VolumeOutput& output ) {
216 : // Calculate position of atom wrt to origin
217 8000 : Vector fpos=input.pbc.distance( Vector(input.refpos[0][0],input.refpos[0][1],input.refpos[0][2]), Vector(input.cpos[0],input.cpos[1],input.cpos[2]) );
218 :
219 : double vcylinder, dcylinder;
220 8000 : if( actioninput.docylinder ) {
221 8000 : HistogramBead bead( actioninput.kerneltype,
222 8000 : actioninput.min, actioninput.max, actioninput.sigma );
223 8000 : vcylinder=bead.calculate( fpos[actioninput.dir[2]], dcylinder );
224 : } else {
225 : vcylinder=1.0;
226 0 : dcylinder=0.0;
227 : }
228 :
229 8000 : const double dd = fpos[actioninput.dir[0]]*fpos[actioninput.dir[0]] + fpos[actioninput.dir[1]]*fpos[actioninput.dir[1]];
230 : double dfunc;
231 8000 : double vswitch = actioninput.switchingFunction.calculateSqr( dd, dfunc );
232 8000 : output.values[0]=vswitch*vcylinder;
233 8000 : output.derivatives[actioninput.dir[0]]=vcylinder*dfunc*fpos[actioninput.dir[0]];
234 8000 : output.derivatives[actioninput.dir[1]]=vcylinder*dfunc*fpos[actioninput.dir[1]];
235 8000 : output.derivatives[actioninput.dir[2]]=vswitch*dcylinder;
236 : // Add derivatives wrt to position of origin atom
237 8000 : output.refders[0][0] = -output.derivatives[0];
238 8000 : output.refders[0][1] = -output.derivatives[1];
239 8000 : output.refders[0][2] = -output.derivatives[2];
240 : // Add virial contribution
241 16000 : output.virial.set( 0, -Tensor(fpos,Vector(output.derivatives[0], output.derivatives[1], output.derivatives[2])) );
242 8000 : }
243 :
244 : }
245 : }
|
