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 : using precision=double;
94 : typedef VolumeIn<precision> VolumeInput;
95 : typedef VolumeOut<precision> VolumeOutput;
96 : double gvol=0.0;
97 : double maxs=0.0;
98 : std::string kerneltype;
99 : std::vector<double> bandwidth;
100 : std::string sfunc_str;
101 : SwitchingFunction sfunc, switchingFunction;
102 : unsigned natoms_in_list=0;
103 : unsigned natoms_per_list=0;
104 : std::vector<std::size_t> nlist;
105 : static void registerKeywords( Keywords& keys );
106 : void parseInput( ActionVolume<VolumeInEnvelope>* action );
107 : void setupRegions( ActionVolume<VolumeInEnvelope>* action, const Pbc& pbc, const std::vector<Vector>& positions );
108 : static void parseAtoms( ActionVolume<VolumeInEnvelope>* action, std::vector<AtomNumber>& atom );
109 : // the assignment is customized (see below): copy construction is the default one
110 2 : VolumeInEnvelope() = default;
111 : VolumeInEnvelope(const VolumeInEnvelope&) = default;
112 1 : VolumeInEnvelope& operator=( const VolumeInEnvelope& m ) {
113 1 : gvol=m.gvol;
114 1 : maxs=m.maxs;
115 1 : kerneltype=m.kerneltype;
116 1 : bandwidth=m.bandwidth;
117 1 : sfunc_str=m.sfunc_str;
118 : std::string errors;
119 1 : sfunc.set(sfunc_str, errors);
120 1 : switchingFunction.set("GAUSSIAN R_0=1.0 NOSTRETCH", errors );
121 1 : natoms_in_list = m.natoms_in_list;
122 1 : natoms_per_list = m.natoms_per_list;
123 1 : return *this;
124 : }
125 : static void calculateNumberInside( const VolumeInput& input, const VolumeInEnvelope& actioninput, VolumeOutput& output );
126 : };
127 :
128 : typedef ActionVolume<VolumeInEnvelope> volenv;
129 : PLUMED_REGISTER_ACTION(volenv,"INENVELOPE_CALC")
130 : char glob_contours[] = "INENVELOPE";
131 : typedef VolumeShortcut<glob_contours> VolumeInEnvelopeShortcut;
132 : PLUMED_REGISTER_ACTION(VolumeInEnvelopeShortcut,"INENVELOPE")
133 :
134 11 : void VolumeInEnvelope::registerKeywords( Keywords& keys ) {
135 22 : keys.setDisplayName("INENVELOPE");
136 11 : keys.add("atoms","FIELD_ATOMS","the atom whose positions we are constructing a field from");
137 11 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
138 11 : keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density esimtation");
139 11 : keys.add("compulsory","CONTOUR","a switching function that tells PLUMED how large the density should be");
140 11 : 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");
141 11 : }
142 :
143 1 : void VolumeInEnvelope::parseInput( ActionVolume<VolumeInEnvelope>* action ) {
144 2 : action->parse("KERNEL",kerneltype);
145 :
146 : std::string errors;
147 2 : action->parse("CONTOUR",sfunc_str);
148 1 : if(sfunc_str.length()==0) {
149 0 : action->error("missing CONTOUR keyword");
150 : }
151 1 : sfunc.set(sfunc_str,errors);
152 1 : if( errors.length()!=0 ) {
153 0 : action->error("problem reading RADIUS keyword : " + errors );
154 : }
155 1 : action->log.printf(" density at atom must be larger than %s \n", ( sfunc.description() ).c_str() );
156 :
157 1 : std::vector<double> pp(3,0.0);
158 1 : bandwidth.resize(3);
159 1 : action->parseVector("BANDWIDTH",bandwidth);
160 1 : action->log.printf(" using %s kernel with bandwidths %f %f %f \n",kerneltype.c_str(),bandwidth[0],bandwidth[1],bandwidth[2] );
161 : std::string errors2;
162 2 : switchingFunction.set("GAUSSIAN R_0=1.0 NOSTRETCH", errors2 );
163 1 : if( errors2.length()!=0 ) {
164 0 : action->error("problem reading switching function description " + errors2);
165 : }
166 : double det=1;
167 4 : for(unsigned i=0; i<bandwidth.size(); ++i) {
168 3 : det*=bandwidth[i]*bandwidth[i];
169 : }
170 1 : gvol=1.0;
171 1 : if( kerneltype=="gaussian" ) {
172 1 : gvol=pow( 2*pi, 0.5*bandwidth.size() ) * pow( det, 0.5 );
173 : } else {
174 0 : action->error("cannot use kernel other than gaussian");
175 : }
176 : double dp2cutoff;
177 1 : action->parse("CUTOFF",dp2cutoff);
178 1 : maxs = sqrt(2*dp2cutoff)*bandwidth[0];
179 3 : for(unsigned j=1; j<bandwidth.size(); ++j) {
180 2 : if( sqrt(2*dp2cutoff)*bandwidth[j]>maxs ) {
181 0 : maxs=sqrt(2*dp2cutoff)*bandwidth[j];
182 : }
183 : }
184 2 : }
185 :
186 1 : void VolumeInEnvelope::parseAtoms( ActionVolume<VolumeInEnvelope>* action, std::vector<AtomNumber>& atoms ) {
187 1 : action->parseAtomList("FIELD_ATOMS",atoms);
188 1 : action->log.printf(" creating density field from atoms : ");
189 9 : for(unsigned i=0; i<atoms.size(); ++i) {
190 8 : action->log.printf("%d ",atoms[i].serial() );
191 : }
192 1 : action->log.printf("\n");
193 1 : }
194 :
195 5 : void VolumeInEnvelope::setupRegions( ActionVolume<VolumeInEnvelope>* action, const Pbc& pbc, const std::vector<Vector>& positions ) {
196 5 : LinkCells mylinks(action->comm);
197 5 : mylinks.setCutoff( maxs );
198 5 : std::vector<unsigned> ltmp_ind(positions.size());
199 45 : for(unsigned i=0; i<ltmp_ind.size(); ++i) {
200 40 : ltmp_ind[i]=i;
201 : }
202 5 : natoms_in_list = (action->copyOutput(0))->getShape()[0];
203 5 : std::vector<unsigned> ind( natoms_in_list );
204 5 : std::vector<unsigned> tind( natoms_in_list );
205 5 : std::vector<Vector> volpos( natoms_in_list );
206 505 : for(unsigned i=0; i<natoms_in_list; ++i) {
207 500 : tind[i] = i;
208 500 : ind[i] = positions.size() + i;
209 500 : volpos[i] = action->getPosition(i);
210 : }
211 5 : mylinks.createNeighborList( natoms_in_list,
212 : make_const_view(volpos),
213 : make_const_view(ind),
214 : make_const_view(tind),
215 : make_const_view(positions),
216 : make_const_view(ltmp_ind),
217 : pbc,
218 5 : natoms_per_list,
219 5 : nlist );
220 5 : }
221 :
222 1000 : void VolumeInEnvelope::calculateNumberInside( const VolumeInput& input, const VolumeInEnvelope& actioninput, VolumeOutput& output ) {
223 : double value=0, dfunc;
224 1000 : std::vector<double> der(3);
225 1000 : Vector tder;
226 :
227 1000 : Tensor vir;
228 : vir.zero();
229 1000 : output.derivatives[0]=output.derivatives[1]=output.derivatives[2]=0;
230 1000 : std::size_t lstart = actioninput.natoms_in_list + input.task_index*(1+actioninput.natoms_per_list);
231 9000 : for(unsigned i=1; i<actioninput.nlist[input.task_index]; ++i) {
232 8000 : unsigned atno = actioninput.nlist[lstart+i];
233 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]) );
234 : double dval=0;
235 32000 : for(unsigned j=0; j<3; ++j) {
236 24000 : der[j] = dist[j]/actioninput.bandwidth[j];
237 24000 : dval += der[j]*der[j];
238 24000 : der[j] = der[j] / actioninput.bandwidth[j];
239 : }
240 8000 : value += actioninput.switchingFunction.calculateSqr( dval, dfunc ) / actioninput.gvol;
241 8000 : double tmp = dfunc / actioninput.gvol;
242 32000 : for(unsigned j=0; j<3; ++j) {
243 24000 : output.derivatives[j] -= tmp*der[j];
244 24000 : output.refders[atno][j] = tmp*der[j];
245 24000 : tder[j]=tmp*der[j];
246 : }
247 8000 : vir -= Tensor( tder, dist );
248 : }
249 : double deriv;
250 1000 : output.values[0] = 1 - actioninput.sfunc.calculate( value, deriv );
251 1000 : output.derivatives[0] *= -deriv*value;
252 1000 : output.derivatives[1] *= -deriv*value;
253 1000 : output.derivatives[2] *= -deriv*value;
254 1000 : output.virial.set( 0, -deriv*value*vir );
255 9000 : for(unsigned i=1; i<actioninput.nlist[input.task_index]; ++i) {
256 8000 : unsigned atno = actioninput.nlist[lstart+i];
257 8000 : output.refders[ atno ][0] *= -deriv*value;
258 8000 : output.refders[ atno ][1] *= -deriv*value;
259 8000 : output.refders[ atno ][2] *= -deriv*value;
260 : }
261 9000 : }
262 :
263 : }
264 : }
|
