# h2plus.R
#
# Find the approximate ground state energy, and the equilibrium
# internuclear separation for the H2+ molecular ion.
#
# Copyright (c) 2015 Krishna Myneni, http://ccreweb.org
#
# This code may be used for any purpose as long as the
# copyright notice above is preserved.
#
# Revisions:
#   2015-01-16  km; first working version.
#   2015-01-17  km; added example and reference, and revised
#                   the notes. Simplification of the integration
#                   needs to be considered.
#   2015-01-22  km; added ref. [2] and revised notes.
#
# Notes:
#
#  1) All units used here are atomic units.
#
#  2) The Hamiltonian for the H2+ ion is given in cylindrical
#     coordinates by,
#
#        H = -hbar^2/(2*mu) * Del^2 
#	     - 1/sqrt(rho^2 + (z+a/2)^2) 
#            - 1/sqrt(rho^2 + (z-a/2)^2)
#	     + 1/a
#
#    where Del^2 is the Laplacian operator, rho and z are
#    the electron coordinates, and the nuclei are on the
#    z-axis at z = +/- a/2. In atomic units, hbar = 1,
#    and m_e =1, where m_e is the mass of the electron.
#    We will take mu = m_e.
#
# 3) We take for the trial wavefunction of the H2+ molecule
#    ground state, the function,
#
#        Rt(rho, z) = exp(-alpha*sqrt(rho^2 + (z-a/2)^2)) 
#                     +
#                     exp(-alpha*sqrt(rho^2 + (z+a/2)^2))
#
#    This function is simply the sum of two ground state
#    hydrogen atom-like wavefunctions, separated by a distance,
#    "a". The function is expressed in cylindrical coordinates.
#
# 4) The variational principle guarantees that, for a continuous
#    trial wavefunction, the minimum energy with respect to the
#    adjustable parameters in the function will always be at or
#    above the actual energy. Our trial wavefunction for the 
#    molecule's ground state is based on the symmetry of the
#    molecule, and has two adjustable parameters: 
#
#       a     -- the internuclear spacing
#       alpha -- the inverse decay length of the wavefunction,
#                or, equivalently, an effective nuclear charge.
#
#    The procedure for finding the ground state energy and
#    the equilibrium separation of the nuclei is to compute
#    <H>, the expectation value of the total energy, for
#    a given pair of parameters (alpha, a), and varying
#    the parameters until the minimum is found for <H>.
#
#    In standard introductory treatments of the H2+ molecule [1],
#    the trial wavefunction is a linear combination of H-atom
#    wavefunctions, i.e. with alpha fixed at 1. Such a trial
#    function gives a less accurate estimate of the ground state
#    energy and the equilibrium internuclear separation.
#    The calculations of ref. [1] may be reproduced by this
#    program by fixing alpha to 1, and varying "a". This program
#    uses the same trial function as the lowest accuracy
#    calculation of ref. [2].
#
#    Computing <H> for general trial wavefunctions requires
#    evaluation of difficult three-dimensional integrals 
#    (two dims in our case; because of the symmetry of the 
#    problem, integration over the coordinate phi is trivial
#    for our trial function). In this program, the 2D integrals
#    are evaluated numerically using the R package, "cubature".
#
# 5) An example of using this code is given at the bottom.
#
# References:
#
#  1. http://farside.ph.utexas.edu/teaching/qmech/Quantum/node129.html
#
#  2. H.W. Jones and B. Etemadi, Accurate ground-state
#     calculations of H2^+ using basis sets of atom-centered
#     Slater-type orbitals, Physical Review A, vol. 47, p. 3430
#     (1993).
#      
require(cubature)

LN_MIN = -log(1e-15)

# The trial wavefunction parameters, alpha and "a", are
# global.

alpha = 2
a = 2

# Some useful functions of rho and z, for defining our
# trial wavefunction and its partial derivatives.
f_m <- function( rho, z ) { sqrt( rho^2 + (z - a/2)^2 ) }
f_p <- function( rho, z ) { sqrt( rho^2 + (z + a/2)^2 ) }

u_m <- function( rho, z ) {
	fm = f_m(rho, z)
	return( exp(-alpha*fm)/fm )
}
u_p <- function( rho, z ) {
	fp = f_p(rho, z)
	return( exp(-alpha*fp)/fp )
}

# Trial wavefunction and its square
Rt <- function( rho, z ) { 
	exp(-alpha*f_m(rho, z)) + exp(-alpha*f_p(rho, z))
}

Rt2 <- function( rho, z ) { (Rt(rho, z))^2  }


# Partial first derivatives of Rt and u_m and u_p
D_rho_Rt <- function( rho, z ) { -alpha*rho*(u_m(rho, z) + u_p(rho,z)) }

D_z_Rt <- function( rho, z ) {
	-alpha*((z-a/2)*u_m(rho, z) + (z+a/2)*u_p(rho, z))
}

D_rho_u_m <- function( rho, z ) {
	fm = f_m(rho, z)
	um = exp(-alpha*fm)/fm
	return( -rho*um*(alpha + 1/fm)/fm )
}

D_rho_u_p <- function( rho, z ) {
	fp = f_p(rho, z)
	up = exp(-alpha*fp)/fp
	return( -rho*up*(alpha + 1/fp)/fp )
}

D_z_u_m <- function( rho, z ) {
	fm = f_m(rho, z)
	um = exp(-alpha*fm)/fm
	return( -(z-a/2)*um*(alpha + 1/fm)/fm )
}

D_z_u_p <- function( rho, z ) {
	fp = f_p(rho, z)
	up = exp(-alpha*fp)/fp
	return( -(z+a/2)*up*(alpha + 1/fp)/fp )
}

# Partial second derivatives of Rt w.r.t. rho and z

D2_rho_Rt <- function( rho, z ) {
	fm = f_m(rho, z)
	fp = f_p(rho, z)
	um = exp(-alpha*fm)/fm
	up = exp(-alpha*fp)/fp
	Drho_um = -rho*um*(alpha+1/fm)/fm
	Drho_up = -rho*up*(alpha+1/fp)/fp
	return( -alpha*(rho*(Drho_um + Drho_up) + um + up) )
}
	
D2_z_Rt <- function( rho, z ) {
	fm = f_m(rho, z)
	fp = f_p(rho, z)
	um = exp(-alpha*fm)/fm
	up = exp(-alpha*fp)/fp
	Dz_um = -(z-a/2)*um*(alpha+1/fm)/fm
	Dz_up = -(z+a/2)*up*(alpha+1/fp)/fp
	return( -alpha*((z-a/2)*Dz_um + (z+a/2)*Dz_up + um + up) )
}


# Integrands for normalization, average kinetic energy, and
#   average potential energy
integ_N <- function(x) { x[1]*Rt2(x[1],x[2]) }

integ_K1 <- function(x) { x[1]*Rt(x[1],x[2])*D2_rho_Rt(x[1],x[2]) }
integ_K2 <- function(x) { x[1]*Rt(x[1],x[2])*D2_z_Rt(x[1],x[2]) }
integ_K3 <- function(x) { Rt(x[1],x[2])*D_rho_Rt(x[1],x[2]) }

integ_V1 <- function(x) { x[1]*Rt2(x[1],x[2]) / f_p(x[1],x[2]) }
integ_V2 <- function(x) { x[1]*Rt2(x[1],x[2]) / f_m(x[1],x[2]) }

Integrate_rho_z <- function( integrand ) {
	rho_max = LN_MIN/alpha
	z_max = rho_max
	Intg = adaptIntegrate( integrand, tol=1e-7,
		lowerLimit = c(0,0), upperLimit = c(rho_max, z_max))$integral
	return( Intg )
}


E_av <- function() {
        # Compute normalization constant for the trial wavefunction
	N = 4*pi*Integrate_rho_z( integ_N )

	# Compute the integrals for <K>, the expectation
	# value of the kinetic energy
	I_K1 = Integrate_rho_z( integ_K1 )
	I_K2 = Integrate_rho_z( integ_K2 )
	I_K3 = Integrate_rho_z( integ_K3 )

	K_av = -2*pi*(I_K1 + I_K2 + I_K3)/N

	# Compute the integrals for <V_el>, the expectation
	# value of the electronic potential energy
        I_V1 = Integrate_rho_z( integ_V1 )
	I_V2 = Integrate_rho_z( integ_V2 )

	V_av = -4*pi*(I_V1 + I_V2)/N

	# Coulomb repulsion energy from nuclei
        V_c = 1/a

	return( c(K_av, V_av, V_c) )
}

# Find the expectation values of energies over a grid of the
# parameter values for "a" and "alpha". The arguments are a
# sequence of alpha values and a sequence of "a" values. The
# matrix of  energies with alpha varying by row, and "a" 
# varying by column is returned.

FindEnergies <- function( alpha_s, a_s ) {
	n_alpha = length(alpha_s)
	n_a = length(a_s)

	energies = matrix( rep(0, n_alpha*n_a), nrow=n_alpha,
			ncol=n_a, byrow=TRUE )

	for (i in 1:n_alpha) {
	   alpha <<- alpha_s[i]
	   for (j in 1:n_a) {	
		a <<- a_s[j]
		energies[i,j] = sum(E_av())
	   }
	}
	return( energies )
}


# Example:
# This example may require a few minutes of computation time.
#
# > alpha_vals = seq(1.0, 1.5, by=0.05)
# > a_vals = seq(1.8, 2.2, by=0.05)
# > emat = FindEnergies(alpha_vals, a_vals)
# > min(emat)
# >
# > require(plot3D)
# > persp3D( x = alpha_vals, y = a_vals, z = emat)