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 : #include "core/ActionRegister.h"
23 : #include "tools/Pbc.h"
24 : #include "tools/SwitchingFunction.h"
25 : #include "ActionVolume.h"
26 : #include "VolumeShortcut.h"
27 :
28 : //+PLUMEDOC VOLUMES INCYLINDER
29 : /*
30 : 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.
31 :
32 : 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:
33 :
34 : ```plumed
35 : f: FIXEDATOM AT=0,0,0
36 : a: INCYLINDER ATOMS=1-100 CENTER=f DIRECTION=Z RADIUS={TANH R_0=1.5} SIGMA=0.1 LOWER=-1.0 UPPER=1.0
37 : PRINT ARG=a FILE=colvar
38 : ```
39 :
40 : The 100 elements of the vector `a` that is returned from the INCYLINDER action in the above input are calculated using:
41 :
42 : $$
43 : 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)
44 : $$
45 :
46 : where
47 :
48 : $$
49 : r_{xy} = \sqrt{ x_i^2 + y_i^2 }
50 : $$
51 :
52 : 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
53 : for the integrals are the values specified using the keywords `LOWER` and `UPPER`. $x_i$, $y_i$ and $z_i$ are then the components
54 : 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,
55 : $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
56 : 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
57 : specified using the `RADIUS` keyword.
58 :
59 : 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:
60 :
61 : ```plumed
62 : f: FIXEDATOM AT=0,0,0
63 : a: INCYLINDER ATOMS=1-100 CENTER=f DIRECTION=X RADIUS={TANH R_0=1.5} SIGMA=0.1 LOWER=-1.0 UPPER=1.0
64 : s: SUM ARG=a PERIODIC=NO
65 : PRINT ARG=s FILE=colvar
66 : ```
67 :
68 : 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.
69 :
70 : 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)
71 : 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
72 : types of calculations, which can still be used with this new version of the command. However, we strongly recommend using the newer syntax.
73 :
74 :
75 : */
76 : //+ENDPLUMEDOC
77 :
78 : //+PLUMEDOC MCOLVAR INCYLINDER_CALC
79 : /*
80 : 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
81 :
82 : \par Examples
83 :
84 : */
85 : //+ENDPLUMEDOC
86 :
87 : namespace PLMD {
88 : namespace volumes {
89 :
90 : class VolumeInCylinder : public ActionVolume {
91 : private:
92 : bool docylinder;
93 : Vector origin;
94 : HistogramBead bead;
95 : std::vector<unsigned> dir;
96 : SwitchingFunction switchingFunction;
97 : public:
98 : static void registerKeywords( Keywords& keys );
99 : explicit VolumeInCylinder (const ActionOptions& ao);
100 : void setupRegions() override;
101 : double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const override;
102 : };
103 :
104 : PLUMED_REGISTER_ACTION(VolumeInCylinder,"INCYLINDER_CALC")
105 : char glob_cylinder[] = "INCYLINDER";
106 : typedef VolumeShortcut<glob_cylinder> VolumeInCylinderShortcut;
107 : PLUMED_REGISTER_ACTION(VolumeInCylinderShortcut,"INCYLINDER")
108 :
109 7 : void VolumeInCylinder::registerKeywords( Keywords& keys ) {
110 7 : ActionVolume::registerKeywords( keys );
111 14 : keys.setDisplayName("INCYLINDER");
112 7 : keys.add("atoms","CENTER","the atom whose vicinity we are interested in examining");
113 7 : keys.add("compulsory","DIRECTION","the direction of the long axis of the cylinder. Must be x, y or z");
114 7 : keys.add("compulsory","RADIUS","a switching function that gives the extent of the cylinder in the plane perpendicular to the direction");
115 7 : keys.add("compulsory","LOWER","0.0","the lower boundary on the direction parallel to the long axis of the cylinder");
116 7 : keys.add("compulsory","UPPER","0.0","the upper boundary on the direction parallel to the long axis of the cylinder");
117 14 : keys.reset_style("SIGMA","optional");
118 14 : keys.linkActionInDocs("RADIUS","LESS_THAN");
119 7 : }
120 :
121 1 : VolumeInCylinder::VolumeInCylinder(const ActionOptions& ao):
122 : Action(ao),
123 : ActionVolume(ao),
124 1 : docylinder(false) {
125 : std::vector<AtomNumber> atom;
126 2 : parseAtomList("CENTER",atom);
127 1 : if( atom.size()!=1 ) {
128 0 : error("should only be one atom specified");
129 : }
130 1 : log.printf(" center of cylinder is at position of atom : %d\n",atom[0].serial() );
131 :
132 : std::string sdir;
133 2 : parse("DIRECTION",sdir);
134 1 : if( sdir=="X") {
135 0 : dir.push_back(1);
136 0 : dir.push_back(2);
137 0 : dir.push_back(0);
138 1 : } else if( sdir=="Y") {
139 0 : dir.push_back(0);
140 0 : dir.push_back(2);
141 0 : dir.push_back(1);
142 1 : } else if( sdir=="Z") {
143 1 : dir.push_back(0);
144 1 : dir.push_back(1);
145 1 : dir.push_back(2);
146 : } else {
147 0 : error(sdir + "is not a valid direction. Should be X, Y or Z");
148 : }
149 1 : log.printf(" cylinder's long axis is along %s axis\n",sdir.c_str() );
150 :
151 : std::string sw, errors;
152 2 : parse("RADIUS",sw);
153 1 : if(sw.length()==0) {
154 0 : error("missing RADIUS keyword");
155 : }
156 1 : switchingFunction.set(sw,errors);
157 1 : if( errors.length()!=0 ) {
158 0 : error("problem reading RADIUS keyword : " + errors );
159 : }
160 1 : log.printf(" radius of cylinder is given by %s \n", ( switchingFunction.description() ).c_str() );
161 :
162 : double min, max;
163 1 : parse("LOWER",min);
164 1 : parse("UPPER",max);
165 1 : if( min!=0.0 || max!=0.0 ) {
166 1 : if( min>max ) {
167 0 : error("minimum of cylinder should be less than maximum");
168 : }
169 1 : docylinder=true;
170 1 : log.printf(" cylinder extends from %f to %f along the %s axis\n",min,max,sdir.c_str() );
171 : bead.isNotPeriodic();
172 2 : bead.setKernelType( getKernelType() );
173 1 : bead.set( min, max, getSigma() );
174 : }
175 :
176 1 : checkRead();
177 1 : requestAtoms(atom);
178 1 : }
179 :
180 20 : void VolumeInCylinder::setupRegions() { }
181 :
182 4667 : double VolumeInCylinder::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const {
183 : // Calculate position of atom wrt to origin
184 4667 : Vector fpos=pbcDistance( getPosition(0), cpos );
185 :
186 : double vcylinder, dcylinder;
187 4667 : if( docylinder ) {
188 4667 : vcylinder=bead.calculate( fpos[dir[2]], dcylinder );
189 : } else {
190 : vcylinder=1.0;
191 0 : dcylinder=0.0;
192 : }
193 :
194 4667 : const double dd = fpos[dir[0]]*fpos[dir[0]] + fpos[dir[1]]*fpos[dir[1]];
195 4667 : double dfunc, vswitch = switchingFunction.calculateSqr( dd, dfunc );
196 4667 : derivatives.zero();
197 4667 : double value=vswitch*vcylinder;
198 4667 : derivatives[dir[0]]=vcylinder*dfunc*fpos[dir[0]];
199 4667 : derivatives[dir[1]]=vcylinder*dfunc*fpos[dir[1]];
200 4667 : derivatives[dir[2]]=vswitch*dcylinder;
201 : // Add derivatives wrt to position of origin atom
202 4667 : refders[0] = -derivatives;
203 : // Add virial contribution
204 4667 : vir -= Tensor(fpos,derivatives);
205 4667 : return value;
206 : }
207 :
208 : }
209 : }
|
