1 /*
2 This file is part of MADNESS.
3 Copyright (C) 2007,2010 Oak Ridge National Laboratory
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2 of the License, or
7 (at your option) any later version.
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12 You should have received a copy of the GNU General Public License
13 along with this program; if not, write to the Free Software
14 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
15
16 For more information please contact:
17 Robert J. Harrison
18 Oak Ridge National Laboratory
19 One Bethel Valley Road
20 P.O. Box 2008, MS-6367
21
22 email: harrisonrj@ornl.gov
23 tel: 865-241-3937
24 fax: 865-572-0680
25
26
27 This proram simulates the effect of surface solute interaction between a colloid and a hydrogen atom
28 */
29
30
31 //#define WORLD_INSTANTIATE_STATIC_TEMPLATES
32 //#include <madness/mra/operator.h>
33 #include "molecularmask.h"
34 #include <madness/mra/nonlinsol.h>
35 #include <madness/mra/mra.h>
36 #include <madness/mra/lbdeux.h>
37 #include <madness/misc/ran.h>
38 #include <madness/tensor/solvers.h>
39 #include <ctime>
40 #include <list>
41 #include <jacob/molecule.h>
42 #include <madness/mra/sdf_shape_3D.h>
43 #include <madness/mra/funcplot.h>
44 #include <madness/constants.h>
45 #include <cmath>
46 #include <vector>
47 using namespace madness;
48 using namespace std;
49 typedef real_function_3d realfunc;
50 //typedef SeparatedConvolution<double,3> operatorT;
51 //typedef std::shared_ptr<operatorT> poperatorT;
52
53 //compute the distance between two points
distance1(const coord_3d & r,const coord_3d & center)54 inline double distance1(const coord_3d& r, const coord_3d& center){
55 double x1 = r[0], y1 = r[1], z1 = r[2];
56 double x2 = center[0], y2 = center[1], z2 = center[2];
57 double xx = x1-x2;
58 double yy = y1-y2;
59 double zz = z1-z2;
60 return sqrt(xx*xx + yy*yy + zz*zz);
61 }
62
nuclear_charge_function(const coord_3d & r)63 double nuclear_charge_function(const coord_3d& r) {
64 const double expnt = 100.0;
65 const double coeff = pow(1.0/constants::pi*expnt,0.5*3);
66 return coeff*exp(-expnt*(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]));
67 }
68
electronic_charge_function(const coord_3d & r)69 double electronic_charge_function(const coord_3d& r) {
70 const double coeff = 1.0/constants::pi;
71 return coeff*exp(-2.0*sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]));
72 }
73
charge_function(const coord_3d & r)74 double charge_function(const coord_3d& r) {
75 return nuclear_charge_function(r) - electronic_charge_function(r);
76 }
77
78
79 class SurfaceMoleculeInteraction {
80 private:
81 const double& d; //distance between reaction surface and molecular center
82 const double& R; //radius of the spheres forming the colloid surface
83 const std::vector<double> charge_center; //defines the center of nuclear charge
84 realfunc& rhot; //total charge density
85 const double& sigma; //surface width
86 World& world;
87 const int& maxiter;
88 //compute the reciprocal of a madness function
89 template <typename T, int NDIM>
90 struct Reciprocal {
operator ()SurfaceMoleculeInteraction::Reciprocal91 void operator()(const Key<NDIM>& key, Tensor<T>& t) const {
92 UNARY_OPTIMIZED_ITERATOR(T, t, *_p0 = 1.0/(*_p0));
93 }
serializeSurfaceMoleculeInteraction::Reciprocal94 template <typename Archive> void serialize(Archive& ar) {}
95 };
96
97 //this is used to perform a binary operation
98 struct Bop {
operator ()SurfaceMoleculeInteraction::Bop99 void operator()(const Key<3>& key,
100 real_tensor rfunc,
101 const real_tensor& func1,
102 const real_tensor& func2) const {
103 ITERATOR(rfunc,
104 double d = func1(IND);
105 double p = func2(IND);
106 rfunc(IND) = d*p;
107 );
108 }
109
110 template <typename Archive>
serializeSurfaceMoleculeInteraction::Bop111 void serialize(Archive& ar) {}
112 };
113 public:
114 //set coordinates of colloid atoms
colloid_coords() const115 std::vector< madness::Vector<double,3> > colloid_coords()const{
116 // const std::vector<double> cc(0.0);// = charge_center;
117 std::vector< madness::Vector<double,3> > c(6); //6 spheres on the colloid surface
118 double sqrttwo = std::sqrt(2.0);
119 double dist= (sqrttwo/2.0)*R;
120 //double x = cc[0], y = cc[1], z = cc[2];
121 double x = 0.0, y = 0.0, z = 0.0;
122 c[0][0]= x - dist, c[0][1]= y - d - dist, c[0][2] = z;
123 c[1][0]= x - dist, c[1][1]= y - d + dist, c[1][2] = z;
124 c[2][0]= x + dist, c[2][1]= y - d + dist, c[2][2] = z;
125 c[3][0]= x + dist, c[3][1]= y - d - dist, c[3][2] = z;
126 c[4][0]= x , c[4][1]= y - d , c[4][2] = z + R;
127 c[5][0]= x , c[5][1]= y - d, c[5][2] = z - R;
128 return c;
129 }
colloid_radii() const130 std::vector<double> colloid_radii()const {
131 int nsphere = colloid_coords().size();
132 std::vector<double> c(nsphere);
133 for(int i=0; i<nsphere; i++)
134 c[i] = R;
135 return c;
136 }
137
138 //make surface charge
make_surfcharge(const realfunc & u,const realfunc & surface,const realfunc & volume) const139 realfunc make_surfcharge(const realfunc& u,const realfunc& surface,const realfunc& volume) const {
140 real_derivative_3d Dx = free_space_derivative<double,3>(rhot.world(), 0);
141 real_derivative_3d Dy = free_space_derivative<double,3>(rhot.world(), 1);
142 real_derivative_3d Dz = free_space_derivative<double,3>(rhot.world(), 2);
143 realfunc Revolume = copy(volume);
144 Revolume.unaryop(Reciprocal<double,3>());
145 coord_3d lo(0.0),hi(0.0);
146 lo[0] = -20.0;
147 hi[0] = 20.0;
148 plot_line("Revolume.dat",401,lo,hi,Revolume);
149 realfunc gradu = Dx(u) + Dy(u) + Dz(u);
150 const double rfourpi = 1.0/(4.0*constants::pi);
151 realfunc ratio_surface_vol = binary_op(Revolume,surface, Bop());
152 realfunc scharge = binary_op(gradu,ratio_surface_vol, Bop());
153 return scharge.scale(rfourpi);
154 }
155 //The perturbed potential near the colloid as a result of the presence of the molecular charge distribution
perturbed_molecular_pot(const realfunc & surface,const realfunc & volume) const156 realfunc perturbed_molecular_pot(const realfunc& surface,const realfunc& volume) const {
157 const bool USE_SOLVER = true;
158 double tol = std::max(1e-4,FunctionDefaults<3>::get_thresh());
159 real_convolution_3d op = madness::CoulombOperator(world, tol*10.0, tol*0.1);
160 realfunc U0 = op(rhot);
161 coord_3d lo(0.0),hi(0.0);
162 lo[0] = -20.0;
163 hi[0] = 20.0;
164 plot_line("initial_pot.dat",401,lo,hi,U0);
165 realfunc charge = make_surfcharge(U0,surface,volume);
166 realfunc U = op(charge);
167 double unorm = U.norm2();
168 if (USE_SOLVER) {
169 madness::NonlinearSolver solver;//(5);
170 realfunc uvec, rvec;
171 if (world.rank() == 0){
172 print("\n\n");//for formating output
173 madness::print(" Computing the perturbed solute potential ");
174 madness::print(" ______________________ \n ");
175
176 madness::print("iteration "," "," residue norm2\n");
177 }
178 realfunc charg = make_surfcharge(U,surface,volume);
179 for (int iter=0; iter<maxiter; iter++) {
180 uvec = U;
181 //coord_3d lo(0.0),hi(0.0);
182 //lo[0] = -20.0;
183 // hi[0] = 20.0;
184 // plot_line("imp_Surfacepot.dat",1001,lo,hi,U);
185
186 rvec = (U -op(charg)).truncate() ;
187 realfunc U_new = solver.update(uvec,rvec);
188 double err = rvec.norm2();
189 if (world.rank()==0)
190 madness::print(" ", iter," " , err);
191 if (err >0.3*unorm){ U = 0.5*U + 0.5*U_new;
192 }
193 else
194 U = U_new;
195 if (err < 10.0*tol) break;
196 }
197 }
198 return U;
199 }
200
SurfaceMoleculeInteraction(const double & d,const double & R,const std::vector<double> & charge_center,realfunc & rhot,const double & sigma,World & world,const int & maxiter)201 SurfaceMoleculeInteraction(const double& d, const double& R,const
202 std::vector<double>& charge_center,
203 realfunc& rhot, const double& sigma, World& world,const int& maxiter)
204 :d(d),R(R), //d and R are in a.u
205 charge_center(0.0),
206 rhot(rhot),sigma(sigma),world(world),maxiter(maxiter){
207 MADNESS_ASSERT(colloid_coords().size()==colloid_radii().size());
208 }
209 };
210
main(int argc,char ** argv)211 int main(int argc, char **argv) {
212 initialize(argc, argv);
213 World world(SafeMPI::COMM_WORLD);
214 startup(world,argc,argv);
215
216 const int k = 6; // wavelet order
217
218 const double thresh = 1e-6; // truncation threshold
219
220 const double L = 50; // box is [-L,L]
221
222 //const int natom = 6; // number of atoms
223
224 double sigma = 0.5; // Surface width
225
226 const double R = 3.2503287; // radius of colloid sphere
227
228 const double d = 6.61404096; // distance between center of charge and coilloid center
229
230 const int maxiter = 25; // maximum iteration
231
232
233 // Function defaults
234
235 FunctionDefaults<3>::set_k(k);
236 FunctionDefaults<3>::set_thresh(thresh);
237 FunctionDefaults<3>::set_cubic_cell(-L, L);
238 FunctionDefaults<3>::set_initial_level(2);
239 //Creat an object
240 //Molecule molecule;
241 const std::vector<double> charge_center(0.0); //defines the center of nuclear charge
242 real_function_3d rhot = real_factory_3d(world).f(charge_function);
243 SurfaceMoleculeInteraction SMI(d,R,charge_center,rhot,sigma,world,maxiter);
244
245 real_functor_3d volume_functor(new MolecularVolumeMask(sigma, SMI.colloid_radii(), SMI.colloid_coords()));
246 real_functor_3d surface_functor(new MolecularSurface(sigma, SMI.colloid_radii(), SMI.colloid_coords()));
247 real_function_3d vol = real_factory_3d(world).functor(volume_functor);
248 real_function_3d surface = real_factory_3d(world).functor(surface_functor).truncate_on_project();
249 realfunc pot = SMI.perturbed_molecular_pot(surface,vol);
250 coord_3d lo(0.0),hi(0.0);
251 lo[0] = -20.0;
252 hi[0] = 20.0;
253 plot_line("colloid_vol.dat",1001,lo,hi,vol);
254 plot_line("colloid_surface.dat",1001,lo,hi,surface);
255 plot_line("colloid_pot.dat",1001,lo,hi,pot);
256
257 print("the volume is", vol.trace());
258 print("the area is", surface.trace());
259 finalize();
260 return 0;
261 }
262