*  Internuclear potential energy calculation for H2+ ion ground level
close all
clear all
set century on
set confirm on
clear
set date german
set decimals to 18
set exact on
set safety off
set status off
set sysmenu off
set talk off
set typeahead to 3000
on shutdown quit
on escape quit
set alternate to "INPE_Res.txt"
set alternate on 
TEXT

  Internuclear potential energy (electronic energy including internuclear
  repulsion) for hydrogen molecular monovalent cation ground state
 (Calculation of the same for the first excited states is too error-prone)

  Formulae used (Quantum Chemistry by Ira N. Levine, Chap. 13, Sec. MO theory of H2+):

  Overlap integral    Sab = (1 + k*R + k*k*R*R/3) * exp((-1)*k*R)
  Coulomb integral    Haa = 0.5*k*k - k - 1/R + (k+1/R) * exp((-2)*k*R)
  Resonance Integral  Hab = (-0.5)*k*k*Sab - k*(2-k)*(1+k*R)*exp((-1)*k*R)
  Internuclear Potn.  U   = (Haa + Hab)/(1 + Sab) + 1/R
  (Electronic energy)


  Value of the orbital exponent k for each value of R is obtained as the
  optimum value that gives the lowest value of U (variation theorem).

  The results of this work will be stored in INPE_Res.txt in this folder.

ENDTEXT
SET ALTERNATE OFF
wait "  Would you include values of Sab, Haa, Hab & U in INPE_Res.txt [Y,N]? ";
to openness
@ 1,0 clear
R = 0.4
SET ALTERNATE ON
If upper(openness) = "Y"
  TEXT

 Internuclear   Orbital      Overlap     Coulomb    Resonance  Electronic
  Distance R   Exponent k   Integral    Integral    Integral    Energy U
 (AtomicUnit) (PureNumber)  Sab(a.u.)   Haa(a.u.)   Hab(a.u.)    (a. u.)

  ENDTEXT
Else
  TEXT

 Internuclear   Orbital      Overlap     Coulomb    Resonance  Electronic
  Distance R   Exponent k   Integral    Integral    Integral    Energy U
 (AtomicUnit) (PureNumber)  Sab(a.u.)   Haa(a.u.)   Hab(a.u.)    (a. u.)

  ENDTEXT
Endif
SET ALTERNATE OFF
DO WHILE R <= 8.0
k = 0.1
EP_MIN = 1000
DO WHILE k <= 2
  NUME = k*k - k - 1/R + (1+k*R)*EXP((-2)*k*R)/R + k*(k-2)*(1+k*R) ;
  *EXP((-1)*k*R)
  DENO = 1 + EXP((-1)*k*R)*(1+k*R+k*k*R*R/3)
  EP = (-0.5)*k*k + NUME/DENO + 1/R
  IF EP < EP_MIN
    k_min = k
    EP_MIN = EP
  ENDIF
  k = k + 0.01
ENDDO
k = k_min - 0.01
k_lim = k_min + 0.01
DO WHILE k <= k_lim
  NUME = k*k - k - 1/R + (1 + k*R)*EXP((-2)*k*R)/R + k*(k-2)*(1 + k*R) ;
  *EXP((-1)*k*R)
  DENO = 1 + EXP((-1)*k*R)*(1 + k*R + k*k*R*R/3)
  EP = (-0.5)*k*k + NUME/DENO + 1/R
  IF EP < EP_MIN
    k_min = k
    EP_MIN = EP
  ENDIF
  k = k + 0.0001
ENDDO
k = k_min
EP = EP_MIN
if k < 1.0
  k = 1.0
endif
Sab = (1 + k*R + k*k*R*R/3) * exp((-1)*k*R)
Haa = 0.5*k*k - k - 1/R + (k+1/R) * exp((-2)*k*R)
Hab = (-0.5)*k*k*Sab - k*(2-k)*(1+k*R) * exp((-1)*k*R)
U = (Haa + Hab)/(1 + Sab) + 1/R
SET ALTERNATE ON
If upper(openness) = "Y"
  ?str(round(R,5),11,5),str(round(k,5),11,5),str(round(Sab,5),11,5),str(round(Haa,5),12,5), ;
  str(round(Hab,5),11,5),str(round(U,5),11,5)
Else
  ?str(round(R,5),11,5),str(round(k,5),11,5)
Endif
SET ALTERNATE OFF
R = R + 0.2
ENDDO
SET ALTERNATE ON
?
SET ALTERNATE OFF
close alternate 
close all
clear all
quit