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 : Each of the base quantities calculated by a multicolvar can can be assigned to a particular point in three
32 : dimensional space. For example, if we have the coordination numbers for all the atoms in the
33 : system each coordination number can be assumed to lie on the position of the central atom.
34 : Because each base quantity can be assigned to a particular point in space we can calculate functions of the
35 : distribution of base quantities in a particular part of the box by using:
36 :
37 : \f[
38 : \overline{s}_{\tau} = \frac{ \sum_i f(s_i) w(x_i,y_i,z_i) }{ \sum_i w(x_i,y_i,z_i) }
39 : \f]
40 :
41 : 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$.
42 : The function \f$ w(x_i,y_i,z_i) \f$ measures whether or not the system is in the subregion of interest. It
43 : is equal to:
44 :
45 : \f[
46 : 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)
47 : \f]
48 :
49 : where \f$K\f$ is one of the kernel functions described on \ref histogrambead and \f$\sigma\f$ is a bandwidth parameter.
50 : 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.
51 :
52 : When AROUND is used with the \ref DENSITY action the number of atoms in the specified region is calculated
53 :
54 : \par Examples
55 :
56 : The following commands tell plumed to calculate the average coordination number for the atoms
57 : that have x (in fractional coordinates) within 2.0 nm of the com of mass c1. The final value will be labeled s.mean.
58 : \plumedfile
59 : COM ATOMS=1-100 LABEL=c1
60 : COORDINATIONNUMBER SPECIES=1-100 R_0=1.0 LABEL=c
61 : AROUND DATA=c ATOM=c1 XLOWER=-2.0 XUPPER=2.0 SIGMA=0.1 MEAN LABEL=s
62 : \endplumedfile
63 :
64 : */
65 : //+ENDPLUMEDOC
66 :
67 : //+PLUMEDOC MCOLVAR AROUND_CALC
68 : /*
69 : 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
70 :
71 : \par Examples
72 :
73 : */
74 : //+ENDPLUMEDOC
75 :
76 : namespace PLMD {
77 : namespace volumes {
78 :
79 : class VolumeAround : public ActionVolume {
80 : private:
81 : Vector origin;
82 : bool dox, doy, doz;
83 : double xlow, xhigh;
84 : double ylow, yhigh;
85 : double zlow, zhigh;
86 : public:
87 : static void registerKeywords( Keywords& keys );
88 : explicit VolumeAround(const ActionOptions& ao);
89 : void setupRegions() override;
90 : double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const override;
91 : };
92 :
93 : PLUMED_REGISTER_ACTION(VolumeAround,"AROUND_CALC")
94 : char glob_around[] = "AROUND";
95 : typedef VolumeShortcut<glob_around> VolumeAroundShortcut;
96 : PLUMED_REGISTER_ACTION(VolumeAroundShortcut,"AROUND")
97 :
98 156 : void VolumeAround::registerKeywords( Keywords& keys ) {
99 156 : ActionVolume::registerKeywords( keys );
100 156 : keys.setDisplayName("AROUND");
101 312 : keys.add("atoms","ORIGIN","the atom whose vicinity we are interested in examining");
102 312 : keys.add("atoms-2","ATOM","an alternative to ORIGIN");
103 312 : 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).");
104 312 : 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).");
105 312 : 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).");
106 312 : 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).");
107 312 : 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).");
108 312 : 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).");
109 156 : }
110 :
111 46 : VolumeAround::VolumeAround(const ActionOptions& ao):
112 : Action(ao),
113 46 : ActionVolume(ao) {
114 : std::vector<AtomNumber> atom;
115 92 : parseAtomList("ORIGIN",atom);
116 46 : if( atom.size()==0 ) {
117 0 : parseAtomList("ATOM",atom);
118 : }
119 46 : if( atom.size()!=1 ) {
120 0 : error("should only be one atom specified");
121 : }
122 46 : log.printf(" boundaries for region are calculated based on positions of atom : %d\n",atom[0].serial() );
123 :
124 46 : dox=true;
125 46 : parse("XLOWER",xlow);
126 46 : parse("XUPPER",xhigh);
127 46 : doy=true;
128 46 : parse("YLOWER",ylow);
129 46 : parse("YUPPER",yhigh);
130 46 : doz=true;
131 46 : parse("ZLOWER",zlow);
132 46 : parse("ZUPPER",zhigh);
133 46 : if( xlow==0.0 && xhigh==0.0 ) {
134 6 : dox=false;
135 : }
136 46 : if( ylow==0.0 && yhigh==0.0 ) {
137 14 : doy=false;
138 : }
139 46 : if( zlow==0.0 && zhigh==0.0 ) {
140 14 : doz=false;
141 : }
142 46 : if( !dox && !doy && !doz ) {
143 0 : error("no subregion defined use XLOWER, XUPPER, YLOWER, YUPPER, ZLOWER, ZUPPER");
144 : }
145 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);
146 46 : checkRead();
147 46 : requestAtoms(atom);
148 46 : }
149 :
150 388 : void VolumeAround::setupRegions() { }
151 :
152 98837 : double VolumeAround::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const {
153 : // Setup the histogram bead
154 98837 : HistogramBead bead;
155 : bead.isNotPeriodic();
156 98837 : bead.setKernelType( getKernelType() );
157 :
158 : // Calculate position of atom wrt to origin
159 98837 : Vector fpos=pbcDistance( getPosition(0), cpos );
160 : double xcontr, ycontr, zcontr, xder, yder, zder;
161 98837 : if( dox ) {
162 66341 : bead.set( xlow, xhigh, getSigma() );
163 66341 : xcontr=bead.calculate( fpos[0], xder );
164 : } else {
165 : xcontr=1.;
166 32496 : xder=0.;
167 : }
168 98837 : if( doy ) {
169 63681 : bead.set( ylow, yhigh, getSigma() );
170 63681 : ycontr=bead.calculate( fpos[1], yder );
171 : } else {
172 : ycontr=1.;
173 35156 : yder=0.;
174 : }
175 98837 : if( doz ) {
176 62589 : bead.set( zlow, zhigh, getSigma() );
177 62589 : zcontr=bead.calculate( fpos[2], zder );
178 : } else {
179 : zcontr=1.;
180 36248 : zder=0.;
181 : }
182 98837 : derivatives[0]=xder*ycontr*zcontr;
183 98837 : derivatives[1]=xcontr*yder*zcontr;
184 98837 : derivatives[2]=xcontr*ycontr*zder;
185 : // Add derivatives wrt to position of origin atom
186 98837 : refders[0] = -derivatives;
187 : // Add virial contribution
188 98837 : vir -= Tensor(fpos,derivatives);
189 98837 : return xcontr*ycontr*zcontr;
190 : }
191 :
192 : }
193 : }
|
