Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 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 caculates 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 concentraed along some part of the \f$z\f$ axis so that you
58 : an interface between a liquid/solid and the vapour. 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 : \verbatim
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 : \endverbatim
65 :
66 : */
67 : //+ENDPLUMEDOC
68 :
69 : namespace PLMD {
70 : namespace multicolvar {
71 :
72 2 : class DistanceFromContour : public MultiColvarBase {
73 : private:
74 : unsigned dir;
75 : bool derivTime;
76 : double rcut2;
77 : double contour;
78 : std::string kerneltype;
79 : std::vector<Value*> pval;
80 : std::vector<double> bw, pos, dirv;
81 : std::vector<double> forces;
82 : std::vector<unsigned> perp_dirs;
83 : vesselbase::ValueVessel* myvalue_vessel;
84 : vesselbase::ValueVessel* myderiv_vessel;
85 : RootFindingBase<DistanceFromContour> mymin;
86 : public:
87 : static void registerKeywords( Keywords& keys );
88 : explicit DistanceFromContour( const ActionOptions& );
89 1293 : bool isDensity() const { return true; }
90 : void calculate();
91 : unsigned getNumberOfQuantities() const ;
92 0 : bool isPeriodic() { return false; }
93 : double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
94 : double getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der );
95 : // We need an apply action as we are using an independent value
96 : void apply();
97 : };
98 :
99 2524 : PLUMED_REGISTER_ACTION(DistanceFromContour,"DISTANCE_FROM_CONTOUR")
100 :
101 2 : void DistanceFromContour::registerKeywords( Keywords& keys ) {
102 2 : MultiColvarBase::registerKeywords( keys );
103 2 : keys.add("compulsory","DATA","The input base multicolvar which is being used to calculate the contour");
104 2 : keys.add("atoms","ATOM","The atom whose perpendicular distance we are calculating from the contour");
105 2 : keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density esimtation");
106 : keys.add("compulsory","KERNEL","gaussian","the kernel function you are using. More details on the kernels available "
107 2 : "in plumed plumed can be found in \\ref kernelfunctions.");
108 2 : keys.add("compulsory","DIR","the direction perpendicular to the contour that you are looking for");
109 2 : keys.add("compulsory","CONTOUR","the value we would like for the contour");
110 2 : }
111 :
112 1 : DistanceFromContour::DistanceFromContour( const ActionOptions& ao ):
113 : Action(ao),
114 : MultiColvarBase(ao),
115 : derivTime(false),
116 : bw(3),
117 : pos(3,0.0),
118 : dirv(3,0.0),
119 : perp_dirs(2),
120 1 : mymin(this)
121 : {
122 : // Read in the multicolvar/atoms
123 1 : std::vector<AtomNumber> all_atoms;
124 1 : bool read2 = parseMultiColvarAtomList("DATA", -1, all_atoms);
125 1 : if( !read2 ) error("missing DATA keyword");
126 1 : bool read1 = parseMultiColvarAtomList("ATOM", -1, all_atoms);
127 1 : if( !read1 ) error("missing ATOM keyword");
128 1 : if( all_atoms.size()!=1 ) error("should only be one atom specified");
129 : // Read in the center of the binding object
130 1 : log.printf(" computing distance of atom %d from contour \n",all_atoms[0].serial() );
131 1 : setupMultiColvarBase( all_atoms ); forces.resize( 3*getNumberOfAtoms() + 9 );
132 1 : if( getNumberOfBaseMultiColvars()!=1 ) error("should only be one input multicolvar");
133 :
134 : // Get the direction
135 2 : std::string ldir; parse("DIR",ldir );
136 1 : if( ldir=="x" ) { dir=0; perp_dirs[0]=1; perp_dirs[1]=2; dirv[0]=1; }
137 1 : else if( ldir=="y" ) { dir=1; perp_dirs[0]=0; perp_dirs[1]=2; dirv[1]=1; }
138 1 : else if( ldir=="z" ) { dir=2; perp_dirs[0]=0; perp_dirs[1]=1; dirv[2]=1; }
139 0 : else error(ldir + " is not a valid direction use x, y or z");
140 :
141 : // Read in details of phase field construction
142 1 : parseVector("BANDWIDTH",bw); parse("KERNEL",kerneltype); parse("CONTOUR",contour);
143 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() );
144 1 : log.printf(" constructing phase field using %s kernels with bandwidth (%f, %f, %f) \n",kerneltype.c_str(), bw[0], bw[1], bw[2] );
145 :
146 : // Now create a task list
147 1 : for(unsigned i=0; i<mybasemulticolvars[0]->getFullNumberOfTasks(); ++i) addTaskToList(i);
148 : // And a cutoff
149 2 : std::vector<double> pp( bw.size(),0 );
150 2 : KernelFunctions kernel( pp, bw, kerneltype, false, 1.0, true );
151 1 : double rcut = kernel.getCutoff( bw[0] );
152 3 : for(unsigned j=1; j<bw.size(); ++j) {
153 2 : if( kernel.getCutoff(bw[j])>rcut ) rcut=kernel.getCutoff(bw[j]);
154 : }
155 1 : rcut2=rcut*rcut;
156 : // Create the value
157 1 : addValueWithDerivatives(); setNotPeriodic();
158 : // Create sum vessels
159 2 : std::string fake_input; std::string deriv_input="COMPONENT=2";
160 1 : if( mybasemulticolvars[0]->isDensity() ) {
161 1 : addVessel( "SUM", fake_input, -1 ); addVessel( "SUM", deriv_input, -1 );
162 : } else {
163 0 : addVessel( "MEAN", fake_input, -1 ); addVessel( "MEAN", deriv_input, -1 );
164 : }
165 : // And convert to a value vessel so we can get the final value
166 1 : myvalue_vessel = dynamic_cast<vesselbase::ValueVessel*>( getPntrToVessel(0) );
167 1 : myderiv_vessel = dynamic_cast<vesselbase::ValueVessel*>( getPntrToVessel(1) );
168 1 : plumed_assert( myvalue_vessel && myderiv_vessel ); resizeFunctions();
169 :
170 : // Create the vector of values that holds the position
171 2 : for(unsigned i=0; i<3; ++i) pval.push_back( new Value() );
172 1 : }
173 :
174 2586 : unsigned DistanceFromContour::getNumberOfQuantities() const {
175 2586 : return 3;
176 : }
177 :
178 137 : void DistanceFromContour::calculate() {
179 : // Check box is orthorhombic
180 137 : if( !getPbc().isOrthorombic() ) error("cell box must be orthorhombic");
181 :
182 : // The nanoparticle is at the origin of our coordinate system
183 137 : pos[0]=pos[1]=pos[2]=0.0;
184 :
185 : // Set bracket as center of mass of membrane in active region
186 137 : double d2, norm=0.0; dirv[dir]=0; deactivateAllTasks();
187 274 : for(unsigned j=0; j<getNumberOfAtoms()-1; ++j) {
188 137 : Vector distance=getSeparation( getPosition(getNumberOfAtoms()-1), getPosition(j) );
189 411 : if( (d2=distance[perp_dirs[0]]*distance[perp_dirs[0]])<rcut2 &&
190 274 : (d2+=distance[perp_dirs[1]]*distance[perp_dirs[1]])<rcut2 ) {
191 137 : dirv[dir]+=distance[dir]; taskFlags[j]=1; norm+=1.0;
192 : }
193 : }
194 137 : dirv[dir] = dirv[dir] / norm;
195 137 : lockContributors(); derivTime=false;
196 :
197 : // Now do a search for the contour
198 137 : mymin.lsearch( dirv, pos, &DistanceFromContour::getDifferenceFromContour );
199 :
200 : // Set the final value
201 137 : setValue( pval[dir]->get() );
202 : // Now calculate the derivatives
203 137 : if( !doNotCalculateDerivatives() ) {
204 137 : Value* ival=myvalue_vessel->getFinalValue(); ival->clearDerivatives();
205 137 : derivTime=true; double prefactor; std::vector<double> der(3); getDifferenceFromContour( pos, der );
206 137 : if( mybasemulticolvars[0]->isDensity() ) prefactor = 1.0 / myderiv_vessel->getOutputValue(); else plumed_error();
207 137 : Value* val=getPntrToValue();
208 137 : for(unsigned i=0; i<val->getNumberOfDerivatives(); ++i) val->setDerivative( i, -prefactor*ival->getDerivative(i) );
209 : }
210 137 : }
211 :
212 1293 : double DistanceFromContour::getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der ) {
213 2586 : std::string min, max;
214 5172 : for(unsigned j=0; j<3; ++j) {
215 3879 : Tools::convert( -0.5*getBox()(j,j), min );
216 3879 : Tools::convert( +0.5*getBox()(j,j), max );
217 3879 : pval[j]->setDomain( min, max ); pval[j]->set( x[j] );
218 : }
219 1293 : runAllTasks();
220 2586 : return myvalue_vessel->getOutputValue() - contour;
221 : }
222 :
223 1293 : double DistanceFromContour::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
224 1293 : Vector distance = getSeparation( getPosition(getNumberOfAtoms()-1), myatoms.getPosition(0) );
225 2586 : std::vector<double> pp(3), der(3,0); for(unsigned j=0; j<3; ++j) pp[j] = distance[j];
226 :
227 : // Now create the kernel and evaluate
228 2586 : KernelFunctions kernel( pp, bw, kerneltype, false, 1.0, true );
229 1293 : double newval = kernel.evaluate( pval, der, true );
230 :
231 1293 : if( mybasemulticolvars[0]->isDensity() ) {
232 1293 : if( !doNotCalculateDerivatives() && derivTime ) {
233 137 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
234 137 : Vector vder; unsigned basen=myvals.getNumberOfDerivatives() - 12;
235 548 : for(unsigned i=0; i<3; ++i) {
236 411 : vder[i]=der[i]; myvals.addDerivative( 1, basen+i, vder[i] );
237 : }
238 137 : myatoms.setValue( 2, der[dir] );
239 137 : addAtomDerivatives( 1, 0, -vder, myatoms );
240 137 : myatoms.addBoxDerivatives( 1, Tensor(vder,distance) );
241 : }
242 1293 : myatoms.setValue( 0, 1.0 ); return newval;
243 : }
244 :
245 : // This does the stuff for averaging
246 0 : myatoms.setValue( 0, newval );
247 :
248 : // This gets the average if we are using a phase field
249 0 : std::vector<double> cvals( mybasemulticolvars[0]->getNumberOfQuantities() );
250 0 : mybasedata[0]->retrieveValueWithIndex( tindex, false, cvals );
251 1293 : return newval*cvals[0]*cvals[1];
252 : }
253 :
254 137 : void DistanceFromContour::apply() {
255 137 : if( getPntrToValue()->applyForce( forces ) ) setForcesOnAtoms( forces );
256 137 : }
257 :
258 : }
259 2523 : }
|
