main.py

import myfuncs as mf
import numpy as np


# Setting up the calculation

input_path = 'HF_inputs_outputs_debugging'
input_name = 'li4h4.input'

out_f = open(input_name.split('.')[0]+'.out','w', encoding='utf-8')
debug_f = open(input_name.split('.')[0]+'.debug','w', encoding='utf-8')


mf.header(out_f,False)
mf.header(debug_f,True)

# Reading the input

if 'extended' not in input_name:
  out_f.write('This is a normal input, not extended. Electronic Integrals will be computed.')
  coordinates, coeff_bs_funct, n_basis_funct, N_elec, potential_NN, atom_number,  atom_type, n_primitives, relation_atom_bf, *_ = mf.read_INPUT(input_path,input_name,out_f) # *_ means i dont care about the other terms of the return, which are necessary elsewhere.
  
  # Computing electronic integrals
  overlap_np, kinetic_np, N_attraction_np, overlap_der_np, kinetic_der_np, N_attraction_der_np = mf.GEN_one_electron_integrals(coordinates, coeff_bs_funct, n_basis_funct, atom_number, relation_atom_bf, atom_type, out_f)
  Two_e_np, Two_e_der_np = mf.GEN_two_electron_integrals(coordinates, coeff_bs_funct, n_basis_funct, atom_number, relation_atom_bf, out_f)

else:
  out_f.write('This is an extended input. We are using precomputed Electronic Integrals.')
  coordinates, coeff_bs_funct, n_basis_funct, N_elec, potential_NN, atom_number,  atom_type,  overlap_np, kinetic_np, N_attraction_np, Two_e_np, overlap_der_np, kinetic_der_np, N_attraction_der_np, Two_e_der_np = mf.read_EXTENDED_INPUT(input_path,input_name,out_f,debug_f)

if 'extended' not in input_name:
  out_f.write('\n=====================\n')
  out_f.write('~~~ SCF procedure ~~~')
  out_f.write('\n=====================')

  debug_f.write('\n=====================\n')
  debug_f.write('~~~ SCF procedure ~~~')
  debug_f.write('\n=====================')
else:
  out_f.write('\n=================================================\n')
  out_f.write('~~~ SCF procedure using precomputed integrals ~~~')
  out_f.write('\n================================================')

  debug_f.write('\n================================================\n')
  debug_f.write('~~~ SCF procedure using precomputed integrals ~~~')
  debug_f.write('\n=================================================')


overlap_inv_sq_np = mf.DIAG_overlap_matrix(overlap_np)   #before eq 21
debug_f.write('\n\n~~~~ INITIAL CALCULATIONS ~~~~\n\n')
debug_f.write('~ Inverse square root of Overlap Matrix S^-1/2 ~\n\n')
mf.print_matrix_2d(overlap_inv_sq_np,debug_f)

Hcore_np = mf.GEN_Hcore_matrix(n_basis_funct,kinetic_np,N_attraction_np)
debug_f.write('\n~ H core Matrix ~\n\n')
mf.print_matrix_2d(Hcore_np,debug_f)


#### Preparing fo first iteration
iter = 0
F_np = Hcore_np

debug_f.write('\n~~ Solving Fock Equations ~~\n')
C_np, epsilon, F_prime_np = mf.solve_Fock_eqs(F_np, overlap_inv_sq_np,iter,debug_f)

out_f.write(f"\n\n~~ INITIAL CALCULATIONS ~~ \n")
out_f.write(f"\n~ Initial Hcore matrix in orthogonal basis (F') ~ \n\n")
mf.print_matrix_2d(F_prime_np,out_f)

out_f.write(f"\n~ Initial Coefficients ~\n\n")
mf.print_matrix_2d(C_np,out_f)

# Generating density matrix
iter += 1
density_np = mf.GEN_density_np(n_basis_funct, N_elec,C_np,iter,debug_f)

out_f.write('\n~ INITIAL DENSITY MATRIX ~\n\n')
mf.print_matrix_2d(density_np, out_f)

debug_f.write('\n~~~~ STARTING SCF PROCEDURE ~~~~\n')
out_f.write('\n~~~~ STARTING SCF PROCEDURE ~~~~\n')
convergence = 1 # So the SCF procedure can start.
while True:
  if abs(convergence) < 1e-6:
    break
  
  debug_f.write(f'\n~~ SCF CYCLE {iter} ~~')
  out_f.write(f'\n~~ SCF CYCLE {iter} ~~')

  F_np = mf.GEN_new_Fock_M(n_basis_funct,Hcore_np,density_np,Two_e_np,iter,debug_f)

  C_np, epsilon, _ = mf.solve_Fock_eqs(F_np, overlap_inv_sq_np,iter,debug_f)

  E_elec = mf.GEN_electronic_energy(n_basis_funct,density_np,Hcore_np,F_np)
  E_tot = E_elec + potential_NN

  debug_f.write(f'\nElectronic Energy = {E_elec} Eh\n')
  out_f.write(f'\nElectronic Energy = {E_elec} Eh\n')

  debug_f.write(f'Total Energy = {E_tot} Eh\n')
  out_f.write(f'Total Energy = {E_tot} Eh\n')

  iter += 1
  old_density_np = density_np
  density_np = mf.GEN_density_np(n_basis_funct, N_elec,C_np,iter,debug_f)

  #convergence check
  convergence = 0
  for mu in range(n_basis_funct):
    for nu in range(n_basis_funct):
      diff_P = old_density_np[mu,nu]-density_np[mu,nu]
      if abs(diff_P) > convergence:
        convergence = diff_P

  out_f.write(f'There is a convergence of {convergence}\n\n')
  debug_f.write(f'There is a convergence of {convergence}\n')

out_f.write('\n~~~~ SCF FINAL RESULTS ~~~~ \n')
debug_f.write('\n~~~~ SCF FINAL RESULTS ~~~~ \n')

out_f.write('\n~ FINAL DENSITY MATRIX ~\n\n')
mf.print_matrix_2d(density_np, out_f)
debug_f.write('\n~ FINAL DENSITY MATRIX ~\n\n')
mf.print_matrix_2d(density_np, debug_f)

out_f.write('\n~ FINAL FOCK MATRIX ~\n\n')
mf.print_matrix_2d(F_np, out_f)
debug_f.write('\n~ FINAL FOCK MATRIX ~\n\n')
mf.print_matrix_2d(F_np, debug_f)

out_f.write('\n~ FINAL COEFFICIENT MATRIX ~\n\n')
mf.print_matrix_2d(C_np, out_f)
debug_f.write('\n~ FINAL COEFFICIENT MATRIX ~\n\n')
mf.print_matrix_2d(C_np, debug_f)

out_f.write('\n~ FINAL EPSILON MATRIX ~\n\n')
out_f.write(str(epsilon))
debug_f.write('\n~ FINAL EPSILON MATRIX ~\n\n')
debug_f.write(str(epsilon))


debug_f.write(f'\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n')
out_f.write(f'\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n')

debug_f.write(f'\nFINAL ELECTRONIC ENERGY = {E_elec} Eh\n')
out_f.write(f'\nFINAL ELECTRONIC ENERGY = {E_elec} Eh\n')

debug_f.write(f'\nFINAL TOTAL ENERGY     = {E_tot} Eh\n')
out_f.write(f'\nFINAL TOTAL ENERGY      = {E_tot} Eh\n')

debug_f.write(f'\nFINAL CONVERGENCE      = {convergence} \n')
out_f.write(f'\nFINAL CONVERGENCE       = {convergence} \n')

debug_f.write(f'\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n')
out_f.write(f'\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n')


# GRADIENT with Pre-Computed Integrals



if 'extended' not in input_name:
  out_f.write('\n==========================\n')
  out_f.write('~~~ Computing Gradient ~~~')
  out_f.write('\n==========================\n')

  debug_f.write('\n==========================\n')
  debug_f.write('~~~ Computing Gradient ~~~')
  debug_f.write('\n==========================\n')
else:
  out_f.write('\n======================================================\n')
  out_f.write('~~~ Computing Gradient using precomputed integrals ~~~')
  out_f.write('\n======================================================\n')

  debug_f.write('\n======================================================\n')
  debug_f.write('~~~ Computing Gradient using precomputed integrals ~~~')
  debug_f.write('\n======================================================\n')


grad_tot_np = mf.GEN_gradient(n_basis_funct, N_elec, epsilon, C_np, atom_number,density_np, overlap_der_np, kinetic_der_np, N_attraction_der_np, Two_e_der_np, coordinates, atom_type, out_f, debug_f)


out_f.close()
debug_f.close()