madness/svpe__molecule__colloid_8h_source.html

 /*

   This file is part of MADNESS.


   Copyright (C) 2007,2010 Oak Ridge National Laboratory


   This program is free software; you can redistribute it and/or modify

   it under the terms of the GNU General Public License as published by

   the Free Software Foundation; either version 2 of the License, or

   (at your option) any later version.


   This program is distributed in the hope that it will be useful,

   but WITHOUT ANY WARRANTY; without even the implied warranty of

   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

   GNU General Public License for more details.


   You should have received a copy of the GNU General Public License

   along with this program; if not, write to the Free Software

   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA


   For more information please contact:


   Robert J. Harrison

   Oak Ridge National Laboratory

   One Bethel Valley Road

   P.O. Box 2008, MS-6367


   email: harrisonrj@ornl.gov

   tel:   865-241-3937

   fax:   865-572-0680


   $Id$

 */

 #ifndef SVPE_MOLECULE_COLLOID_H

 #define SVPE_MOLECULE_COLLOID_H


 #include <madness/mra/mra.h>

 #include <madness/mra/operator.h>

 #include <madness/mra/funcplot.h>

 #include <madness/tensor/solvers.h>

 #include <examples/molecularmask.h>

 #include <examples/nonlinsol.h>

 #include <madness/constants.h>

 #include <vector>


 using namespace madness;

 using namespace std;


 //set coordinates of colloid atoms

 std::vector< madness::Vector<double,3> > colloid_coords(const double&d,const double R,const std::vector<double> cc) {

     int nsphere = 1;

     std::vector< madness::Vector<double,3> > c(nsphere, madness::Vector<double,3>(0.0)); //6 spheres on the colloid surface

     /* double sqrttwo = std::sqrt(2.0);

        double dist= (sqrttwo/2.0)*R */

     double x = cc[0], y =  cc[1], z = cc[2];

 /*  c[0][0]= x - R,  c[0][1]= y - d - 2.0*R,   c[0][2] = z; //A  base sphere

     c[1][0]= x ,     c[1][1]= y - d - R,       c[1][2] = z;     //B base sphere

     c[2][0]= x + R,  c[2][1]= y - d + 2.0*R,   c[2][2] = z;  //C  base sphere

     c[3][0]= x ,     c[3][1]= y - d + 3.0*R,   c[3][2] = z;     //D base spheres

     c[4][0]= x ,     c[4][1]= y - d - 2.0*R ,  c[4][2] = z - R; //E bottom sphere

     c[5][0]= x ,     c[5][1]= y - d - 2.0*R,   c[5][2] = z + R;*/ // top sphere

     c[0][0]= x, c[0][1]= y - R - d ,    c[0][2] = z;// for a single sphere

     return c;

 }

 /*std::vector< madness::Vector<double,3> > colloid_coordss(const double&d,const double R,std::vector<double> cc) {

     std::vector< madness::Vector<double,3> > c(16,madness::Vector<double,3>(0,0)); //16 spheres on the colloid surface

     double sqrttwo = std::sqrt(2.0);

     double dist= (R/sqrttwo);

     double x = cc[0], y =  cc[1], z = cc[2];

     c[0][0]=  x - dist,           c[0][1]= y - dist,           c[0][2] = z-d-R; //A

     c[1][0]=  x - dist,           c[1][1]= y + dist,           c[1][2] = z-d-R; //B

     c[2][0]=  x + dist,           c[2][1]= y + dist,           c[2][2] = z-d-R;  //C

     c[3][0]=  x + dist,           c[3][1]= y - dist,           c[3][2] = z-d-R;   //D

     c[4][0]=  x + dist,           c[4][1]= y + 2.0*R - dist,   c[4][2] = z-d-R; //C'

     c[5][0]=  x + dist + 2.0*R,   c[5][1]= y + dist,           c[5][2] = z-d-R; //C"

     c[6][0]=  x + 2.0*R + dist ,  c[6][1]= y + 2.0*R + dist,   c[6][2] = z-d-R;//C'"

     c[7][0]=  x - dist ,          c[7][1]= y + 2.0*R + dist,   c[7][2] = z-d-R; //B'

     c[8][0]=  x - 2.0*R - dist ,  c[8][1]= y + dist,           c[8][2] = z-d-R; //B"

     c[9][0]=  x - 2.0*R - dist ,  c[9][1]= y + 2.0*R + dist,   c[9][2] = z-d-R; //B"'

     c[10][0]=  x - dist ,         c[10][1]= y - 2.0*R - dist,  c[10][2] = z-d-R; //A'

     c[11][0]=  x - 2.0*R - dist , c[11][1]= y  - 2.0*R - dist, c[11][2] = z-d-R; //A"

     c[12][0]=  x - 2.0*R - dist , c[12][1]= y - dist,          c[12][2] = z-d-R; //A"'

     c[13][0]=  x + 2.0*R + dist , c[13][1]= y - dist,          c[13][2] = z-d-R; //D'

     c[14][0]=  x + 2.0*R + dist , c[14][1]= y  - 2.0*R - dist, c[14][2] = z-d-R; //D"

     c[15][0]=  x + dist ,         c[15][1]= y - 2.0*R - dist,  c[15][2] = z-d-R; //D'"

     return c;

 }*/


 //colloid radii

 //if number of colloid spheres changes don't forget to change it here

 std::vector<double> colloid_radii(const double& R) {

     int nsphere = 1; //number of colloid spheres

     std::vector<double> c(nsphere,0.0);

     for(int i=0; i<nsphere; i++)

         c[i] = R;

     return c;

 }


 class SVPEColloidSolver {

     World& world;

     const double thresh;

     const double minlen;

     const double sigma;                 //< Width of surface layer

     const double epsilon_0;             //< Interior dielectric

     const double epsilon_1;             //< Exterior dielectric

     real_convolution_3d op;             //< Coulomb operator (1/r ... no 4pi)

     //std::vector<real_convolution_3d_ptr> gop; //< Gradient of the Coulomb operator

     vector_real_function_3d dlog; //< Log-derivative of the dielectric

     real_function_3d rdielectric; //< Reciprocal of the dielectric

     real_function_3d volume; //< volume function of the sphere/colloid

     double L;                            //<half lenght of the simulation volume

     //    vector_real_function_3d Eind; //<surface induced electric

     const int maxiter;

     static const double cutrho = 1e-12;  //cutoff value of the surface charge

     //this is used to perform a pointwise multiplication

     struct Bop {

         void operator()(const Key<3>& key,

                         real_tensor rfunc,

                         const real_tensor& func1,

                         const real_tensor& func2) const {

             ITERATOR(rfunc,

                      double d = func1(IND);

                      double p = func2(IND);

                      rfunc(IND) = d*p;

                      );

         }

         template <typename Archive>

         void serialize(Archive& ar) {}

     };

     real_function_3d func_pdt(const real_function_3d& func1,real_function_3d& func2) const {

         return binary_op(func1,func2, Bop());

     }

 public:

     SVPEColloidSolver(World& world,

                       double sigma, double epsilon_0, double epsilon_1,

                       const vector_real& atomic_radii, const vector_coord_3d& atomic_coords,

                       const double minlen,double L, const int maxiter)

         :world(world)

         ,thresh(FunctionDefaults<3>::get_thresh())

         , minlen(minlen)

         , sigma(sigma)

         , epsilon_0(epsilon_0)

         , epsilon_1(epsilon_1)

         , op(CoulombOperator(world, minlen, thresh))

         , dlog(3)

         ,L(L)

         , maxiter(maxiter)

     {

         MADNESS_ASSERT(atomic_radii.size() == atomic_coords.size());//check on the consistency

         // Functors for mask related quantities

         real_functor_3d rdielectric_functor(new MolecularVolumeExponentialSwitchReciprocal(sigma,epsilon_0,epsilon_1,atomic_radii,atomic_coords));

         real_functor_3d gradx_functor(new MolecularVolumeExponentialSwitchLogGrad(sigma, epsilon_0, epsilon_1, atomic_radii, atomic_coords,0));

         real_functor_3d grady_functor(new MolecularVolumeExponentialSwitchLogGrad(sigma, epsilon_0, epsilon_1, atomic_radii, atomic_coords,1));

         real_functor_3d gradz_functor(new MolecularVolumeExponentialSwitchLogGrad(sigma, epsilon_0, epsilon_1, atomic_radii, atomic_coords,2));

         real_functor_3d volume_functor(new MolecularVolumeMask(sigma,atomic_radii, atomic_coords));

         // Make the actual functions

         // const double rfourpi = 1.0/(4.0*constants::pi);

         rdielectric = real_factory_3d(world).functor(rdielectric_functor).nofence();

         volume = real_factory_3d(world).functor(volume_functor).nofence();

         dlog[0] = real_factory_3d(world).functor(gradx_functor).nofence().truncate_on_project();

         dlog[1] = real_factory_3d(world).functor(grady_functor).nofence().truncate_on_project();

         dlog[2] = real_factory_3d(world).functor(gradz_functor).truncate_on_project(); // FENCE

         // scale(world, dlog, rfourpi);

         rdielectric.truncate(false);

         truncate(world, dlog);

     }


     // Given the full Coulomb potential computes the surface charge

     real_function_3d make_surface_charge(const real_function_3d& u) const {

         real_derivative_3d Dx = free_space_derivative<double,3>(world, 0);

         real_derivative_3d Dy = free_space_derivative<double,3>(world, 1);

         real_derivative_3d Dz = free_space_derivative<double,3>(world, 2);

         //double fac = -1.0/(4.0*constants::pi);

         real_function_3d dx =  Dx(u);

         real_function_3d dy =  Dy(u);

         real_function_3d dz =  Dz(u);

         real_function_3d sc = func_pdt(dlog[0],dx).truncate() + func_pdt(dlog[1],dy).truncate() + func_pdt(dlog[2],dz).truncate();

         //coord_3d lo,hi;

         //lo[1]=-L, hi[1]=L;

         //plot_line("colloid_surf_charge.dat", 1001, lo, hi, sc);

         return sc; //.scale(fac);

     }

     // Given the full Laplace potential compute the surface charge

     real_function_3d make_Laplace_surface_charge(const real_function_3d& u,std::vector<double> E=std::vector<double>(3,0.0)) const {

         real_derivative_3d Dx = free_space_derivative<double,3>(world, 0);

         real_derivative_3d Dy = free_space_derivative<double,3>(world, 1);

         real_derivative_3d Dz = free_space_derivative<double,3>(world, 2);

         //double fac = -1.0/(4.0*constants::pi);

         real_function_3d dx = -E[0] + Dx(u);

         real_function_3d dy = -E[1] + Dy(u);

         real_function_3d dz = -E[2] + Dz(u);

         real_function_3d sc = func_pdt(dlog[0],dx).truncate() + func_pdt(dlog[1],dy).truncate() + func_pdt(dlog[2],dz).truncate();

         return sc;

     }

     //computes components of the the electric field due to the surface charge

     vector_real_function_3d make_electric_field(const real_function_3d& u) const {

         // throw "not correct";

         vector_real_function_3d E(3);// = madness::vector_factory(0.0, 0.0 ,0.0);

         real_derivative_3d Dx = free_space_derivative<double,3>(u.world(), 0);

         real_derivative_3d Dy = free_space_derivative<double,3>(u.world(), 1);

         real_derivative_3d Dz = free_space_derivative<double,3>(u.world(), 2);

         double fac = -1.0/(4.0*constants::pi);

         real_function_3d Sigmax =(Dx(u)).scale(fac), Sigmay = (Dy(u)).scale(fac),Sigmaz = (Dz(u)).scale(fac); //excess charge on colloid surface

         coord_3d lo(0.0),hi(0.0);

         lo[1]=-L, hi[1]=L;

         plot_line("ab_sigma.dat", 1001, hi, lo, Sigmax, Sigmay, Sigmaz);

         E[0] = Dx(op(Sigmax)), E[1] = Dy(op(Sigmay)), E[2] = Dz(op(Sigmaz)) ;

         return E;

     }


     real_function_3d solve(const real_function_3d& rho) const {

         const double fac = -1.0/(4.0*constants::pi);

         real_function_3d charge = (rdielectric*rho).scale(-1.0);


         // Initial guess is constant dielectric

         real_function_3d u0 = op(charge);//.truncate();

         //real_function_3d u = uguess.is_initialized() ? uguess : u0;

         real_function_3d u = u0; // op(rho);  ?? RJH

         double unorm = u.norm2();

         if (world.rank()==0)

             print("UNORM IS", unorm);

         NonlinearSolver solver;

         if (world.rank()==0){

             print("\n\n");//for formating output

             madness::print("              Computing the Perturbed Potential (due to molecule)        ");

             madness::print("                         Near the Colloid Surface                         ");

             madness::print("                          ______________________                       \n ");

         }

         for (int iter=0; iter<maxiter; iter++) {

             double start = wall_time();

             real_function_3d surface_charge = make_surface_charge(u);

             real_function_3d r = (u - u0 - op(surface_charge).scale(fac)).truncate(.032*FunctionDefaults<3>::get_thresh());

             double sigtot = surface_charge.trace()*fac;

             //surface_charge.clear();

             // u0.clear();

             real_function_3d unew = solver.update(u,r);

             double change = (unew-u).norm2();//volume.trace(); //times 1000 when not in solvent

             if(world.rank()==0){

                 print("iter", iter, "change", change,

                       "soln(10.0)", u(coord_3d(10.0)),

                       "surface charge", sigtot,"used",wall_time()-start);


             }

             // Step restriction

             if (change > 0.3*unorm)

                 u = 0.5*unew + 0.5*u;

             else

                 u = unew;

             //if (change < std::max(1e-4,10.0*thresh/volume.trace())) break;

             if (change < 10.0*thresh) break;

         }

         if(world.rank()==0)

             print("\n\n");//format printing

         real_function_3d urxn = op(make_surface_charge(u)).scale(fac);//u - op(rho);

         coord_3d lo,hi;

         lo[1]=-L, hi[1]=L;

         plot_line("colloid_rxn_pot.dat", 1001, lo, hi, urxn);

         plot_line("colloid_surfcharge.dat", 1001, lo, hi, make_surface_charge(urxn).scale(fac));

         //real_tensor cell(3,2);

         // cell(_,0) = -L;

         //cell(_,1) =  L;

         // plotdx(volume, "testpot.dx", cell);

         return urxn;

     }

     real_function_3d solve_Laplace(std::vector<double>E) const {

         // Initial guess is constant dielectric

         const double fac = -1.0/(4.0*constants::pi);

         real_function_3d u(world);//  guess pot is zero;

         double unorm = 0.0;//u.norm2();

         NonlinearSolver solver;

         //print for formating

         if (world.rank()==0){

             print("\n\n");//for formating output

             madness::print("            Computing the Perturbed Potential (due to external field)           ");

             madness::print("                           ______________________                            \n ");

         }

         for (int iter=0; iter<maxiter; iter++) {

             double start = wall_time();

             real_function_3d surface_charge = make_Laplace_surface_charge(u,E);

             real_function_3d r = (u - op(surface_charge).scale(fac)).truncate(.032*FunctionDefaults<3>::get_thresh());

             double sigtot = surface_charge.trace()*fac;

             //surface_charge.clear();

             real_function_3d unew = solver.update(u, r);

             double change = (unew-u).norm2();//volume.trace();///(8.0*std::pow(L,3.0));

             if (world.rank()==0){

                 print("iter", iter, "change", change,

                       "soln(10.0)", u(coord_3d(10.0)),

                       "surface charge", sigtot,"used",wall_time()-start);

             }

             // Step restriction

             if (change > 0.3*unorm)

                 u = 0.5*unew + 0.5*u;

             else

                 u = unew;


             if (change < std::min(1e-4,10.0*thresh)) break;

         }

         if (world.rank()==0)

             print("\n\n");

         coord_3d lo,hi;

         lo[1]=-L, hi[1]= L;

         plot_line("laplace_surfcharge.dat", 1001, lo, hi, make_Laplace_surface_charge(u,E).scale(fac));

         plot_line("laplace_pot.dat", 1001, lo, hi, u);

         return u;

     }

     //calculate the average reaction field(\int C(r)F_r(r)d \tau/\int C(r)d\tau

     //the mask is that of the molecule because the average field is that felt by the molecule

     double ave_rxn_field(const real_function_3d& u,const real_function_3d& mask)const {

         real_function_3d  pdtx = mask*make_electric_field(u)[0];

         real_function_3d  pdty = mask*make_electric_field(u)[1];

         real_function_3d  pdtz = mask*make_electric_field(u)[2];

         double numx = pdtx.trace(), numy = pdty.trace(), numz = pdtz.trace();

         double denominator = mask.trace();

         double Favx = numx/denominator, Favy = numy/denominator,Favz = numz/denominator;

         return std::sqrt(std::pow(Favx,2.0) + std::pow(Favy,2.0) + std::pow(Favz,2.0));

     }

 };


 #endif


thresh
const double thresh
Definition: dielectric.cc:185

madness::NonlinearSolverND
A simple Krylov-subspace nonlinear equation solver.
Definition: nonlinsol.h:91

ITERATOR
#define ITERATOR(t, exp)
Definition: tensor_macros.h:249

SVPEColloidSolver::solve
real_function_3d solve(const real_function_3d &rho) const
Solve for the full Coulomb potential using the free-particle GF.
Definition: svpe_molecule_colloid.h:211

madness::tr1::shptr::shared_ptr
Definition: shared_ptr_bits.h:38

madness::constants::pi
const double pi
Mathematical constant pi.
Definition: constants.h:44

SVPEColloidSolver::solve_Laplace
real_function_3d solve_Laplace(std::vector< double >E) const
Solve for the full Coulomb potential using the free-particle GF.
Definition: svpe_molecule_colloid.h:266

mra.h
Main include file for MADNESS and defines Function interface.

madness::Function::truncate
Function< T, NDIM > & truncate(double tol=0.0, bool fence=true)
Truncate the function with optional fence. Compresses with fence if not compressed.
Definition: mra.h:577

IND
#define IND
Definition: tensor_macros.h:204

madness::FunctionDefaults
FunctionDefaults holds default paramaters as static class members.
Definition: funcdefaults.h:175

epsilon_1
const double epsilon_1
Definition: dielectric.cc:190

madness::Vector< double, 3 >

madness::Derivative
Implements derivatives operators with variety of boundary conditions on simulation domain...
Definition: derivative.h:272

sigma
const double sigma
Definition: dielectric.cc:187

L
const double L
Definition: 3dharmonic.cc:123

madness::truncate
void truncate(World &world, std::vector< Function< T, NDIM > > &v, double tol=0.0, bool fence=true)
Truncates a vector of functions.
Definition: vmra.h:194

constants.h
Defines common mathematical and physical constants.

madness::SeparatedConvolution
Convolutions in separated form (including Gaussian)
Definition: operator.h:117

colloid_coords
std::vector< madness::Vector< double, 3 > > colloid_coords(const double &d, const double R, const std::vector< double > cc)
Definition: svpe_molecule_colloid.h:49

SVPEColloidSolver::ave_rxn_field
double ave_rxn_field(const real_function_3d &u, const real_function_3d &mask) const
Definition: svpe_molecule_colloid.h:309

std
Definition: mpreal.h:3066

madness::coord_3d
Vector< double, 3 > coord_3d
Definition: funcplot.h:634

epsilon_0
const double epsilon_0
Definition: dielectric.cc:189

MolecularVolumeMask
Definition: molecularmask.h:131

colloid_radii
std::vector< double > colloid_radii(const double &R)
Definition: svpe_molecule_colloid.h:90

solvers.h
Defines interfaces for optimization and non-linear equation solvers.

MolecularVolumeExponentialSwitchReciprocal
Computes the reciprocal of MolecularVolumeExponentialSwitch.
Definition: molecularmask.h:233

madness::vector_real_function_3d
std::vector< real_function_3d > vector_real_function_3d
Definition: functypedefs.h:79

madness::vector_coord_3d
std::vector< Vector< double, 3 > > vector_coord_3d
Definition: functypedefs.h:58

operator.h
Implements most functionality of separated operators.

madness::Function::world
World & world() const
Returns the world.
Definition: mra.h:622

SVPEColloidSolver::make_Laplace_surface_charge
real_function_3d make_Laplace_surface_charge(const real_function_3d &u, std::vector< double > E=std::vector< double >(3, 0.0)) const
Definition: svpe_molecule_colloid.h:183

madness::Function::norm2
double norm2() const
Returns the 2-norm of the function ... global sum ... works in either basis.
Definition: mra.h:650

madness::scale
void scale(World &world, std::vector< Function< T, NDIM > > &v, const std::vector< Q > &factors, bool fence=true)
Scales inplace a vector of functions by distinct values.
Definition: vmra.h:290

madness::NonlinearSolverND::update
Function< double, NDIM > update(const Function< double, NDIM > &u, const Function< double, NDIM > &r, const double rcondtol=1e-8, const double cabsmax=1000.0)
Computes next trial solution vector.
Definition: nonlinsol.h:111

nonlinsol.h
Example implementation of Krylov-subspace nonlinear equation solver.

madness::Function::trace
T trace() const
Returns global value of int(f(x),x) ... global comm required.
Definition: mra.h:1034

mpfr::min
const mpreal min(const mpreal &x, const mpreal &y)
Definition: mpreal.h:2675

cc
const double cc
Definition: navstokes_cosines.cc:108

madness::real_tensor
Tensor< double > real_tensor
Definition: functypedefs.h:40

sqrt
tensorT sqrt(const tensorT &s, double tol=1e-8)
Computes matrix square root (not used any more?)
Definition: DFcode/moldft.cc:446

madness::World
A parallel world with full functionality wrapping an MPI communicator.
Definition: worldfwd.h:416

madness::archive::serialize
enable_if_c< is_serializable< T >::value &&is_output_archive< Archive >::value >::type serialize(const Archive &ar, const T *t, unsigned int n)
Definition: archive.h:615

madness::norm2
double norm2(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norm of a vector of functions.
Definition: vmra.h:325

SVPEColloidSolver::make_electric_field
vector_real_function_3d make_electric_field(const real_function_3d &u) const
Definition: svpe_molecule_colloid.h:195

maxiter
const int maxiter
Definition: gygi_soltion.cc:69

madness::Function
A multiresolution adaptive numerical function.
Definition: derivative.h:61

SVPEColloidSolver
Definition: svpe_molecule_colloid.h:98

SVPEColloidSolver::make_surface_charge
real_function_3d make_surface_charge(const real_function_3d &u) const
Definition: svpe_molecule_colloid.h:168

madness::World::rank
ProcessID rank() const
Returns the process rank in this world (same as MPI_Comm_rank()))
Definition: worldfwd.h:526

madness::plot_line
void plot_line(const char *filename, int npt, const Vector< double, NDIM > &lo, const Vector< double, NDIM > &hi, const Function< T, NDIM > &f)
Generates ASCII file tabulating f(r) at npoints along line r=lo,...,hi.
Definition: funcplot.h:388

molecularmask.h

madness::vector_real
std::vector< double > vector_real
Definition: functypedefs.h:53

funcplot.h
Defines/implements plotting interface for functions.

op
Tensor< double > op(const Tensor< double > &x)
Definition: kain.cc:508

madness::print
void print(const A &a)
Print a single item to std::cout terminating with new line.
Definition: print.h:122

madness::real_factory_3d
FunctionFactory< double, 3 > real_factory_3d
Definition: functypedefs.h:93

MolecularVolumeExponentialSwitchLogGrad
Returns the requested component of the derivative of the log of MolecularVolumeExponentialSwitch.
Definition: molecularmask.h:295

madness
Holds machinery to set up Functions/FuncImpls using various Factories and Interfaces.
Definition: chem/atomutil.cc:45

c
const double c
Definition: gfit.cc:200

madness::binary_op
Function< TENSOR_RESULT_TYPE(L, R), NDIM > binary_op(const Function< L, NDIM > &left, const Function< R, NDIM > &right, const opT &op, bool fence=true)
Generate new function = op(left,right) where op acts on the function values.
Definition: mra.h:1593

madness::wall_time
double wall_time()
Returns the wall time in seconds relative to arbitrary origin.
Definition: world.cc:248

madness::Key
Key is the index for a node of the 2^NDIM-tree.
Definition: key.h:69

SVPEColloidSolver::SVPEColloidSolver
SVPEColloidSolver(World &world, double sigma, double epsilon_0, double epsilon_1, const vector_real &atomic_radii, const vector_coord_3d &atomic_coords, const double minlen, double L, const int maxiter)
Definition: svpe_molecule_colloid.h:133