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 "MultiColvarBase.h"
23 : #include "AtomValuePack.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/KernelFunctions.h"
26 : #include "tools/RootFindingBase.h"
27 : #include "vesselbase/ValueVessel.h"
28 :
29 : //+PLUMEDOC COLVAR DISTANCE_FROM_CONTOUR
30 : /*
31 : Calculate the perpendicular distance from a Willard-Chandler dividing surface.
32 :
33 : Suppose that you have calculated a multicolvar. By doing so you have calculated a
34 : set of colvars, \f$s_i\f$, and each of these colvars has a well defined position in
35 : space \f$(x_i,y_i,z_i)\f$. You can use this information to calculate a phase-field
36 : model of the colvar density using:
37 :
38 : \f[
39 : 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]
40 : \f]
41 :
42 : In this expression \f$\sigma_x, \sigma_y\f$ and \f$\sigma_z\f$ are bandwidth parameters and
43 : \f$K\f$ is one of the \ref kernelfunctions. This is what is done within \ref MULTICOLVARDENS
44 :
45 : The Willard-Chandler surface is a surface of constant density in the above phase field \f$p(x,y,z)\f$.
46 : In other words, it is a set of points, \f$(x',y',z')\f$, in your box which have:
47 :
48 : \f[
49 : p(x',y',z') = \rho
50 : \f]
51 :
52 : where \f$\rho\f$ is some target density. This action calculates the distance projected on the \f$x, y\f$ or
53 : \f$z\f$ axis between the position of some test particle and this surface of constant field density.
54 :
55 : \par Examples
56 :
57 : In this example atoms 2-100 are assumed to be concentrated along some part of the \f$z\f$ axis so that you
58 : an interface between a liquid/solid and the vapor. The quantity dc measures the distance between the
59 : surface at which the density of 2-100 atoms is equal to 0.2 and the position of the test particle atom 1.
60 :
61 : \plumedfile
62 : dens: DENSITY SPECIES=2-100
63 : dc: DISTANCE_FROM_CONTOUR DATA=dens ATOM=1 BANDWIDTH=0.5,0.5,0.5 DIR=z CONTOUR=0.2
64 : \endplumedfile
65 :
66 : */
67 : //+ENDPLUMEDOC
68 :
69 : namespace PLMD {
70 : namespace multicolvar {
71 :
72 : class DistanceFromContour : public MultiColvarBase {
73 : private:
74 : unsigned dir;
75 : bool derivTime;
76 : double rcut2;
77 : double contour;
78 : double pbc_param;
79 : std::string kerneltype;
80 : std::vector<std::unique_ptr<Value>> pval;
81 : std::vector<double> bw, pos1, pos2, dirv, dirv2;
82 : std::vector<double> forces;
83 : std::vector<unsigned> perp_dirs;
84 : vesselbase::ValueVessel* myvalue_vessel;
85 : vesselbase::ValueVessel* myderiv_vessel;
86 : RootFindingBase<DistanceFromContour> mymin;
87 : public:
88 : static void registerKeywords( Keywords& keys );
89 : explicit DistanceFromContour( const ActionOptions& );
90 3115 : bool isDensity() const override {
91 3115 : return true;
92 : }
93 : void calculate() override;
94 : unsigned getNumberOfQuantities() const override;
95 0 : bool isPeriodic() override {
96 0 : return false;
97 : }
98 : double compute( const unsigned& tindex, AtomValuePack& myatoms ) const override;
99 : double getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der );
100 : // We need an apply action as we are using an independent value
101 : void apply() override;
102 : };
103 :
104 13787 : PLUMED_REGISTER_ACTION(DistanceFromContour,"DISTANCE_FROM_CONTOUR")
105 :
106 5 : void DistanceFromContour::registerKeywords( Keywords& keys ) {
107 5 : MultiColvarBase::registerKeywords( keys );
108 10 : keys.addOutputComponent("dist1","default","the distance between the reference atom and the nearest contour");
109 10 : keys.addOutputComponent("dist2","default","the distance between the reference atom and the other contour");
110 10 : keys.addOutputComponent("qdist","default","the differentiable (squared) distance between the two contours (see above)");
111 10 : keys.addOutputComponent("thickness","default","the distance between the two contours on the line from the reference atom");
112 10 : keys.add("compulsory","DATA","The input base multicolvar which is being used to calculate the contour");
113 10 : keys.add("atoms","ATOM","The atom whose perpendicular distance we are calculating from the contour");
114 10 : keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density estimation");
115 10 : keys.add("compulsory","KERNEL","gaussian","the kernel function you are using. More details on the kernels available "
116 : "in plumed plumed can be found in \\ref kernelfunctions.");
117 10 : keys.add("compulsory","DIR","the direction perpendicular to the contour that you are looking for");
118 10 : keys.add("compulsory","CONTOUR","the value we would like for the contour");
119 10 : keys.add("compulsory","TOLERANCE","0.1","this parameter is used to manage periodic boundary conditions. The problem "
120 : "here is that we can be between contours even when we are not within the membrane "
121 : "because of periodic boundary conditions. When we are in the contour, however, we "
122 : "should have it so that the sums of the absolute values of the distances to the two "
123 : "contours is approximately the distance between the two contours. There can be numerical errors in these calculations, however, so "
124 : "we specify a small tolerance here");
125 5 : }
126 :
127 1 : DistanceFromContour::DistanceFromContour( const ActionOptions& ao ):
128 : Action(ao),
129 : MultiColvarBase(ao),
130 1 : derivTime(false),
131 1 : bw(3),
132 1 : pos1(3,0.0),
133 1 : pos2(3,0.0),
134 1 : dirv(3,0.0),
135 1 : dirv2(3,0.0),
136 1 : perp_dirs(2),
137 2 : mymin(this) {
138 : // Read in the multicolvar/atoms
139 : std::vector<AtomNumber> all_atoms;
140 1 : parse("TOLERANCE",pbc_param);
141 1 : bool read2 = parseMultiColvarAtomList("DATA", -1, all_atoms);
142 1 : if( !read2 ) {
143 0 : error("missing DATA keyword");
144 : }
145 1 : bool read1 = parseMultiColvarAtomList("ATOM", -1, all_atoms);
146 1 : if( !read1 ) {
147 0 : error("missing ATOM keyword");
148 : }
149 1 : if( all_atoms.size()!=1 ) {
150 0 : error("should only be one atom specified");
151 : }
152 : // Read in the center of the binding object
153 1 : log.printf(" computing distance of atom %d from contour \n",all_atoms[0].serial() );
154 1 : setupMultiColvarBase( all_atoms );
155 1 : forces.resize( 3*getNumberOfAtoms() + 9 );
156 1 : if( getNumberOfBaseMultiColvars()!=1 ) {
157 0 : error("should only be one input multicolvar");
158 : }
159 :
160 : // Get the direction
161 : std::string ldir;
162 2 : parse("DIR",ldir );
163 1 : if( ldir=="x" ) {
164 0 : dir=0;
165 0 : perp_dirs[0]=1;
166 0 : perp_dirs[1]=2;
167 0 : dirv[0]=1;
168 0 : dirv2[0]=-1;
169 1 : } else if( ldir=="y" ) {
170 0 : dir=1;
171 0 : perp_dirs[0]=0;
172 0 : perp_dirs[1]=2;
173 0 : dirv[1]=1;
174 0 : dirv2[1]=-1;
175 1 : } else if( ldir=="z" ) {
176 1 : dir=2;
177 1 : perp_dirs[0]=0;
178 1 : perp_dirs[1]=1;
179 1 : dirv[2]=1;
180 1 : dirv2[2]=-1;
181 : } else {
182 0 : error(ldir + " is not a valid direction use x, y or z");
183 : }
184 :
185 : // Read in details of phase field construction
186 1 : parseVector("BANDWIDTH",bw);
187 1 : parse("KERNEL",kerneltype);
188 1 : parse("CONTOUR",contour);
189 1 : log.printf(" searching for contour in %s direction at %f in phase field for multicolvar %s \n",ldir.c_str(), contour, mybasemulticolvars[0]->getLabel().c_str() );
190 1 : log.printf(" constructing phase field using %s kernels with bandwidth (%f, %f, %f) \n",kerneltype.c_str(), bw[0], bw[1], bw[2] );
191 :
192 : // Now create a task list
193 2 : for(unsigned i=0; i<mybasemulticolvars[0]->getFullNumberOfTasks(); ++i) {
194 1 : addTaskToList(i);
195 : }
196 : // And a cutoff
197 1 : std::vector<double> pp( bw.size(),0 );
198 2 : KernelFunctions kernel( pp, bw, kerneltype, "DIAGONAL", 1.0 );
199 1 : double rcut = kernel.getCutoff( bw[0] );
200 3 : for(unsigned j=1; j<bw.size(); ++j) {
201 2 : if( kernel.getCutoff(bw[j])>rcut ) {
202 0 : rcut=kernel.getCutoff(bw[j]);
203 : }
204 : }
205 1 : rcut2=rcut*rcut;
206 : // Create the values
207 1 : addComponent("thickness");
208 1 : componentIsNotPeriodic("thickness");
209 1 : addComponent("dist1");
210 1 : componentIsNotPeriodic("dist1");
211 1 : addComponent("dist2");
212 1 : componentIsNotPeriodic("dist2");
213 1 : addComponentWithDerivatives("qdist");
214 2 : componentIsNotPeriodic("qdist");
215 : // Create sum vessels
216 : std::string fake_input;
217 1 : std::string deriv_input="COMPONENT=2";
218 1 : if( mybasemulticolvars[0]->isDensity() ) {
219 1 : addVessel( "SUM", fake_input, -1 );
220 2 : addVessel( "SUM", deriv_input, -1 );
221 : } else {
222 0 : addVessel( "MEAN", fake_input, -1 );
223 0 : addVessel( "MEAN", deriv_input, -1 );
224 : }
225 : // And convert to a value vessel so we can get the final value
226 1 : myvalue_vessel = dynamic_cast<vesselbase::ValueVessel*>( getPntrToVessel(0) );
227 1 : myderiv_vessel = dynamic_cast<vesselbase::ValueVessel*>( getPntrToVessel(1) );
228 1 : plumed_assert( myvalue_vessel && myderiv_vessel );
229 1 : resizeFunctions();
230 :
231 : // Create the vector of values that holds the position
232 4 : for(unsigned i=0; i<3; ++i) {
233 3 : pval.emplace_back( Tools::make_unique<Value>() );
234 : }
235 1 : }
236 :
237 6230 : unsigned DistanceFromContour::getNumberOfQuantities() const {
238 6230 : return 3;
239 : }
240 :
241 137 : void DistanceFromContour::calculate() {
242 : // Check box is orthorhombic
243 137 : if( !getPbc().isOrthorombic() ) {
244 0 : error("cell box must be orthorhombic");
245 : }
246 :
247 : // The nanoparticle is at the origin of our coordinate system
248 137 : pos1[0]=pos1[1]=pos1[2]=0.0;
249 137 : pos2[0]=pos2[1]=pos2[2]=0.0;
250 :
251 : // Set bracket as center of mass of membrane in active region
252 137 : deactivateAllTasks();
253 137 : Vector myvec = getSeparation( getPosition(getNumberOfAtoms()-1), getPosition(0) );
254 137 : pos2[dir]=myvec[dir];
255 137 : taskFlags[0]=1;
256 137 : double mindist = myvec.modulo2();
257 137 : for(unsigned j=1; j<getNumberOfAtoms()-1; ++j) {
258 0 : Vector distance=getSeparation( getPosition(getNumberOfAtoms()-1), getPosition(j) );
259 : double d2;
260 0 : if( (d2=distance[perp_dirs[0]]*distance[perp_dirs[0]])<rcut2 &&
261 0 : (d2+=distance[perp_dirs[1]]*distance[perp_dirs[1]])<rcut2 ) {
262 0 : d2+=distance[dir]*distance[dir];
263 0 : if( d2<mindist && std::fabs(distance[dir])>epsilon ) {
264 0 : pos2[dir]=distance[dir];
265 : mindist = d2;
266 : }
267 0 : taskFlags[j]=1;
268 : }
269 : }
270 137 : lockContributors();
271 137 : derivTime=false;
272 : // pos1 position of the nanoparticle, in the first time
273 : // pos2 is the position of the closer atom in the membrane with respect the nanoparticle
274 : // fa = distance between pos1 and the contour
275 : // fb = distance between pos2 and the contour
276 137 : std::vector<double> faked(3);
277 137 : double fa = getDifferenceFromContour( pos1, faked );
278 137 : double fb = getDifferenceFromContour( pos2, faked );
279 137 : if( fa*fb>0 ) {
280 0 : unsigned maxtries = std::floor( ( getBox()(dir,dir) ) / bw[dir] );
281 0 : for(unsigned i=0; i<maxtries; ++i) {
282 0 : double sign=(pos2[dir]>0)? -1 : +1; // If the nanoparticle is inside the membrane push it out
283 0 : pos1[dir] += sign*bw[dir];
284 0 : fa = getDifferenceFromContour( pos1, faked );
285 0 : if( fa*fb<0 ) {
286 : break;
287 : }
288 : // if fa*fb is less than zero the new pos 1 is outside the contour
289 : }
290 : }
291 : // Set direction for contour search
292 137 : dirv[dir] = pos2[dir] - pos1[dir];
293 : // Bracket for second root starts in center of membrane
294 137 : double fc = getDifferenceFromContour( pos2, faked );
295 137 : if( fc*fb>0 ) {
296 : // first time is true, because fc=fb
297 : // push pos2 from its initial position inside the membrane towards the second contourn
298 137 : unsigned maxtries = std::floor( ( getBox()(dir,dir) ) / bw[dir] );
299 230 : for(unsigned i=0; i<maxtries; ++i) {
300 230 : double sign=(dirv[dir]>0)? +1 : -1;
301 230 : pos2[dir] += sign*bw[dir];
302 230 : fc = getDifferenceFromContour( pos2, faked );
303 230 : if( fc*fb<0 ) {
304 : break;
305 : }
306 : }
307 137 : dirv2[dir] = ( pos1[dir] + dirv[dir] ) - pos2[dir];
308 : }
309 :
310 : // Now do a search for the two contours
311 137 : mymin.lsearch( dirv, pos1, &DistanceFromContour::getDifferenceFromContour );
312 : // Save the first value
313 137 : Vector root1;
314 137 : root1.zero();
315 137 : root1[dir] = pval[dir]->get();
316 137 : mymin.lsearch( dirv2, pos2, &DistanceFromContour::getDifferenceFromContour );
317 : // Calculate the separation between the two roots using PBC
318 137 : Vector root2;
319 137 : root2.zero();
320 137 : root2[dir]=pval[dir]->get();
321 137 : Vector sep = getSeparation( root1, root2 );
322 137 : double spacing = std::fabs( sep[dir] );
323 137 : plumed_assert( spacing>epsilon );
324 137 : getPntrToComponent("thickness")->set( spacing );
325 :
326 : // Make sure the sign is right
327 137 : double predir=(root1[dir]*root2[dir]<0)? -1 : 1;
328 : // This deals with periodic boundary conditions - if we are inside the membrane the sum of the absolute
329 : // distances from the contours should add up to the spacing. When this is not the case we must be outside
330 : // the contour
331 : // if( predir==-1 && (std::fabs(root1[dir])+std::fabs(root2[dir]))>(spacing+pbc_param) ) predir=1;
332 : // Set the final value to root that is closest to the "origin" = position of atom
333 137 : if( std::fabs(root1[dir])<std::fabs(root2[dir]) ) {
334 137 : getPntrToComponent("dist1")->set( predir*std::fabs(root1[dir]) );
335 274 : getPntrToComponent("dist2")->set( std::fabs(root2[dir]) );
336 : } else {
337 0 : getPntrToComponent("dist1")->set( predir*std::fabs(root2[dir]) );
338 0 : getPntrToComponent("dist2")->set( std::fabs(root1[dir]) );
339 : }
340 137 : getPntrToComponent("qdist")->set( root2[dir]*root1[dir] );
341 :
342 : // Now calculate the derivatives
343 137 : if( !doNotCalculateDerivatives() ) {
344 137 : Value* ival=myvalue_vessel->getFinalValue();
345 : ival->clearDerivatives();
346 137 : std::vector<double> root1v(3);
347 548 : for(unsigned i=0; i<3; ++i) {
348 411 : root1v[i]=root1[i];
349 : }
350 137 : derivTime=true;
351 137 : std::vector<double> der(3);
352 137 : getDifferenceFromContour( root1v, der );
353 : double prefactor;
354 137 : if( mybasemulticolvars[0]->isDensity() ) {
355 137 : prefactor = root2[dir] / myderiv_vessel->getOutputValue();
356 : } else {
357 0 : plumed_error();
358 : }
359 137 : Value* val=getPntrToComponent("qdist");
360 2192 : for(unsigned i=0; i<val->getNumberOfDerivatives(); ++i) {
361 2055 : val->setDerivative( i, -prefactor*ival->getDerivative(i) );
362 : }
363 : ival->clearDerivatives();
364 137 : std::vector<double> root2v(3);
365 548 : for(unsigned i=0; i<3; ++i) {
366 411 : root2v[i]=root2[i];
367 : }
368 137 : getDifferenceFromContour( root2v, der );
369 137 : if( mybasemulticolvars[0]->isDensity() ) {
370 137 : prefactor = root1[dir] / myderiv_vessel->getOutputValue();
371 : } else {
372 0 : plumed_error();
373 : }
374 2192 : for(unsigned i=0; i<val->getNumberOfDerivatives(); ++i) {
375 2055 : val->addDerivative( i, -prefactor*ival->getDerivative(i) );
376 : }
377 : }
378 137 : }
379 :
380 3115 : double DistanceFromContour::getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der ) {
381 : std::string min, max;
382 12460 : for(unsigned j=0; j<3; ++j) {
383 9345 : Tools::convert( -0.5*getBox()(j,j), min );
384 9345 : Tools::convert( +0.5*getBox()(j,j), max );
385 9345 : pval[j]->setDomain( min, max );
386 9345 : pval[j]->set( x[j] );
387 : }
388 3115 : runAllTasks();
389 6230 : return myvalue_vessel->getOutputValue() - contour;
390 : }
391 :
392 3115 : double DistanceFromContour::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
393 3115 : Vector distance = getSeparation( getPosition(getNumberOfAtoms()-1), myatoms.getPosition(0) );
394 3115 : std::vector<double> pp(3), der(3,0);
395 12460 : for(unsigned j=0; j<3; ++j) {
396 9345 : pp[j] = distance[j];
397 : }
398 :
399 : // convert the pointer once
400 3115 : auto pval_ptr=Tools::unique2raw(pval);
401 :
402 : // Now create the kernel and evaluate
403 3115 : KernelFunctions kernel( pp, bw, kerneltype, "DIAGONAL", 1.0 );
404 3115 : kernel.normalize( pval_ptr );
405 3115 : double newval = kernel.evaluate( pval_ptr, der, true );
406 :
407 3115 : if( mybasemulticolvars[0]->isDensity() ) {
408 3115 : if( !doNotCalculateDerivatives() && derivTime ) {
409 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
410 274 : Vector vder;
411 274 : unsigned basen=myvals.getNumberOfDerivatives() - 12;
412 1096 : for(unsigned i=0; i<3; ++i) {
413 822 : vder[i]=der[i];
414 822 : myvals.addDerivative( 1, basen+i, vder[i] );
415 : }
416 274 : myatoms.setValue( 2, der[dir] );
417 274 : addAtomDerivatives( 1, 0, -vder, myatoms );
418 274 : myatoms.addBoxDerivatives( 1, Tensor(vder,distance) );
419 : }
420 : myatoms.setValue( 0, 1.0 );
421 3115 : return newval;
422 : }
423 :
424 : // This does the stuff for averaging
425 : myatoms.setValue( 0, newval );
426 :
427 : // This gets the average if we are using a phase field
428 0 : std::vector<double> cvals( mybasemulticolvars[0]->getNumberOfQuantities() );
429 0 : mybasedata[0]->retrieveValueWithIndex( tindex, false, cvals );
430 0 : return newval*cvals[0]*cvals[1];
431 : }
432 :
433 137 : void DistanceFromContour::apply() {
434 274 : if( getPntrToComponent("qdist")->applyForce( forces ) ) {
435 0 : setForcesOnAtoms( forces );
436 : }
437 137 : }
438 :
439 : }
440 : }
|
