Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2023 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 "DistanceFromContourBase.h"
23 : #include "core/ActionRegister.h"
24 :
25 : //+PLUMEDOC COLVAR DISTANCE_FROM_CONTOUR
26 : /*
27 : Calculate the perpendicular distance from a Willard-Chandler dividing surface.
28 :
29 : As discussed in the documentation for the [contour](module_contour.md) module, you
30 : can generate a continuous representation for the density as a function of positions for a set
31 : of $N$ atoms with positions $(x_i,y_i,z_i)$ using:
32 :
33 : $$
34 : p(x,y,x) = \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]
35 : $$
36 :
37 : In this expression $\sigma_x, \sigma_y$ and $\sigma_z$ are bandwidth parameters and
38 : $K$ is one of a Gaussian kernel function.
39 :
40 : The Willard-Chandler surface is defined a surface of constant density in the above field $p(x,y,z)$.
41 : In other words, it is a set of points, $(x',y',z')$, in your box which have:
42 :
43 : $$
44 : p(x',y',z') = \rho
45 : $$
46 :
47 : where $\rho$ is some target density. This action calculates the distance projected on the $x, y$ or
48 : $z$ axis between the position of some test particle and this surface of constant field density.
49 :
50 : ## Examples
51 :
52 : In this example atoms 2-100 are assumed to be concentrated along some part of the $z$ axis so that you
53 : an interface between a liquid/solid and the vapor. The quantity `dc.dist1` measures the projection on $z$ of the distance
54 : between the position of atom 1 and the nearest point at which density of atoms 2-100 is equal to 0.2.
55 :
56 : ```plumed
57 : dc: DISTANCE_FROM_CONTOUR POSITIONS=2-100 ATOM=1 BANDWIDTH=0.5,0.5,0.5 DIR=z CONTOUR=0.2
58 : PRINT ARG=dc.dist1 FILE=colvar
59 : ```
60 :
61 : Notice that, as discussed in the paper cited below, if you are running with periodic boundary conditions there will be two
62 : isocontours in the box where the density is equal to 0.2. If you wish to find the distance betwene atom 1 and the second
63 : closest of these two contours you would print `dc.dist2`. `dc.thickness` tells you the difference between `dc.dist1` and
64 : `dc.dist2` and `dc.qdist` is the quantity with continuous derivatives that is introduced in the paper cited below.
65 :
66 : PLUMED also contains an experimental implementation that allows you to find the density from a isosurface in a density that is calculated using:
67 :
68 : $$
69 : p(x,y,x) = \sum_{i=1}^N w_i K\left[\frac{x-x_i}{\sigma_x},\frac{y-y_i}{\sigma_y},\frac{z-z_i}{\sigma_z} \right]
70 : $$
71 :
72 : where $w_i$ is a non-constant weight that is ascribed to each of the points. The following illustrates how this functionality can be used
73 : to find the distance from an isocontour in a phase field that describes the average value of the coordination number at each point in the box:
74 :
75 : ```plumed
76 : mat: CONTACT_MATRIX GROUPA=2-100 GROUPB=101-1000 SWITCH={RATIONAL R_0=0.2}
77 : ones: ONES SIZE=900
78 : cc: MATRIX_VECTOR_PRODUCT ARG=mat,ones
79 : dc: DISTANCE_FROM_CONTOUR ARG=cc POSITIONS=2-100 ATOM=1 BANDWIDTH=0.5,0.5,0.5 DIR=z CONTOUR=0.2
80 : PRINT ARG=dc.dist1 FILE=colvar
81 : ```
82 :
83 : Notice that, although you can calculate derivatives for the first example input above, you __cannot__ calculate derivatives if you use the ARG keyword
84 : that appears in the second example input above.
85 :
86 : */
87 : //+ENDPLUMEDOC
88 :
89 : namespace PLMD {
90 : namespace contour {
91 :
92 : class DistanceFromContour : public DistanceFromContourBase {
93 : private:
94 : unsigned dir;
95 : double pbc_param;
96 : std::vector<double> pos1, pos2, dirv, dirv2;
97 : std::vector<unsigned> perp_dirs;
98 : std::vector<Vector> atom_deriv;
99 : public:
100 : static void registerKeywords( Keywords& keys );
101 : explicit DistanceFromContour( const ActionOptions& );
102 : void calculate() override;
103 : void evaluateDerivatives( const Vector& root1, const double& root2 );
104 : };
105 :
106 : PLUMED_REGISTER_ACTION(DistanceFromContour,"DISTANCE_FROM_CONTOUR")
107 :
108 3 : void DistanceFromContour::registerKeywords( Keywords& keys ) {
109 3 : DistanceFromContourBase::registerKeywords( keys );
110 6 : keys.addOutputComponent("dist1","default","scalar","the distance between the reference atom and the nearest contour");
111 6 : keys.addOutputComponent("dist2","default","scalar","the distance between the reference atom and the other contour");
112 6 : keys.addOutputComponent("qdist","default","scalar","the differentiable (squared) distance between the two contours (see above)");
113 6 : keys.addOutputComponent("thickness","default","scalar","the distance between the two contours on the line from the reference atom");
114 3 : keys.add("compulsory","DIR","the direction perpendicular to the contour that you are looking for");
115 3 : keys.add("compulsory","TOLERANCE","0.1","this parameter is used to manage periodic boundary conditions. The problem "
116 : "here is that we can be between contours even when we are not within the membrane "
117 : "because of periodic boundary conditions. When we are in the contour, however, we "
118 : "should have it so that the sums of the absolute values of the distances to the two "
119 : "contours is approximately the distance between the two contours. There can be numerical errors in these calculations, however, so "
120 : "we specify a small tolerance here");
121 3 : keys.addDOI("10.1021/acs.jpcb.8b03661");
122 3 : }
123 :
124 1 : DistanceFromContour::DistanceFromContour( const ActionOptions& ao ):
125 : Action(ao),
126 : DistanceFromContourBase(ao),
127 2 : pos1(3,0.0),
128 1 : pos2(3,0.0),
129 1 : dirv(3,0.0),
130 1 : dirv2(3,0.0),
131 1 : perp_dirs(2),
132 2 : atom_deriv(getNumberOfAtoms()-1) {
133 : // Get the direction
134 : std::string ldir;
135 2 : parse("DIR",ldir );
136 1 : if( ldir=="x" ) {
137 0 : dir=0;
138 0 : perp_dirs[0]=1;
139 0 : perp_dirs[1]=2;
140 0 : dirv[0]=1;
141 0 : dirv2[0]=-1;
142 1 : } else if( ldir=="y" ) {
143 0 : dir=1;
144 0 : perp_dirs[0]=0;
145 0 : perp_dirs[1]=2;
146 0 : dirv[1]=1;
147 0 : dirv2[1]=-1;
148 1 : } else if( ldir=="z" ) {
149 1 : dir=2;
150 1 : perp_dirs[0]=0;
151 1 : perp_dirs[1]=1;
152 1 : dirv[2]=1;
153 1 : dirv2[2]=-1;
154 : } else {
155 0 : error(ldir + " is not a valid direction use x, y or z");
156 : }
157 :
158 : // Read in the tolerance for the pbc parameter
159 2 : parse("TOLERANCE",pbc_param);
160 :
161 : std::vector<std::size_t> shape;
162 : // Create the values
163 1 : addComponent("thickness", shape );
164 1 : componentIsNotPeriodic("thickness");
165 1 : addComponent("dist1", shape );
166 1 : componentIsNotPeriodic("dist1");
167 1 : addComponent("dist2", shape );
168 1 : componentIsNotPeriodic("dist2");
169 1 : addComponentWithDerivatives("qdist", shape );
170 2 : componentIsNotPeriodic("qdist");
171 1 : }
172 :
173 137 : void DistanceFromContour::calculate() {
174 : // Check box is orthorhombic
175 137 : if( !getPbc().isOrthorombic() ) {
176 0 : error("cell box must be orthorhombic");
177 : }
178 :
179 : // The nanoparticle is at the origin of our coordinate system
180 137 : pos1[0]=pos1[1]=pos1[2]=0.0;
181 137 : pos2[0]=pos2[1]=pos2[2]=0.0;
182 :
183 : // Set bracket as center of mass of membrane in active region
184 137 : Vector myvec = pbcDistance( getPosition(getNumberOfAtoms()-1), getPosition(0) );
185 137 : pos2[dir]=myvec[dir];
186 137 : nactive=1;
187 137 : active_list[0]=0;
188 : double d2, mindist = myvec.modulo2();
189 137 : for(unsigned j=1; j<getNumberOfAtoms()-1; ++j) {
190 0 : Vector distance=pbcDistance( getPosition(getNumberOfAtoms()-1), getPosition(j) );
191 0 : if( (d2=distance[perp_dirs[0]]*distance[perp_dirs[0]])<rcut2 &&
192 0 : (d2+=distance[perp_dirs[1]]*distance[perp_dirs[1]])<rcut2 ) {
193 0 : d2+=distance[dir]*distance[dir];
194 0 : if( d2<mindist && fabs(distance[dir])>epsilon ) {
195 0 : pos2[dir]=distance[dir];
196 : mindist = d2;
197 : }
198 0 : active_list[nactive]=j;
199 0 : nactive++;
200 : }
201 : }
202 : // pos1 position of the nanoparticle, in the first time
203 : // pos2 is the position of the closer atom in the membrane with respect the nanoparticle
204 : // fa = distance between pos1 and the contour
205 : // fb = distance between pos2 and the contour
206 137 : std::vector<double> faked(3);
207 137 : double fa = getDifferenceFromContour( pos1, faked );
208 137 : double fb = getDifferenceFromContour( pos2, faked );
209 137 : if( fa*fb>0 ) {
210 0 : unsigned maxtries = std::floor( ( getBox()(dir,dir) ) / bw[dir] );
211 0 : for(unsigned i=0; i<maxtries; ++i) {
212 0 : double sign=(pos2[dir]>0)? -1 : +1; // If the nanoparticle is inside the membrane push it out
213 0 : pos1[dir] += sign*bw[dir];
214 0 : fa = getDifferenceFromContour( pos1, faked );
215 0 : if( fa*fb<0 ) {
216 : break;
217 : }
218 : // if fa*fb is less than zero the new pos 1 is outside the contour
219 : }
220 : }
221 : // Set direction for contour search
222 137 : dirv[dir] = pos2[dir] - pos1[dir];
223 : // Bracket for second root starts in center of membrane
224 137 : double fc = getDifferenceFromContour( pos2, faked );
225 137 : if( fc*fb>0 ) {
226 : // first time is true, because fc=fb
227 : // push pos2 from its initial position inside the membrane towards the second contourn
228 137 : unsigned maxtries = std::floor( ( getBox()(dir,dir) ) / bw[dir] );
229 230 : for(unsigned i=0; i<maxtries; ++i) {
230 230 : double sign=(dirv[dir]>0)? +1 : -1;
231 230 : pos2[dir] += sign*bw[dir];
232 230 : fc = getDifferenceFromContour( pos2, faked );
233 230 : if( fc*fb<0 ) {
234 : break;
235 : }
236 : }
237 137 : dirv2[dir] = ( pos1[dir] + dirv[dir] ) - pos2[dir];
238 : }
239 :
240 : // Now do a search for the two contours
241 137 : findContour( dirv, pos1 );
242 : // Save the first value
243 : Vector root1;
244 : root1.zero();
245 137 : root1[dir] = pval[dir];
246 137 : findContour( dirv2, pos2 );
247 : // Calculate the separation between the two roots using PBC
248 : Vector root2;
249 : root2.zero();
250 137 : root2[dir] = pval[dir];
251 : Vector sep = pbcDistance( root1, root2 );
252 137 : double spacing = fabs( sep[dir] );
253 137 : plumed_assert( spacing>epsilon );
254 137 : getPntrToComponent("thickness")->set( spacing );
255 :
256 : // Make sure the sign is right
257 137 : double predir=(root1[dir]*root2[dir]<0)? -1 : 1;
258 : // This deals with periodic boundary conditions - if we are inside the membrane the sum of the absolute
259 : // distances from the contours should add up to the spacing. When this is not the case we must be outside
260 : // the contour
261 : // if( predir==-1 && (fabs(root1[dir])+fabs(root2[dir]))>(spacing+pbc_param) ) predir=1;
262 : // Set the final value to root that is closest to the "origin" = position of atom
263 137 : if( fabs(root1[dir])<fabs(root2[dir]) ) {
264 137 : getPntrToComponent("dist1")->set( predir*fabs(root1[dir]) );
265 274 : getPntrToComponent("dist2")->set( fabs(root2[dir]) );
266 : } else {
267 0 : getPntrToComponent("dist1")->set( predir*fabs(root2[dir]) );
268 0 : getPntrToComponent("dist2")->set( fabs(root1[dir]) );
269 : }
270 137 : getPntrToComponent("qdist")->set( root2[dir]*root1[dir] );
271 :
272 : // Now calculate the derivatives
273 137 : if( !doNotCalculateDerivatives() ) {
274 137 : evaluateDerivatives( root1, root2[dir] );
275 137 : evaluateDerivatives( root2, root1[dir] );
276 : }
277 137 : }
278 :
279 274 : void DistanceFromContour::evaluateDerivatives( const Vector& root1, const double& root2 ) {
280 274 : if( getNumberOfArguments()>0 ) {
281 0 : plumed_merror("derivatives for phase field distance from contour have not been implemented yet");
282 : }
283 :
284 : Vector origind;
285 : origind.zero();
286 : Tensor vir;
287 : vir.zero();
288 : double sumd = 0;
289 274 : std::vector<double> pp(3), ddd(3,0);
290 548 : for(unsigned i=0; i<nactive; ++i) {
291 274 : /* double newval = */evaluateKernel( getPosition(active_list[i]), root1, ddd );
292 274 : Vector distance = pbcDistance( getPosition(getNumberOfAtoms()-1), getPosition(active_list[i]) );
293 :
294 274 : if( getNumberOfArguments()==1 ) {
295 : } else {
296 274 : sumd += ddd[dir];
297 1096 : for(unsigned j=0; j<3; ++j) {
298 822 : atom_deriv[i][j] = -ddd[j];
299 : }
300 : origind += -atom_deriv[i];
301 548 : vir -= Tensor(atom_deriv[i],distance);
302 : }
303 : }
304 :
305 : // Add derivatives to atoms involved
306 274 : Value* val=getPntrToComponent("qdist");
307 274 : double prefactor = root2 / sumd;
308 548 : for(unsigned i=0; i<nactive; ++i) {
309 274 : val->addDerivative( 3*active_list[i] + 0, -prefactor*atom_deriv[i][0] );
310 274 : val->addDerivative( 3*active_list[i] + 1, -prefactor*atom_deriv[i][1] );
311 274 : val->addDerivative( 3*active_list[i] + 2, -prefactor*atom_deriv[i][2] );
312 : }
313 :
314 : // Add derivatives to atoms at origin
315 274 : unsigned nbase = 3*(getNumberOfAtoms()-1);
316 274 : val->addDerivative( nbase, -prefactor*origind[0] );
317 : nbase++;
318 274 : val->addDerivative( nbase, -prefactor*origind[1] );
319 : nbase++;
320 274 : val->addDerivative( nbase, -prefactor*origind[2] );
321 : nbase++;
322 :
323 : // Add derivatives to virial
324 1096 : for(unsigned i=0; i<3; ++i)
325 3288 : for(unsigned j=0; j<3; ++j) {
326 2466 : val->addDerivative( nbase, -prefactor*vir(i,j) );
327 : nbase++;
328 : }
329 274 : }
330 :
331 : }
332 : }
|
