Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-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 "tools/LinkCells.h"
26 : #include "ActionVolume.h"
27 : #include "VolumeShortcut.h"
28 : #include <memory>
29 :
30 : //+PLUMEDOC VOLUMES INENVELOPE
31 : /*
32 : This quantity can be used to calculate functions of the distribution of collective variables for the atoms that lie in a region where the density of a certain type of atom is high.
33 :
34 : This collective variable can be used to determine whether atoms are within region where the density
35 : of a particular atom type is high. This is achieved by calculating the following function at the point where
36 : each atom is located $(x_i,y_i,z_i)$:
37 :
38 : $$
39 : w_i = 1 - \sigma\left[ \sum_{i=1}^N K\left( \frac{x-x_i}{\sigma_x},\frac{y-y_i}{\sigma_y},\frac{z-z_i}{\sigma_z} \right) \right]
40 : $$
41 :
42 : Here $\sigma$ is one of the switching functions that is discussed in the documentation for the action [LESS_THAN](LESS_THAN.md) and $K$ is
43 : one of the kernel functions that is discussed in the documentation for the action [BETWEEN](BETWEEN.md). The sum runs over the atoms
44 : specified using the FIELD_ATOMS keyword and a $w_j$ value is calculated for each of the central atoms of the input
45 : multicolvar.
46 :
47 : The input below shows how this action works in practise. This input calculates a density field from the positions of atoms 1-14400. A vector which has
48 : as many elements as atoms that were specified using the ATOMS keyword. The $i$th element of this vector is calculated using the expression above with $(x_i,y_i,z_i)$
49 : being the position of the $i$th atom that was specified using that ATOMS keyword.
50 :
51 : ```plumed
52 : fi: INENVELOPE ...
53 : ATOMS=14401-74134:3 FIELD_ATOMS=1-14400
54 : CONTOUR={RATIONAL D_0=2.0 R_0=1.0}
55 : KERNEL=gaussian BANDWIDTH=0.1,0.1,0.1
56 : CUTOFF=6.25
57 : ...
58 : PRINT ARG=fi FILE=colvar
59 : ```
60 :
61 : This particular action was developed with the intention of determining whether water molecules had penetrated a membrane or not. The FIELD_ATOMS were thus the atoms of the
62 : lipid molecules that made up the membrane and the ATOMS were the oxygens of the water molecules. The vector that is output by this action can be used in all the ways that the
63 : vector that is output by the [AROUND](AROUND.md) action is used. In other words, this action can be used to calculate the number of water molecules in the membrane or the average
64 : values for a symmetry function for those atoms that are within the membrane. You can also use this action to calculate the number of atoms that are not in the membrane by using
65 : an input like the one shown below:
66 :
67 : ```plumed
68 : fi: INENVELOPE ...
69 : ATOMS=14401-74134:3 FIELD_ATOMS=1-14400
70 : CONTOUR={RATIONAL D_0=2.0 R_0=1.0}
71 : BANDWIDTH=0.1,0.1,0.1
72 : OUTSIDE
73 : ...
74 : s: SUM ARG=fi PERIODIC=NO
75 : PRINT ARG=s FILE=colvar
76 : ```
77 :
78 : !!! note ""
79 :
80 : As with [AROUND](AROUND.md) there was syntax for caclulating the average values of order parameters for those atoms that are inside/outside the membrane, which can
81 : still be used with this new version of the command. However, the same calculations can be performed in later versions of the code with a better syntax. We strongly
82 : recommend using the newer syntax but if you are interested in the old syntax you can find more information in the old syntax section of the documentation for [AROUND](AROUND.md).
83 : The documentation for that action tells you how that old syntax worked and how you can achieve the same results using the new syntax.
84 :
85 : */
86 : //+ENDPLUMEDOC
87 :
88 : namespace PLMD {
89 : namespace volumes {
90 :
91 : class VolumeInEnvelope {
92 : public:
93 : double gvol, maxs;
94 : std::string kerneltype;
95 : std::vector<double> bandwidth;
96 : std::string sfunc_str;
97 : SwitchingFunction sfunc, switchingFunction;
98 : unsigned natoms_in_list;
99 : unsigned natoms_per_list;
100 : std::vector<std::size_t> nlist;
101 : static void registerKeywords( Keywords& keys );
102 : void parseInput( ActionVolume<VolumeInEnvelope>* action );
103 : void setupRegions( ActionVolume<VolumeInEnvelope>* action, const Pbc& pbc, const std::vector<Vector>& positions );
104 : static void parseAtoms( ActionVolume<VolumeInEnvelope>* action, std::vector<AtomNumber>& atom );
105 1 : VolumeInEnvelope& operator=( const VolumeInEnvelope& m ) {
106 1 : gvol=m.gvol;
107 1 : maxs=m.maxs;
108 1 : kerneltype=m.kerneltype;
109 1 : bandwidth=m.bandwidth;
110 1 : sfunc_str=m.sfunc_str;
111 : std::string errors;
112 1 : sfunc.set(sfunc_str, errors);
113 1 : switchingFunction.set("GAUSSIAN R_0=1.0 NOSTRETCH", errors );
114 1 : natoms_in_list = m.natoms_in_list;
115 1 : natoms_per_list = m.natoms_per_list;
116 1 : return *this;
117 : }
118 : static void calculateNumberInside( const VolumeInput& input, const VolumeInEnvelope& actioninput, VolumeOutput& output );
119 : };
120 :
121 : typedef ActionVolume<VolumeInEnvelope> volenv;
122 : PLUMED_REGISTER_ACTION(volenv,"INENVELOPE_CALC")
123 : char glob_contours[] = "INENVELOPE";
124 : typedef VolumeShortcut<glob_contours> VolumeInEnvelopeShortcut;
125 : PLUMED_REGISTER_ACTION(VolumeInEnvelopeShortcut,"INENVELOPE")
126 :
127 7 : void VolumeInEnvelope::registerKeywords( Keywords& keys ) {
128 14 : keys.setDisplayName("INENVELOPE");
129 7 : keys.add("atoms","FIELD_ATOMS","the atom whose positions we are constructing a field from");
130 7 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
131 7 : keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density esimtation");
132 7 : keys.add("compulsory","CONTOUR","a switching funciton that tells PLUMED how large the density should be");
133 7 : keys.add("compulsory","CUTOFF","6.25","the cutoff at which to stop evaluating the kernel functions is set equal to sqrt(2*x)*bandwidth in each direction where x is this number");
134 7 : }
135 :
136 1 : void VolumeInEnvelope::parseInput( ActionVolume<VolumeInEnvelope>* action ) {
137 2 : action->parse("KERNEL",kerneltype);
138 :
139 : std::string errors;
140 2 : action->parse("CONTOUR",sfunc_str);
141 1 : if(sfunc_str.length()==0) {
142 0 : action->error("missing CONTOUR keyword");
143 : }
144 1 : sfunc.set(sfunc_str,errors);
145 1 : if( errors.length()!=0 ) {
146 0 : action->error("problem reading RADIUS keyword : " + errors );
147 : }
148 1 : action->log.printf(" density at atom must be larger than %s \n", ( sfunc.description() ).c_str() );
149 :
150 1 : std::vector<double> pp(3,0.0);
151 1 : bandwidth.resize(3);
152 1 : action->parseVector("BANDWIDTH",bandwidth);
153 1 : action->log.printf(" using %s kernel with bandwidths %f %f %f \n",kerneltype.c_str(),bandwidth[0],bandwidth[1],bandwidth[2] );
154 : std::string errors2;
155 2 : switchingFunction.set("GAUSSIAN R_0=1.0 NOSTRETCH", errors2 );
156 1 : if( errors2.length()!=0 ) {
157 0 : action->error("problem reading switching function description " + errors2);
158 : }
159 : double det=1;
160 4 : for(unsigned i=0; i<bandwidth.size(); ++i) {
161 3 : det*=bandwidth[i]*bandwidth[i];
162 : }
163 1 : gvol=1.0;
164 1 : if( kerneltype=="gaussian" ) {
165 1 : gvol=pow( 2*pi, 0.5*bandwidth.size() ) * pow( det, 0.5 );
166 : } else {
167 0 : action->error("cannot use kernel other than gaussian");
168 : }
169 : double dp2cutoff;
170 1 : action->parse("CUTOFF",dp2cutoff);
171 1 : maxs = sqrt(2*dp2cutoff)*bandwidth[0];
172 3 : for(unsigned j=1; j<bandwidth.size(); ++j) {
173 2 : if( sqrt(2*dp2cutoff)*bandwidth[j]>maxs ) {
174 0 : maxs=sqrt(2*dp2cutoff)*bandwidth[j];
175 : }
176 : }
177 1 : }
178 :
179 1 : void VolumeInEnvelope::parseAtoms( ActionVolume<VolumeInEnvelope>* action, std::vector<AtomNumber>& atoms ) {
180 1 : action->parseAtomList("FIELD_ATOMS",atoms);
181 1 : action->log.printf(" creating density field from atoms : ");
182 9 : for(unsigned i=0; i<atoms.size(); ++i) {
183 8 : action->log.printf("%d ",atoms[i].serial() );
184 : }
185 1 : action->log.printf("\n");
186 1 : }
187 :
188 5 : void VolumeInEnvelope::setupRegions( ActionVolume<VolumeInEnvelope>* action, const Pbc& pbc, const std::vector<Vector>& positions ) {
189 5 : LinkCells mylinks(action->comm);
190 5 : mylinks.setCutoff( maxs );
191 5 : std::vector<unsigned> ltmp_ind(positions.size());
192 45 : for(unsigned i=0; i<ltmp_ind.size(); ++i) {
193 40 : ltmp_ind[i]=i;
194 : }
195 5 : natoms_in_list = (action->copyOutput(0))->getShape()[0];
196 5 : std::vector<unsigned> ind( natoms_in_list );
197 5 : std::vector<unsigned> tind( natoms_in_list );
198 5 : std::vector<Vector> volpos( natoms_in_list );
199 505 : for(unsigned i=0; i<natoms_in_list; ++i) {
200 500 : tind[i] = i;
201 500 : ind[i] = positions.size() + i;
202 500 : volpos[i] = action->getPosition(i);
203 : }
204 5 : mylinks.createNeighborList( natoms_in_list,
205 : make_const_view(volpos),
206 : make_const_view(ind),
207 : make_const_view(tind),
208 : make_const_view(positions),
209 : make_const_view(ltmp_ind),
210 : pbc,
211 5 : natoms_per_list,
212 5 : nlist );
213 5 : }
214 :
215 1000 : void VolumeInEnvelope::calculateNumberInside( const VolumeInput& input, const VolumeInEnvelope& actioninput, VolumeOutput& output ) {
216 : double value=0, dfunc;
217 1000 : std::vector<double> der(3);
218 : Vector tder;
219 :
220 : Tensor vir;
221 : vir.zero();
222 1000 : output.derivatives[0]=output.derivatives[1]=output.derivatives[2]=0;
223 1000 : std::size_t lstart = actioninput.natoms_in_list + input.task_index*(1+actioninput.natoms_per_list);
224 9000 : for(unsigned i=1; i<actioninput.nlist[input.task_index]; ++i) {
225 8000 : unsigned atno = actioninput.nlist[lstart+i];
226 8000 : Vector dist = input.pbc.distance( Vector(input.cpos[0],input.cpos[1],input.cpos[2]), Vector(input.refpos[atno][0], input.refpos[atno][1], input.refpos[atno][2]) );
227 : double dval=0;
228 32000 : for(unsigned j=0; j<3; ++j) {
229 24000 : der[j] = dist[j]/actioninput.bandwidth[j];
230 24000 : dval += der[j]*der[j];
231 24000 : der[j] = der[j] / actioninput.bandwidth[j];
232 : }
233 8000 : value += actioninput.switchingFunction.calculateSqr( dval, dfunc ) / actioninput.gvol;
234 8000 : double tmp = dfunc / actioninput.gvol;
235 32000 : for(unsigned j=0; j<3; ++j) {
236 24000 : output.derivatives[j] -= tmp*der[j];
237 24000 : output.refders[atno][j] = tmp*der[j];
238 24000 : tder[j]=tmp*der[j];
239 : }
240 16000 : vir -= Tensor( tder, dist );
241 : }
242 : double deriv;
243 1000 : output.values[0] = 1 - actioninput.sfunc.calculate( value, deriv );
244 1000 : output.derivatives[0] *= -deriv*value;
245 1000 : output.derivatives[1] *= -deriv*value;
246 1000 : output.derivatives[2] *= -deriv*value;
247 1000 : output.virial.set( 0, -deriv*value*vir );
248 9000 : for(unsigned i=1; i<actioninput.nlist[input.task_index]; ++i) {
249 8000 : unsigned atno = actioninput.nlist[lstart+i];
250 8000 : output.refders[ atno ][0] *= -deriv*value;
251 8000 : output.refders[ atno ][1] *= -deriv*value;
252 8000 : output.refders[ atno ][2] *= -deriv*value;
253 : }
254 1000 : }
255 :
256 : }
257 : }
|
