1 /*
2   This file is part of MADNESS.
3 
4   Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6   This program is free software; you can redistribute it and/or modify
7   it under the terms of the GNU General Public License as published by
8   the Free Software Foundation; either version 2 of the License, or
9   (at your option) any later version.
10 
11   This program is distributed in the hope that it will be useful,
12   but WITHOUT ANY WARRANTY; without even the implied warranty of
13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14   GNU General Public License for more details.
15 
16   You should have received a copy of the GNU General Public License
17   along with this program; if not, write to the Free Software
18   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20   For more information please contact:
21 
22   Robert J. Harrison
23   Oak Ridge National Laboratory
24   One Bethel Valley Road
25   P.O. Box 2008, MS-6367
26 
27   email: harrisonrj@ornl.gov
28   tel:   865-241-3937
29   fax:   865-572-0680
30 
31   $Id$
32 */
33 /*!
34   \file he.cc
35   \brief Solves the Hartree-Fock equations for the helium atom
36   \defgroup examplehehf Hartree-Fock equations for the helium atom
37   \ingroup examples
38 
39   The source is <a href=http://code.google.com/p/m-a-d-n-e-s-s/source/browse/local/trunk/src/apps/examples/he.cc>here</a>.
40 
41   \par Points of interest
42   - application of the Coulomb and Helmholtz Green's functions
43   - smoothing of the potential and initial guess
44   - manual evaluation of the solution along a line
45 
46   \par Background
47 
48   The Hartree-Fock wave function is computed for the helium atom
49   in three dimensions without using spherical symmetry.
50 
51   The atomic orbital is an eigenfunction of the Fock operator
52   \f{eqnarray*}{
53      \hat{F} \phi(r) &=& \epsilon \phi(r) \\
54      \hat{F} &=& -\frac{1}{2} \nabla^2 - \frac{2}{r} + u(r) \\
55      u(r) &=& \int \frac{\rho(s)}{| r - s |} d^3s \\
56      \rho(r) &=& \phi(r)^2
57   \f}
58   that depends upon the orbital via the Coulomb potential (\f$ u(r) \f$)
59   arising from the electronic density (\f$ \rho(r) \f$).
60 
61   \par Implementation
62 
63   Per the usual MADNESS practice, the equation is rearranged into
64   integral form
65   \f[
66       \phi = - 2 G_{\mu} * ( V \phi)
67   \f]
68   where \f$ \mu = \sqrt{-2E} \f$ and \f$G_{\mu}\f$ is the Green's function
69   for the Helmholtz equation
70   \f[
71       \left( - \nabla^2 + \mu^2  \right) G(r,r^{\prime}) = \delta(r,r^{\prime})
72   \f]
73 
74   The initial guess is \f$ \exp(-2r) \f$, which is normalized before use.
75   Each iteration proceeds by
76   - computing the density by squaring the orbital,
77   - computing the Coulomb potential by applying the Coulomb Green's function,
78   - multiplying the orbital by the total potential,
79   - updating the orbital by applying the Helmholtz Green's function,
80   - updating the energy according to a second-order accurate estimate
81     (see the initial MRQC paper), and finally
82   - normalizing the new wave function.
83 
84   The kinetic energy operator is denoted by \f$T~=~-\frac{1}{2}
85   \nabla^2 \f$.  Thus, one would expect to compute the kinetic energy
86   with respect to a wave function \f$ \phi \f$ by \f$ <T \phi, \phi >
87   \f$.  In this particular example, the wave functions goes to zero
88   smoothly before the boundary so we apply the product rule for differentiation
89   and the kinetic energy is equal to \f$ < \nabla \phi, \nabla \phi >
90   \f$.
91 
92 
93 */
94 
95 //#define WORLD_INSTANTIATE_STATIC_TEMPLATES
96 #include <madness/mra/mra.h>
97 #include <madness/mra/operator.h>
98 
99 using namespace madness;
100 
101 static const double L = 32.0;   // box size
102 static const long k = 8;        // wavelet order
103 static const double thresh = 1e-6; // precision
104 
guess(const coord_3d & r)105 static double guess(const coord_3d& r) {
106     const double x=r[0], y=r[1], z=r[2];
107     return 6.0*exp(-2.0*sqrt(x*x+y*y+z*z+1e-4));
108 }
109 
V(const coord_3d & r)110 static double V(const coord_3d& r) {
111     const double x=r[0], y=r[1], z=r[2];
112     return -2.0/(sqrt(x*x+y*y+z*z+1e-8));
113 }
114 
iterate(World & world,real_function_3d & V,real_function_3d & psi,double & eps)115 void iterate(World& world, real_function_3d& V, real_function_3d& psi, double& eps) {
116     real_function_3d Vpsi = (V*psi);
117     Vpsi.scale(-2.0).truncate();
118     real_convolution_3d op = BSHOperator3D(world, sqrt(-2*eps), 0.001, 1e-6);
119     real_function_3d tmp = apply(op,Vpsi).truncate();
120     double norm = tmp.norm2();
121     real_function_3d r = tmp-psi;
122     double rnorm = r.norm2();
123     double eps_new = eps - 0.5*inner(Vpsi,r)/(norm*norm);
124     if (world.rank() == 0) {
125         print("norm=",norm," eps=",eps," err(psi)=",rnorm," err(eps)=",eps_new-eps);
126     }
127     psi = tmp.scale(1.0/norm);
128     eps = eps_new;
129 }
130 
main(int argc,char ** argv)131 int main(int argc, char** argv) {
132     initialize(argc, argv);
133     World world(SafeMPI::COMM_WORLD);
134     startup(world,argc,argv);
135     std::cout.precision(6);
136 
137     FunctionDefaults<3>::set_k(k);
138     FunctionDefaults<3>::set_thresh(thresh);
139     FunctionDefaults<3>::set_truncate_mode(1);
140     FunctionDefaults<3>::set_cubic_cell(-L/2,L/2);
141 
142     real_function_3d Vnuc = real_factory_3d(world).f(V).truncate_mode(0);
143     real_function_3d psi  = real_factory_3d(world).f(guess);
144     psi.scale(1.0/psi.norm2());
145     real_convolution_3d op = CoulombOperator(world, 0.001, 1e-6);
146 
147     double eps = -1.0;
148     for (int iter=0; iter<10; iter++) {
149         real_function_3d rho = square(psi).truncate();
150         real_function_3d potential = Vnuc + apply(op,rho).truncate();
151         iterate(world, potential, psi, eps);
152     }
153 
154     double kinetic_energy = 0.0;
155     for (int axis=0; axis<3; axis++) {
156         real_derivative_3d D = free_space_derivative<double,3>(world, axis);
157         real_function_3d dpsi = D(psi);
158         kinetic_energy += inner(dpsi,dpsi);
159     }
160 
161     real_function_3d rho = square(psi).truncate();
162     double two_electron_energy = inner(apply(op,rho),rho);
163     double nuclear_attraction_energy = 2.0*inner(Vnuc*psi,psi);
164     double total_energy = kinetic_energy + nuclear_attraction_energy + two_electron_energy;
165 
166     // Manually tabluate the orbital along a line ... probably easier
167     // to use the lineplot routine
168     coord_3d r(0.0);
169     psi.reconstruct();
170     for (int i=0; i<201; i++) {
171         r[2] = -L/2 + L*i/200.0;
172         print(r[2], psi(r));
173     }
174 
175     if (world.rank() == 0) {
176         print("            Kinetic energy ", kinetic_energy);
177         print(" Nuclear attraction energy ", nuclear_attraction_energy);
178         print("       Two-electron energy ", two_electron_energy);
179         print("              Total energy ", total_energy);
180         print("                    Virial ", (nuclear_attraction_energy + two_electron_energy) / kinetic_energy);
181     }
182 
183     world.gop.fence();
184 
185     finalize();
186     return 0;
187 }
188