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