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 : Each of the base quantities calculated by a multicolvar can can be assigned to a particular point in three
33 : dimensional space. For example, if we have the coordination numbers for all the atoms in the
34 : system each coordination number can be assumed to lie on the position of the central atom.
35 : Because each base quantity can be assigned to a particular point in space we can calculate functions of the
36 : distribution of base quantities in a particular part of the box by using:
37 :
38 : \f[
39 : \overline{s}_{\tau} = \frac{ \sum_i f(s_i) \sigma(r_{xy}) }{ \sum_i \sigma(r_{xy}) }
40 : \f]
41 :
42 : where the sum is over the collective variables, \f$s_i\f$, each of which can be thought to be at \f$ (x_i,y_i,z_i)\f$.
43 : The function \f$\sigma\f$ is a \ref switchingfunction that acts on the distance between the point at which the
44 : collective is located \f$(x_i,y_i,z_i)\f$ and the position of the atom that was specified using the ORIGIN keyword
45 : projected in the xy plane if DIRECTION=z is used. In other words:
46 : \f[
47 : r_{xy} = sqrt{ ( x_i - x_0)^2 + ( y_i - y_0)^2 }
48 : \f]
49 : In short this function, \f$\sigma(r_{xy})\f$, measures whether or not the CV is within a cylinder that
50 : runs along the axis specified using the DIRECTION keyword and that is centered on the position of the atom specified using
51 : ORIGIN.
52 :
53 : The function \f$(s_i)\f$ can be any of the usual LESS_THAN, MORE_THAN, WITHIN etc that are used in all other multicolvars.
54 :
55 : When INCYLINDER is used with the \ref DENSITY action the number of atoms in the specified region is calculated
56 :
57 : \par Examples
58 :
59 : The input below can be use to calculate the average coordination numbers for those atoms that are within a cylindrical tube
60 : of radius 1.5 nm that is centered on the position of atom 101 and that has its long axis parallel to the z-axis.
61 :
62 : \plumedfile
63 : c1: COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=0.1}
64 : d2: INCYLINDER ATOM=101 DATA=c1 DIRECTION=Z RADIUS={TANH R_0=1.5} SIGMA=0.1 LOWER=-0.1 UPPER=0.1 MEAN
65 : PRINT ARG=d2.* FILE=colvar
66 : \endplumedfile
67 :
68 : */
69 : //+ENDPLUMEDOC
70 :
71 : //+PLUMEDOC MCOLVAR INCYLINDER_CALC
72 : /*
73 : 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
74 :
75 : \par Examples
76 :
77 : */
78 : //+ENDPLUMEDOC
79 :
80 : namespace PLMD {
81 : namespace volumes {
82 :
83 : class VolumeInCylinder : public ActionVolume {
84 : private:
85 : bool docylinder;
86 : Vector origin;
87 : HistogramBead bead;
88 : std::vector<unsigned> dir;
89 : SwitchingFunction switchingFunction;
90 : public:
91 : static void registerKeywords( Keywords& keys );
92 : explicit VolumeInCylinder (const ActionOptions& ao);
93 : void setupRegions() override;
94 : double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const override;
95 : };
96 :
97 : PLUMED_REGISTER_ACTION(VolumeInCylinder,"INCYLINDER_CALC")
98 : char glob_cylinder[] = "INCYLINDER";
99 : typedef VolumeShortcut<glob_cylinder> VolumeInCylinderShortcut;
100 : PLUMED_REGISTER_ACTION(VolumeInCylinderShortcut,"INCYLINDER")
101 :
102 7 : void VolumeInCylinder::registerKeywords( Keywords& keys ) {
103 7 : ActionVolume::registerKeywords( keys );
104 7 : keys.setDisplayName("INCYLINDER");
105 14 : keys.add("atoms","CENTER","the atom whose vicinity we are interested in examining");
106 14 : keys.add("compulsory","DIRECTION","the direction of the long axis of the cylinder. Must be x, y or z");
107 14 : keys.add("compulsory","RADIUS","a switching function that gives the extent of the cylinder in the plane perpendicular to the direction");
108 14 : keys.add("compulsory","LOWER","0.0","the lower boundary on the direction parallel to the long axis of the cylinder");
109 14 : keys.add("compulsory","UPPER","0.0","the upper boundary on the direction parallel to the long axis of the cylinder");
110 14 : keys.reset_style("SIGMA","optional");
111 7 : }
112 :
113 1 : VolumeInCylinder::VolumeInCylinder(const ActionOptions& ao):
114 : Action(ao),
115 : ActionVolume(ao),
116 1 : docylinder(false) {
117 : std::vector<AtomNumber> atom;
118 2 : parseAtomList("CENTER",atom);
119 1 : if( atom.size()!=1 ) {
120 0 : error("should only be one atom specified");
121 : }
122 1 : log.printf(" center of cylinder is at position of atom : %d\n",atom[0].serial() );
123 :
124 : std::string sdir;
125 2 : parse("DIRECTION",sdir);
126 1 : if( sdir=="X") {
127 0 : dir.push_back(1);
128 0 : dir.push_back(2);
129 0 : dir.push_back(0);
130 1 : } else if( sdir=="Y") {
131 0 : dir.push_back(0);
132 0 : dir.push_back(2);
133 0 : dir.push_back(1);
134 1 : } else if( sdir=="Z") {
135 1 : dir.push_back(0);
136 1 : dir.push_back(1);
137 1 : dir.push_back(2);
138 : } else {
139 0 : error(sdir + "is not a valid direction. Should be X, Y or Z");
140 : }
141 1 : log.printf(" cylinder's long axis is along %s axis\n",sdir.c_str() );
142 :
143 : std::string sw, errors;
144 2 : parse("RADIUS",sw);
145 1 : if(sw.length()==0) {
146 0 : error("missing RADIUS keyword");
147 : }
148 1 : switchingFunction.set(sw,errors);
149 1 : if( errors.length()!=0 ) {
150 0 : error("problem reading RADIUS keyword : " + errors );
151 : }
152 1 : log.printf(" radius of cylinder is given by %s \n", ( switchingFunction.description() ).c_str() );
153 :
154 : double min, max;
155 1 : parse("LOWER",min);
156 1 : parse("UPPER",max);
157 1 : if( min!=0.0 || max!=0.0 ) {
158 1 : if( min>max ) {
159 0 : error("minimum of cylinder should be less than maximum");
160 : }
161 1 : docylinder=true;
162 1 : log.printf(" cylinder extends from %f to %f along the %s axis\n",min,max,sdir.c_str() );
163 : bead.isNotPeriodic();
164 2 : bead.setKernelType( getKernelType() );
165 1 : bead.set( min, max, getSigma() );
166 : }
167 :
168 1 : checkRead();
169 1 : requestAtoms(atom);
170 1 : }
171 :
172 20 : void VolumeInCylinder::setupRegions() { }
173 :
174 4667 : double VolumeInCylinder::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const {
175 : // Calculate position of atom wrt to origin
176 4667 : Vector fpos=pbcDistance( getPosition(0), cpos );
177 :
178 : double vcylinder, dcylinder;
179 4667 : if( docylinder ) {
180 4667 : vcylinder=bead.calculate( fpos[dir[2]], dcylinder );
181 : } else {
182 : vcylinder=1.0;
183 0 : dcylinder=0.0;
184 : }
185 :
186 4667 : const double dd = fpos[dir[0]]*fpos[dir[0]] + fpos[dir[1]]*fpos[dir[1]];
187 4667 : double dfunc, vswitch = switchingFunction.calculateSqr( dd, dfunc );
188 4667 : derivatives.zero();
189 4667 : double value=vswitch*vcylinder;
190 4667 : derivatives[dir[0]]=vcylinder*dfunc*fpos[dir[0]];
191 4667 : derivatives[dir[1]]=vcylinder*dfunc*fpos[dir[1]];
192 4667 : derivatives[dir[2]]=vswitch*dcylinder;
193 : // Add derivatives wrt to position of origin atom
194 4667 : refders[0] = -derivatives;
195 : // Add virial contribution
196 4667 : vir -= Tensor(fpos,derivatives);
197 4667 : return value;
198 : }
199 :
200 : }
201 : }
|
