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 ATOMS=14401-74134:3 FIELD_ATOMS=1-14400 CONTOUR={RATIONAL D_0=2.0 R_0=1.0} BANDWIDTH=0.1,0.1,0.1
53 : PRINT ARG=fi FILE=colvar
54 : ```
55 :
56 : 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
57 : 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
58 : 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
59 : values for a symmetry function for those atoms that are within the membrane.
60 :
61 : */
62 : //+ENDPLUMEDOC
63 :
64 : //+PLUMEDOC VOLUMES INENVELOPE_CALC
65 : /*
66 : 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.
67 :
68 : \par Examples
69 :
70 : */
71 : //+ENDPLUMEDOC
72 :
73 : namespace PLMD {
74 : namespace volumes {
75 :
76 : class VolumeInEnvelope : public ActionVolume {
77 : private:
78 : LinkCells mylinks;
79 : double gvol;
80 : std::vector<std::unique_ptr<Value>> pos;
81 : std::vector<Vector> ltmp_pos;
82 : std::vector<unsigned> ltmp_ind;
83 : std::vector<double> bandwidth;
84 : SwitchingFunction sfunc, switchingFunction;
85 : public:
86 : static void registerKeywords( Keywords& keys );
87 : explicit VolumeInEnvelope(const ActionOptions& ao);
88 : void setupRegions() override;
89 : double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const override;
90 : };
91 :
92 : PLUMED_REGISTER_ACTION(VolumeInEnvelope,"INENVELOPE_CALC")
93 : char glob_contours[] = "INENVELOPE";
94 : typedef VolumeShortcut<glob_contours> VolumeInEnvelopeShortcut;
95 : PLUMED_REGISTER_ACTION(VolumeInEnvelopeShortcut,"INENVELOPE")
96 :
97 7 : void VolumeInEnvelope::registerKeywords( Keywords& keys ) {
98 7 : ActionVolume::registerKeywords( keys );
99 7 : keys.remove("SIGMA");
100 14 : keys.setDisplayName("INENVELOPE");
101 7 : keys.add("atoms","FIELD_ATOMS","the atom whose positions we are constructing a field from");
102 7 : keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density esimtation");
103 7 : keys.add("compulsory","CONTOUR","a switching funciton that tells PLUMED how large the density should be");
104 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");
105 7 : }
106 :
107 1 : VolumeInEnvelope::VolumeInEnvelope(const ActionOptions& ao):
108 : Action(ao),
109 : ActionVolume(ao),
110 1 : mylinks(comm) {
111 : std::vector<AtomNumber> atoms;
112 1 : parseAtomList("FIELD_ATOMS",atoms);
113 1 : log.printf(" creating density field from atoms : ");
114 9 : for(unsigned i=0; i<atoms.size(); ++i) {
115 8 : log.printf("%d ",atoms[i].serial() );
116 : }
117 1 : log.printf("\n");
118 1 : ltmp_ind.resize( atoms.size() );
119 1 : ltmp_pos.resize( atoms.size() );
120 9 : for(unsigned i=0; i<atoms.size(); ++i) {
121 8 : ltmp_ind[i]=i;
122 : }
123 :
124 : std::string sw, errors;
125 2 : parse("CONTOUR",sw);
126 1 : if(sw.length()==0) {
127 0 : error("missing CONTOUR keyword");
128 : }
129 1 : sfunc.set(sw,errors);
130 1 : if( errors.length()!=0 ) {
131 0 : error("problem reading RADIUS keyword : " + errors );
132 : }
133 1 : log.printf(" density at atom must be larger than %s \n", ( sfunc.description() ).c_str() );
134 :
135 1 : std::vector<double> pp(3,0.0);
136 1 : bandwidth.resize(3);
137 1 : parseVector("BANDWIDTH",bandwidth);
138 2 : log.printf(" using %s kernel with bandwidths %f %f %f \n",getKernelType().c_str(),bandwidth[0],bandwidth[1],bandwidth[2] );
139 : std::string errors2;
140 2 : switchingFunction.set("GAUSSIAN R_0=1.0 NOSTRETCH", errors2 );
141 1 : if( errors2.length()!=0 ) {
142 0 : error("problem reading switching function description " + errors2);
143 : }
144 : double det=1;
145 4 : for(unsigned i=0; i<bandwidth.size(); ++i) {
146 3 : det*=bandwidth[i]*bandwidth[i];
147 : }
148 1 : gvol=1.0;
149 1 : if( getKernelType()=="gaussian" ) {
150 1 : gvol=pow( 2*pi, 0.5*bandwidth.size() ) * pow( det, 0.5 );
151 : } else {
152 0 : error("cannot use kernel other than gaussian");
153 : }
154 : double dp2cutoff;
155 1 : parse("CUTOFF",dp2cutoff);
156 1 : double maxs = sqrt(2*dp2cutoff)*bandwidth[0];
157 3 : for(unsigned j=1; j<bandwidth.size(); ++j) {
158 2 : if( sqrt(2*dp2cutoff)*bandwidth[j]>maxs ) {
159 0 : maxs=sqrt(2*dp2cutoff)*bandwidth[j];
160 : }
161 : }
162 1 : checkRead();
163 1 : requestAtoms(atoms);
164 1 : mylinks.setCutoff( maxs );
165 1 : }
166 :
167 5 : void VolumeInEnvelope::setupRegions() {
168 45 : for(unsigned i=0; i<ltmp_ind.size(); ++i) {
169 40 : ltmp_pos[i]=getPosition(i);
170 : }
171 5 : mylinks.buildCellLists( ltmp_pos, ltmp_ind, getPbc() );
172 5 : }
173 :
174 1000 : double VolumeInEnvelope::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const {
175 1000 : unsigned ncells_required=0, natoms=1;
176 1000 : std::vector<unsigned> cells_required( mylinks.getNumberOfCells() ), indices( 1 + getNumberOfAtoms() );
177 1000 : mylinks.addRequiredCells( mylinks.findMyCell( cpos ), ncells_required, cells_required );
178 1000 : indices[0]=getNumberOfAtoms();
179 1000 : mylinks.retrieveAtomsInCells( ncells_required, cells_required, natoms, indices );
180 : double value=0;
181 1000 : std::vector<double> der(3);
182 1000 : Vector tder;
183 :
184 : // convert pointer once
185 1000 : auto pos_ptr=Tools::unique2raw(pos);
186 9000 : for(unsigned i=1; i<natoms; ++i) {
187 8000 : Vector dist = pbcDistance( cpos, getPosition( indices[i] ) );
188 : double dval=0;
189 32000 : for(unsigned j=0; j<3; ++j) {
190 24000 : der[j] = dist[j]/bandwidth[j];
191 24000 : dval += der[j]*der[j];
192 24000 : der[j] = der[j] / bandwidth[j];
193 : }
194 : double dfunc;
195 8000 : value += switchingFunction.calculateSqr( dval, dfunc ) / gvol;
196 8000 : double tmp = dfunc / gvol;
197 32000 : for(unsigned j=0; j<3; ++j) {
198 24000 : derivatives[j] -= tmp*der[j];
199 24000 : refders[ indices[i] ][j] += tmp*der[j];
200 24000 : tder[j]=tmp*der[j];
201 : }
202 8000 : vir -= Tensor( tder, dist );
203 : }
204 1000 : double deriv, fval = sfunc.calculate( value, deriv );
205 1000 : derivatives *= -deriv*value;
206 1000 : vir *= -deriv*value;
207 9000 : for(unsigned i=1; i<natoms; ++i) {
208 8000 : refders[ indices[i] ] *= -deriv*value;
209 : }
210 2000 : return 1.0 - fval;
211 : }
212 :
213 : }
214 : }
|
