#define calculator
calculator = LCAOCalculator()

#initialize arrays
energies = numpy.array([])
old_configuration = None

#setup distances where h2 energy is evaluated
distances = numpy.linspace(0.7,0.8,6)

#calculate total energy as function of distance
for d in distances:
    #setup H-H geometry
    elements = [ Hydrogen, Hydrogen]
    cartesian_coordinates = [[ 0.0,0.0,0.0],
                       [ 0.0,0.0,d]]*Angstrom
    molecule_configuration = MoleculeConfiguration(
        elements=elements,
        cartesian_coordinates=cartesian_coordinates
        )
    # set the calculator
    if old_configuration == None:
        molecule_configuration.setCalculator(calculator())
    else:
        molecule_configuration.setCalculator(calculator(),
                initial_state=old_configuration)
    #store the scf state
    old_configuration = molecule_configuration
    #calculate total energy
    total_energy = TotalEnergy(molecule_configuration)
    #append the total energy to energies
    energies = numpy.append(energies,total_energy.evaluate().inUnitsOf(eV))

#get the energy of atomic hydrogen
elements = [ Hydrogen]
cartesian_coordinates = [[ 0.0,0.0,0.0]]*Ang
molecule_configuration = MoleculeConfiguration(
        elements=elements,
        cartesian_coordinates=cartesian_coordinates
        )
molecule_configuration.setCalculator(calculator())
h_energy = TotalEnergy(molecule_configuration)

#calculate energies relative to atomic hydrogen
energies = energies - 2.*h_energy.evaluate().inUnitsOf(eV)

#fit polynomium to data
(a,b,c) = numpy.polyfit(distances,energies,2)

#print out the results
print()
print('distance(Angstrom) Energy (eV) polynomial_fit (eV)')
for i in range(len(distances)):
    d=distances[i]
    energy_fit = a*d*d+b*d+c
    print(' %8.2f         %8.4f         %8.4f '% (d,energies[i],energy_fit))

#calculate the bonding distance (experimental 0.742* Angstrom)
d_min = -b/(2.*a)
print('bonding distance     = %8.3f Angstrom'% d_min)

# calculate vibrational frequency (experimental 4395 cm-1)
k = 2.*a*eV/(Ang*Ang)
effective_mass = Hydrogen.atomicMass()/math.sqrt(2)
omega = (k/effective_mass)**(0.5)
wave_number = (omega/speed_of_light/2./math.pi).inUnitsOf(Meter**-1)/100.
print('vibrational frequency = %8.f cm-1 '% (wave_number))

#calculate H-H bond dissociation energy (experimental 436 kJ/mol)
e_min = (a*d_min*d_min+b*d_min+c)*eV
#add zero point energy
e_min = e_min + 0.5 * hbar * omega
#convert to kj/mol
e_min_kjmol = (e_min*avogadro_number).inUnitsOf(Joule/Units.Mol)/1000.
print('dissociation energy   = %8.f kJ/mol'% e_min_kjmol)