9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 19:43:32 +01:00

Merge branch 'dev-stable' of https://github.com/AbdAmmar/qp2 into dev-stable

This commit is contained in:
Abdallah Ammar 2024-03-07 07:35:03 +01:00
commit 0c2f912497
13 changed files with 336 additions and 248 deletions

View File

@ -224,14 +224,18 @@ def write_ezfio(res, filename):
exponent += [p.expo for p in b.prim]
ang_mom.append(str.count(s, "z"))
shell_prim_num.append(len(b.prim))
shell_index += [nshell_tot+1] * len(b.prim)
shell_index += [nshell_tot] * len(b.prim)
shell_num = len(ang_mom)
assert(shell_index[0] = 1)
assert(shell_index[-1] = shell_num)
# ~#~#~#~#~ #
# W r i t e #
# ~#~#~#~#~ #
ezfio.set_basis_basis("Read from ResultsFile")
ezfio.set_basis_shell_num(len(ang_mom))
ezfio.set_basis_shell_num(shell_num)
ezfio.set_basis_basis_nucleus_index(nucl_index)
ezfio.set_basis_prim_num(len(coefficient))

View File

@ -31,11 +31,12 @@ subroutine run_pt2_slave(thread,iproc,energy)
double precision, intent(in) :: energy(N_states_diag)
integer, intent(in) :: thread, iproc
if (N_det > 100000 ) then
call run_pt2_slave_large(thread,iproc,energy)
else
call run_pt2_slave_small(thread,iproc,energy)
endif
! if (N_det > 100000 ) then
! call run_pt2_slave_large(thread,iproc,energy)
! else
! call run_pt2_slave_small(thread,iproc,energy)
! endif
end
subroutine run_pt2_slave_small(thread,iproc,energy)
@ -178,15 +179,12 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
type(pt2_type) :: pt2_data
integer :: n_tasks, k, N
integer :: i_generator, subset
integer :: ifirst
integer :: bsize ! Size of selection buffers
logical :: sending
double precision :: time_shift
PROVIDE global_selection_buffer global_selection_buffer_lock
call random_number(time_shift)
time_shift = time_shift*15.d0
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
@ -198,15 +196,13 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
zmq_socket_push = new_zmq_push_socket(thread)
ifirst = 0
b%N = 0
buffer_ready = .False.
n_tasks = 1
sending = .False.
done = .False.
double precision :: time0, time1
call wall_time(time0)
time0 = time0+time_shift
do while (.not.done)
integer, external :: get_tasks_from_taskserver
@ -233,28 +229,29 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
ASSERT (b%N == bsize)
endif
double precision :: time0, time1
call wall_time(time0)
call pt2_alloc(pt2_data,N_states)
b%cur = 0
call select_connected(i_generator,energy,pt2_data,b,subset,pt2_F(i_generator))
call wall_time(time1)
integer, external :: tasks_done_to_taskserver
if (tasks_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,n_tasks) == -1) then
done = .true.
endif
call sort_selection_buffer(b)
call wall_time(time1)
! if (time1-time0 > 15.d0) then
call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending)
call omp_set_lock(global_selection_buffer_lock)
global_selection_buffer%mini = b%mini
call merge_selection_buffers(b,global_selection_buffer)
if (ifirst /= 0 ) then
b%cur=0
else
ifirst = 1
endif
call omp_unset_lock(global_selection_buffer_lock)
call wall_time(time0)
! endif
call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending)
if ( iproc == 1 .or. i_generator < 100 .or. done) then
if ( iproc == 1 ) then
call omp_set_lock(global_selection_buffer_lock)
call push_pt2_results_async_send(zmq_socket_push, (/i_generator/), (/pt2_data/), global_selection_buffer, (/task_id/), 1,sending)
global_selection_buffer%cur = 0

View File

@ -7,17 +7,20 @@ subroutine run_selection_slave(thread, iproc, energy)
double precision, intent(in) :: energy(N_states)
integer, intent(in) :: thread, iproc
integer :: rc, i
integer :: worker_id, task_id(1), ctask, ltask
character*(512) :: task
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_socket_push
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer(ZMQ_PTR), external :: new_zmq_push_socket
integer(ZMQ_PTR) :: zmq_socket_push
type(selection_buffer) :: buf, buf2
type(pt2_type) :: pt2_data
logical :: done, buffer_ready
type(pt2_type) :: pt2_data
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order

View File

@ -11,7 +11,7 @@ subroutine ZMQ_selection(N_in, pt2_data)
integer, external :: omp_get_thread_num
type(pt2_type), intent(inout) :: pt2_data
PROVIDE psi_det psi_coef N_det qp_max_mem N_states pt2_F s2_eig N_det_generators
! PROVIDE psi_det psi_coef N_det qp_max_mem N_states pt2_F s2_eig N_det_generators
N = max(N_in,1)
N = min(N, (elec_alpha_num * (mo_num-elec_alpha_num))**2)
@ -61,7 +61,6 @@ subroutine ZMQ_selection(N_in, pt2_data)
ipos=1
task = ' '
do i= 1, N_det_generators
do j=1,pt2_F(i)
write(task(ipos:ipos+30),'(I9,1X,I9,1X,I9,''|'')') j, i, N

View File

@ -8,8 +8,13 @@ subroutine provide_all_three_ints_bi_ortho()
END_DOC
implicit none
double precision :: t1, t2
PROVIDE ao_two_e_integrals_in_map
print *, ' start provide_all_three_ints_bi_ortho'
call wall_time(t1)
if(three_body_h_tc) then
if(three_e_3_idx_term) then
@ -32,6 +37,9 @@ subroutine provide_all_three_ints_bi_ortho()
endif
call wall_time(t2)
print *, ' end provide_all_three_ints_bi_ortho after (min) = ', (t2-t1)/60.d0
return
end
@ -83,8 +91,11 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hmono, htwoe, hthree, htot
integer :: degree
PROVIDE pure_three_body_h_tc
hmono = 0.d0
htwoe = 0.d0
htot = 0.d0

View File

@ -15,6 +15,8 @@
implicit none
double precision :: hmono, htwoe, htot, hthree
PROVIDE N_int
PROVIDE HF_bitmask
PROVIDE mo_l_coef mo_r_coef
call diag_htilde_mu_mat_bi_ortho_slow(N_int, HF_bitmask, hmono, htwoe, htot)

View File

@ -19,6 +19,7 @@ subroutine single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe,
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
double precision, intent(out) :: hmono, htwoe, hthree, htot
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
integer :: degree,exc(0:2,2,2)
@ -44,27 +45,28 @@ subroutine single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe,
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
call get_single_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc, 1, h1, p1, h2, p2, s1, s2)
call get_single_excitation_from_fock_tc(key_i, key_j, h1, p1, s1, phase, hmono, htwoe, hthree, htot)
call get_single_excitation_from_fock_tc(Nint, key_i, key_j, h1, p1, s1, phase, hmono, htwoe, hthree, htot)
end
! ---
subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, hmono, htwoe, hthree, htot)
subroutine get_single_excitation_from_fock_tc(Nint, key_i, key_j, h, p, spin, phase, hmono, htwoe, hthree, htot)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer, intent(in) :: h, p, spin
double precision, intent(in) :: phase
integer(bit_kind), intent(in) :: key_i(N_int,2), key_j(N_int,2)
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hmono, htwoe, hthree, htot
integer(bit_kind) :: differences(N_int,2)
integer(bit_kind) :: hole(N_int,2)
integer(bit_kind) :: partcl(N_int,2)
integer :: occ_hole(N_int*bit_kind_size,2)
integer :: occ_partcl(N_int*bit_kind_size,2)
integer(bit_kind) :: differences(Nint,2)
integer(bit_kind) :: hole(Nint,2)
integer(bit_kind) :: partcl(Nint,2)
integer :: occ_hole(Nint*bit_kind_size,2)
integer :: occ_partcl(Nint*bit_kind_size,2)
integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
integer :: i0,i
double precision :: buffer_c(mo_num),buffer_x(mo_num)
@ -74,7 +76,7 @@ subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, h
buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h)
enddo
do i = 1, N_int
do i = 1, Nint
differences(i,1) = xor(key_i(i,1), ref_closed_shell_bitmask(i,1))
differences(i,2) = xor(key_i(i,2), ref_closed_shell_bitmask(i,2))
hole (i,1) = iand(differences(i,1), ref_closed_shell_bitmask(i,1))
@ -83,8 +85,8 @@ subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, h
partcl (i,2) = iand(differences(i,2), key_i(i,2))
enddo
call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int)
call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int)
call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, Nint)
call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, Nint)
hmono = mo_bi_ortho_tc_one_e(p,h)
htwoe = fock_op_2_e_tc_closed_shell(p,h)
@ -122,7 +124,7 @@ subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, h
hthree = 0.d0
if (three_body_h_tc .and. elec_num.gt.2 .and. three_e_4_idx_term) then
call three_comp_fock_elem(key_i, h, p, spin, hthree)
call three_comp_fock_elem(Nint, key_i, h, p, spin, hthree)
endif
htwoe = htwoe * phase
@ -134,24 +136,27 @@ end
! ---
subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree)
implicit none
integer,intent(in) :: h_fock,p_fock,ispin_fock
integer(bit_kind), intent(in) :: key_i(N_int,2)
double precision, intent(out) :: hthree
integer :: nexc(2),i,ispin,na,nb
integer(bit_kind) :: hole(N_int,2)
integer(bit_kind) :: particle(N_int,2)
integer :: occ_hole(N_int*bit_kind_size,2)
integer :: occ_particle(N_int*bit_kind_size,2)
integer :: n_occ_ab_hole(2),n_occ_ab_particle(2)
integer(bit_kind) :: det_tmp(N_int,2)
subroutine three_comp_fock_elem(Nint, key_i, h_fock, p_fock, ispin_fock, hthree)
implicit none
integer, intent(in) :: Nint
integer, intent(in) :: h_fock, p_fock, ispin_fock
integer(bit_kind), intent(in) :: key_i(Nint,2)
double precision, intent(out) :: hthree
integer :: nexc(2),i,ispin,na,nb
integer(bit_kind) :: hole(Nint,2)
integer(bit_kind) :: particle(Nint,2)
integer :: occ_hole(Nint*bit_kind_size,2)
integer :: occ_particle(Nint*bit_kind_size,2)
integer :: n_occ_ab_hole(2),n_occ_ab_particle(2)
integer(bit_kind) :: det_tmp(Nint,2)
nexc(1) = 0
nexc(2) = 0
!! Get all the holes and particles of key_i with respect to the ROHF determinant
do i=1,N_int
do i = 1, Nint
hole(i,1) = xor(key_i(i,1),ref_bitmask(i,1))
hole(i,2) = xor(key_i(i,2),ref_bitmask(i,2))
particle(i,1) = iand(hole(i,1),key_i(i,1))
@ -161,13 +166,14 @@ subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree)
nexc(1) = nexc(1) + popcnt(hole(i,1))
nexc(2) = nexc(2) + popcnt(hole(i,2))
enddo
integer :: tmp(2)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(particle, occ_particle, tmp, N_int)
call bitstring_to_list_ab(particle, occ_particle, tmp, Nint)
ASSERT (tmp(1) == nexc(1)) ! Number of particles alpha
ASSERT (tmp(2) == nexc(2)) ! Number of particle beta
!DIR$ FORCEINLINE
call bitstring_to_list_ab(hole, occ_hole, tmp, N_int)
call bitstring_to_list_ab(hole, occ_hole, tmp, Nint)
ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha
ASSERT (tmp(2) == nexc(2)) ! Number of holes beta
@ -183,13 +189,16 @@ subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree)
nb = elec_num_tab(iand(ispin,1)+1)
do i = 1, nexc(ispin)
!DIR$ FORCEINLINE
call fock_ac_tc_operator( occ_particle(i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, N_int,na,nb)
call fock_ac_tc_operator( occ_particle(i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, Nint, na, nb)
!DIR$ FORCEINLINE
call fock_a_tc_operator ( occ_hole (i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, N_int,na,nb)
call fock_a_tc_operator ( occ_hole (i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, Nint, na, nb)
enddo
enddo
end
! ---
subroutine fock_ac_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,Nint,na,nb)
use bitmasks
implicit none
@ -365,12 +374,18 @@ subroutine fock_a_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,N
end
! ---
BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num)]
implicit none
BEGIN_DOC
! Closed-shell part of the Fock operator for the TC operator
END_DOC
implicit none
PROVIDE N_int
integer :: h0,p0,h,p,k0,k,i
integer :: n_occ_ab(2)
integer :: occ(N_int*bit_kind_size,2)
@ -382,6 +397,7 @@ BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num)
fock_op_2_e_tc_closed_shell = -1000.d0
call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int)
do i = 1, N_int
key_virt(i,1) = full_ijkl_bitmask(i)
key_virt(i,2) = full_ijkl_bitmask(i)
@ -477,8 +493,10 @@ BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num)
END_PROVIDER
! ---
subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
BEGIN_DOC
! <key_j |H_tilde | key_i> for single excitation ONLY FOR ONE- AND TWO-BODY TERMS
!!
@ -493,6 +511,7 @@ subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
double precision, intent(out) :: htot
double precision :: hmono, htwoe
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
@ -517,30 +536,38 @@ subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
call get_single_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc,1,h1,p1,h2,p2,s1,s2)
call get_single_excitation_from_fock_tc_no_3e(key_i,key_j,h1,p1,s1,phase,hmono,htwoe,htot)
call get_single_excitation_from_fock_tc_no_3e(Nint, key_i, key_j, h1, p1, s1, phase, hmono, htwoe, htot)
end
! ---
subroutine get_single_excitation_from_fock_tc_no_3e(Nint, key_i, key_j, h, p, spin, phase, hmono, htwoe, htot)
subroutine get_single_excitation_from_fock_tc_no_3e(key_i,key_j,h,p,spin,phase,hmono,htwoe,htot)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer, intent(in) :: h, p, spin
double precision, intent(in) :: phase
integer(bit_kind), intent(in) :: key_i(N_int,2), key_j(N_int,2)
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hmono,htwoe,htot
integer(bit_kind) :: differences(N_int,2)
integer(bit_kind) :: hole(N_int,2)
integer(bit_kind) :: partcl(N_int,2)
integer :: occ_hole(N_int*bit_kind_size,2)
integer :: occ_partcl(N_int*bit_kind_size,2)
integer(bit_kind) :: differences(Nint,2)
integer(bit_kind) :: hole(Nint,2)
integer(bit_kind) :: partcl(Nint,2)
integer :: occ_hole(Nint*bit_kind_size,2)
integer :: occ_partcl(Nint*bit_kind_size,2)
integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
integer :: i0,i
double precision :: buffer_c(mo_num), buffer_x(mo_num)
do i = 1, mo_num
buffer_c(i) = tc_2e_3idx_coulomb_integrals(i,p,h)
buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h)
enddo
do i = 1, N_int
do i = 1, Nint
differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1))
differences(i,2) = xor(key_i(i,2),ref_closed_shell_bitmask(i,2))
hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1))
@ -548,10 +575,12 @@ subroutine get_single_excitation_from_fock_tc_no_3e(key_i,key_j,h,p,spin,phase,h
partcl(i,1) = iand(differences(i,1),key_i(i,1))
partcl(i,2) = iand(differences(i,2),key_i(i,2))
enddo
call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int)
call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int)
call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, Nint)
call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, Nint)
hmono = mo_bi_ortho_tc_one_e(p,h)
htwoe = fock_op_2_e_tc_closed_shell(p,h)
! holes :: direct terms
do i0 = 1, n_occ_ab_hole(1)
i = occ_hole(i0,1)

View File

@ -7,6 +7,10 @@ program tc_bi_ortho
!
END_DOC
implicit none
PROVIDE N_int
my_grid_becke = .True.
PROVIDE tc_grid1_a tc_grid1_r
my_n_pt_r_grid = tc_grid1_r
@ -66,6 +70,15 @@ subroutine routine_diag()
! provide overlap_bi_ortho
! provide htilde_matrix_elmt_bi_ortho
if(noL_standard) then
PROVIDE noL_0e
PROVIDE noL_1e
PROVIDE noL_2e
endif
PROVIDE htilde_matrix_elmt_bi_ortho
return
if(N_states .eq. 1) then
print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1)

View File

@ -13,16 +13,34 @@ BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)]
implicit none
integer :: i, j
double precision :: t1, t2
double precision :: htot
call provide_all_three_ints_bi_ortho
PROVIDE N_int
PROVIDE psi_det
PROVIDE three_e_3_idx_term
if(noL_standard) then
PROVIDE noL_0e
PROVIDE noL_1e
PROVIDE noL_2e
endif
print *, ' PROVIDING htilde_matrix_elmt_bi_ortho ...'
call wall_time(t1)
call provide_all_three_ints_bi_ortho()
i = 1
j = 1
call htilde_mu_mat_opt_bi_ortho_tot(psi_det(1,1,j), psi_det(1,1,i), N_int, htot)
!$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j, htot) &
!$OMP PARALLEL &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(i, j, htot) &
!$OMP SHARED (N_det, psi_det, N_int, htilde_matrix_elmt_bi_ortho)
!$OMP DO
do i = 1, N_det
do j = 1, N_det
! < J |Htilde | I >
@ -31,7 +49,11 @@ BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)]
htilde_matrix_elmt_bi_ortho(j,i) = htot
enddo
enddo
!$OMP END PARALLEL DO
!$OMP END DO
!$OMP END PARALLEL
call wall_time(t2)
print *, ' wall time for htilde_matrix_elmt_bi_ortho (min) =', (t2-t1)/60.d0
END_PROVIDER

View File

@ -186,6 +186,7 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
type(pt2_type) :: pt2_data
integer :: n_tasks, k, N
integer :: i_generator, subset
integer :: ifirst
integer :: bsize ! Size of selection buffers
logical :: sending
@ -202,6 +203,7 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
zmq_socket_push = new_zmq_push_socket(thread)
ifirst = 0
b%N = 0
buffer_ready = .False.
n_tasks = 1
@ -250,7 +252,11 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
call omp_set_lock(global_selection_buffer_lock)
global_selection_buffer%mini = b%mini
call merge_selection_buffers(b,global_selection_buffer)
if (ifirst /= 0 ) then
b%cur=0
else
ifirst = 1
endif
call omp_unset_lock(global_selection_buffer_lock)
if ( iproc == 1 ) then
call omp_set_lock(global_selection_buffer_lock)

View File

@ -20,4 +20,5 @@ subroutine routine
call diagonalize_CI
print*,'N_det = ',N_det
call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
call print_mol_properties
end

View File

@ -14,5 +14,6 @@ end
subroutine run
implicit none
call print_mol_properties
print *, psi_energy + nuclear_repulsion
end