PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-26 15:12:19 +02:00
Code Releases Activity

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

Browse Source
This commit is contained in:
AbdAmmar 2024-03-11 10:22:07 +01:00
commit c3c7c69853
19 changed files with 497 additions and 258 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
bin/qp_convert_output_to_ezfio
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))
@ -309,10 +313,19 @@ def write_ezfio(res, filename):
MoMatrix = []
sym0 = [i.sym for i in res.mo_sets[MO_type]]
sym = [i.sym for i in res.mo_sets[MO_type]]
sym = [i.sym for i in res.mo_sets[MO_type]]
for i in range(len(sym)):
sym[MOmap[i]] = sym0[i]
irrep = {}
for i in sym:
irrep[i] = 0
for i, j in enumerate(irrep.keys()):
irrep[j] = i+1
sym = [ irrep[k] for k in sym ]
MoMatrix = []
for i in range(len(MOs)):
m = MOs[i]
@ -329,6 +342,7 @@ def write_ezfio(res, filename):
ezfio.set_mo_basis_mo_num(mo_num)
ezfio.set_mo_basis_mo_coef(MoMatrix)
ezfio.set_mo_basis_mo_occ(OccNum)
ezfio.set_mo_basis_mo_symmetry(sym)
print("OK")

45
plugins/local/cipsi_tc_bi_ortho/run_pt2_slave.irp.f
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
call run_pt2_slave_large(thread,iproc,energy)
! 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 omp_set_lock(global_selection_buffer_lock)
global_selection_buffer%mini = b%mini
call merge_selection_buffers(b,global_selection_buffer)
b%cur=0
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
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)
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

29
plugins/local/cipsi_tc_bi_ortho/run_selection_slave.irp.f
View File

@ -5,19 +5,22 @@ subroutine run_selection_slave(thread, iproc, energy)
implicit none
double precision, intent(in) :: energy(N_states)
integer, intent(in) :: thread, iproc
double precision, intent(in) :: energy(N_states)
integer, intent(in) :: thread, iproc
integer :: rc, i
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), external :: new_zmq_push_socket
type(selection_buffer) :: buf, buf2
type(pt2_type) :: pt2_data
logical :: done, buffer_ready
integer :: worker_id, task_id(1), ctask, ltask
character*(512) :: task
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
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
@ -64,7 +67,7 @@ subroutine run_selection_slave(thread, iproc, energy)
stop '-1'
end if
end if
call select_connected(i_generator, energy, pt2_data, buf,subset, pt2_F(i_generator))
call select_connected(i_generator, energy, pt2_data, buf, subset, pt2_F(i_generator))
endif
integer, external :: task_done_to_taskserver

3
plugins/local/cipsi_tc_bi_ortho/zmq_selection.irp.f
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

23
plugins/local/tc_bi_ortho/slater_tc_opt.irp.f
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
@ -99,7 +110,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
if(degree == 0) then
call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot)
else if (degree == 1) then
call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i , hmono, htwoe, hthree, htot)
call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
else if(degree == 2) then
call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
endif
@ -111,7 +122,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
if(degree == 0) then
call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot)
else if (degree == 1) then
call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i , hmono, htwoe, hthree, htot)
call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
else if(degree == 2) then
call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
else
@ -149,16 +160,16 @@ subroutine htilde_mu_mat_opt_bi_ortho_no_3e(key_j, key_i, Nint, htot)
double precision, intent(out) :: htot
integer :: degree
htot = 0.d0
htot = 0.d0
call get_excitation_degree(key_i, key_j, degree, Nint)
if(degree.gt.2) return
if(degree == 0)then
if(degree == 0) then
call diag_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_i,htot)
else if (degree == 1)then
else if (degree == 1) then
call single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint,key_j, key_i , htot)
else if(degree == 2)then
else if(degree == 2) then
call double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
endif

2
plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f
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)

419
plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f
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
@ -181,15 +187,18 @@ subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree)
do ispin=1,2
na = elec_num_tab(ispin)
nb = elec_num_tab(iand(ispin,1)+1)
do i=1,nexc(ispin)
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,111 +374,118 @@ 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
integer :: h0,p0,h,p,k0,k,i
integer :: n_occ_ab(2)
integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab_virt(2)
integer :: occ_virt(N_int*bit_kind_size,2)
integer(bit_kind) :: key_test(N_int)
integer(bit_kind) :: key_virt(N_int,2)
double precision :: accu
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)
key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1))
key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2))
enddo
call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int)
! docc ---> virt single excitations
do h0 = 1, n_occ_ab(1)
h=occ(h0,1)
do p0 = 1, n_occ_ab_virt(1)
p = occ_virt(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
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)
integer :: n_occ_ab_virt(2)
integer :: occ_virt(N_int*bit_kind_size,2)
integer(bit_kind) :: key_test(N_int)
integer(bit_kind) :: key_virt(N_int,2)
double precision :: accu
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)
key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1))
key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2))
enddo
enddo
do h0 = 1, n_occ_ab_virt(1)
h = occ_virt(h0,1)
do p0 = 1, n_occ_ab(1)
p=occ(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int)
! docc ---> virt single excitations
do h0 = 1, n_occ_ab(1)
h = occ(h0,1)
do p0 = 1, n_occ_ab_virt(1)
p = occ_virt(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
enddo
enddo
! virt ---> virt single excitations
do h0 = 1, n_occ_ab_virt(1)
h=occ_virt(h0,1)
do p0 = 1, n_occ_ab_virt(1)
p = occ_virt(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
do h0 = 1, n_occ_ab_virt(1)
h = occ_virt(h0,1)
do p0 = 1, n_occ_ab(1)
p = occ(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
enddo
enddo
do h0 = 1, n_occ_ab_virt(1)
h = occ_virt(h0,1)
do p0 = 1, n_occ_ab_virt(1)
p=occ_virt(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
! virt ---> virt single excitations
do h0 = 1, n_occ_ab_virt(1)
h=occ_virt(h0,1)
do p0 = 1, n_occ_ab_virt(1)
p = occ_virt(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
enddo
! docc ---> docc single excitations
do h0 = 1, n_occ_ab(1)
h=occ(h0,1)
do p0 = 1, n_occ_ab(1)
p = occ(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
do h0 = 1, n_occ_ab_virt(1)
h = occ_virt(h0,1)
do p0 = 1, n_occ_ab_virt(1)
p=occ_virt(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
enddo
do h0 = 1, n_occ_ab(1)
h = occ(h0,1)
do p0 = 1, n_occ_ab(1)
p=occ(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
! docc ---> docc single excitations
do h0 = 1, n_occ_ab(1)
h=occ(h0,1)
do p0 = 1, n_occ_ab(1)
p = occ(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
enddo
do h0 = 1, n_occ_ab(1)
h = occ(h0,1)
do p0 = 1, n_occ_ab(1)
p=occ(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
enddo
! do i = 1, mo_num
! write(*,'(100(F10.5,X))')fock_op_2_e_tc_closed_shell(:,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
!!
@ -492,8 +510,9 @@ subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
implicit none
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
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
integer :: degree,exc(0:2,2,2)
@ -517,75 +536,85 @@ 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)
end
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) :: h,p,spin
double precision, intent(in) :: phase
integer(bit_kind), intent(in) :: key_i(N_int,2), key_j(N_int,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 :: 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
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))
hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2))
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)
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)
htwoe -= buffer_c(i)
enddo
do i0 = 1, n_occ_ab_hole(2)
i = occ_hole(i0,2)
htwoe -= buffer_c(i)
enddo
! holes :: exchange terms
do i0 = 1, n_occ_ab_hole(spin)
i = occ_hole(i0,spin)
htwoe += buffer_x(i)
enddo
! particles :: direct terms
do i0 = 1, n_occ_ab_partcl(1)
i = occ_partcl(i0,1)
htwoe += buffer_c(i)
enddo
do i0 = 1, n_occ_ab_partcl(2)
i = occ_partcl(i0,2)
htwoe += buffer_c(i)
enddo
! particles :: exchange terms
do i0 = 1, n_occ_ab_partcl(spin)
i = occ_partcl(i0,spin)
htwoe -= buffer_x(i)
enddo
htwoe = htwoe * phase
hmono = hmono * phase
htot = htwoe + hmono
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)
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(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hmono,htwoe,htot
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, 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))
hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2))
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, 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)
htwoe -= buffer_c(i)
enddo
do i0 = 1, n_occ_ab_hole(2)
i = occ_hole(i0,2)
htwoe -= buffer_c(i)
enddo
! holes :: exchange terms
do i0 = 1, n_occ_ab_hole(spin)
i = occ_hole(i0,spin)
htwoe += buffer_x(i)
enddo
! particles :: direct terms
do i0 = 1, n_occ_ab_partcl(1)
i = occ_partcl(i0,1)
htwoe += buffer_c(i)
enddo
do i0 = 1, n_occ_ab_partcl(2)
i = occ_partcl(i0,2)
htwoe += buffer_c(i)
enddo
! particles :: exchange terms
do i0 = 1, n_occ_ab_partcl(spin)
i = occ_partcl(i0,spin)
htwoe -= buffer_x(i)
enddo
htwoe = htwoe * phase
hmono = hmono * phase
htot = htwoe + hmono
end

13
plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f
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)

28
plugins/local/tc_bi_ortho/tc_hmat.irp.f
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
PROVIDE N_int
PROVIDE psi_det
PROVIDE three_e_3_idx_term
call provide_all_three_ints_bi_ortho
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

70
plugins/local/tc_scf/routines_rotates.irp.f
View File

@ -103,7 +103,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
double precision, allocatable :: stmp(:,:), T(:,:), Snew(:,:), smat2(:,:)
double precision, allocatable :: mo_l_coef_tmp(:,:), mo_r_coef_tmp(:,:), mo_l_coef_new(:,:)
E_thr = 1d-8
E_thr = 1d-04
E_old = TC_HF_energy
allocate(mo_l_coef_old(ao_num,mo_num), mo_r_coef_old(ao_num,mo_num))
mo_r_coef_old = mo_r_coef
@ -164,10 +164,42 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
allocate(mo_r_coef_tmp(ao_num,n_degen), mo_l_coef_tmp(ao_num,n_degen), mo_l_coef_new(ao_num,n_degen))
allocate(T(n_degen,n_degen), Snew(n_degen,n_degen))
do j = 1, n_degen
mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j))
mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j))
enddo
print*,'Right orbitals before'
do j = 1, n_degen
write(*,'(100(F16.10,X))') mo_r_coef_new(1:ao_num,list_degen(i,j))
enddo
print*,'Left orbitals before'
do j = 1, n_degen
write(*,'(100(F16.10,X))')mo_l_coef(1:ao_num,list_degen(i,j))
enddo
if(angle_left_right(list_degen(i,1)).gt.80.d0.and.n_degen==2)then
integer :: i_list, j_list
i_list = list_degen(i,1)
j_list = list_degen(i,2)
print*,'Huge angle !!! == ',angle_left_right(list_degen(i,1)),angle_left_right(list_degen(i,2))
print*,'i_list = ',i_list
print*,'i_list = ',j_list
print*,'Swapping left/right orbitals'
call print_strong_overlap(i_list, j_list)
mo_r_coef_tmp(1:ao_num,1) = mo_r_coef_new(1:ao_num,i_list)
mo_r_coef_tmp(1:ao_num,2) = mo_l_coef(1:ao_num,i_list)
mo_l_coef_tmp(1:ao_num,1) = mo_l_coef(1:ao_num,j_list)
mo_l_coef_tmp(1:ao_num,2) = mo_r_coef_new(1:ao_num,j_list)
else
do j = 1, n_degen
print*,'i_list = ',list_degen(i,j)
mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j))
mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j))
enddo
endif
print*,'Right orbitals '
do j = 1, n_degen
write(*,'(100(F16.10,X))')mo_r_coef_tmp(1:ao_num,j)
enddo
print*,'Left orbitals '
do j = 1, n_degen
write(*,'(100(F16.10,X))')mo_l_coef_tmp(1:ao_num,j)
enddo
! Orthogonalization of right functions
print *, ' Orthogonalization of RIGHT functions'
print *, ' ------------------------------------'
@ -445,3 +477,31 @@ subroutine sort_by_tc_fock
end
subroutine print_strong_overlap(i_list, j_list)
implicit none
integer, intent(in) :: i_list,j_list
double precision :: o_i, o_j,o_ij
double precision :: s_mat_r(2,2),s_mat_l(2,2)
o_i = dsqrt(overlap_mo_r(i_list, i_list))
o_j = dsqrt(overlap_mo_r(j_list, j_list))
o_ij = overlap_mo_r(j_list, i_list)
s_mat_r(1,1) = o_i*o_i
s_mat_r(2,1) = o_ij/(o_i * o_j)
s_mat_r(2,2) = o_j*o_j
s_mat_r(1,2) = s_mat_r(2,1)
print*,'Right overlap matrix '
write(*,'(2(F10.5,X))')s_mat_r(1:2,1)
write(*,'(2(F10.5,X))')s_mat_r(1:2,2)
o_i = dsqrt(overlap_mo_l(i_list, i_list))
o_j = dsqrt(overlap_mo_l(j_list, j_list))
o_ij = overlap_mo_l(j_list, i_list)
s_mat_l(1,1) = o_i*o_i
s_mat_l(2,1) = o_ij/(o_i * o_j)
s_mat_l(2,2) = o_j*o_j
s_mat_l(1,2) = s_mat_l(2,1)
print*,'Left overlap matrix '
write(*,'(2(F10.5,X))')s_mat_l(1:2,1)
write(*,'(2(F10.5,X))')s_mat_l(1:2,2)
end

8
src/cipsi/run_pt2_slave.irp.f
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)
b%cur=0
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)

2
src/cipsi/run_selection_slave.irp.f
View File

@ -65,7 +65,7 @@ subroutine run_selection_slave(thread,iproc,energy)
stop '-1'
end if
end if
call select_connected(i_generator,energy,pt2_data,buf,subset,pt2_F(i_generator))
call select_connected(i_generator, energy, pt2_data, buf, subset, pt2_F(i_generator))
endif
integer, external :: task_done_to_taskserver

6
src/mo_basis/EZFIO.cfg
View File

@ -32,6 +32,12 @@ doc: |MO| occupation numbers
interface: ezfio
size: (mo_basis.mo_num)
[mo_symmetry]
type: integer
doc: MOs with the same integer belong to the same irrep.
interface: ezfio
size: (mo_basis.mo_num)
[mo_class]
type: MO_class
doc: [ Core | Inactive | Active | Virtual | Deleted ], as defined by :ref:`qp_set_mo_class`

18
src/mo_one_e_ints/spread_dipole_mo.irp.f
View File

@ -58,3 +58,21 @@ END_PROVIDER
)
END_PROVIDER
BEGIN_PROVIDER [double precision, mo_spread_centered_x, (mo_num, mo_num) ]
&BEGIN_PROVIDER [double precision, mo_spread_centered_y, (mo_num, mo_num) ]
&BEGIN_PROVIDER [double precision, mo_spread_centered_z, (mo_num, mo_num) ]
BEGIN_DOC
! array of the integrals of MO_i * (x^2 - <MO_i|x|MO_j>^2) MO_j = MO_i x^2 MO_j - (MO_i x MO_j)^2
! array of the integrals of MO_i * (y^2 - <MO_i|y|MO_j>^2) MO_j = MO_i y^2 MO_j - (MO_i y MO_j)^2
! array of the integrals of MO_i * (z^2 - <MO_i|z|MO_j>^2) MO_j = MO_i z^2 MO_j - (MO_i z MO_j)^2
END_DOC
implicit none
integer :: i,j
do i = 1, mo_num
do j = 1, mo_num
mo_spread_centered_x(j,i) = mo_spread_x(j,i) - mo_dipole_x(j,i)**2
mo_spread_centered_y(j,i) = mo_spread_y(j,i) - mo_dipole_y(j,i)**2
mo_spread_centered_z(j,i) = mo_spread_z(j,i) - mo_dipole_z(j,i)**2
enddo
enddo
END_PROVIDER

23
src/mol_properties/multi_s_dipole_moment.irp.f
View File

@ -91,3 +91,26 @@ BEGIN_PROVIDER [double precision, multi_s_dipole_moment, (N_states, N_states)]
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, multi_s_x_dipole_moment_eigenvec, (N_states, N_states)]
&BEGIN_PROVIDER [double precision, multi_s_y_dipole_moment_eigenvec, (N_states, N_states)]
&BEGIN_PROVIDER [double precision, multi_s_z_dipole_moment_eigenvec, (N_states, N_states)]
&BEGIN_PROVIDER [double precision, multi_s_x_dipole_moment_eigenval, (N_states)]
&BEGIN_PROVIDER [double precision, multi_s_y_dipole_moment_eigenval, (N_states)]
&BEGIN_PROVIDER [double precision, multi_s_z_dipole_moment_eigenval, (N_states)]
implicit none
PROVIDE multi_s_x_dipole_moment multi_s_y_dipole_moment multi_s_z_dipole_moment
call lapack_diag(multi_s_x_dipole_moment_eigenval(1), multi_s_x_dipole_moment_eigenvec(1,1), multi_s_x_dipole_moment(1,1), N_states, N_states)
call lapack_diag(multi_s_y_dipole_moment_eigenval(1), multi_s_y_dipole_moment_eigenvec(1,1), multi_s_y_dipole_moment(1,1), N_states, N_states)
call lapack_diag(multi_s_z_dipole_moment_eigenval(1), multi_s_z_dipole_moment_eigenvec(1,1), multi_s_z_dipole_moment(1,1), N_states, N_states)
END_PROVIDER
! ---

1
src/tools/diagonalize_h.irp.f
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

3
src/tools/print_energy.irp.f
View File

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

27
src/utils/constants.include.F
View File

@ -18,3 +18,30 @@ double precision, parameter :: c_4_3 = 4.d0/3.d0
double precision, parameter :: c_1_3 = 1.d0/3.d0
double precision, parameter :: sq_op5 = dsqrt(0.5d0)
double precision, parameter :: dlog_2pi = dlog(2.d0*dacos(-1.d0))
! physical constants and units conversion factors
double precision, parameter :: k_boltzman_si = 1.38066d-23 ! K k^-1
double precision, parameter :: k_boltzman_au = 3.1667d-6 ! Hartree k^-1
double precision, parameter :: k_boltzman_m1_au = 315795.26d0 ! Hartree^-1 k
double precision, parameter :: bohr_radius_si = 0.529177d-10 ! m
double precision, parameter :: bohr_radius_cm = 0.529177d-8 ! cm
double precision, parameter :: bohr_radius_angs = 0.529177d0 ! Angstrom
double precision, parameter :: electronmass_si = 9.10953d-31 ! Kg
double precision, parameter :: electronmass_uma = 5.4858d-4 ! uma
double precision, parameter :: electronvolt_si = 1.6021892d-19 ! J
double precision, parameter :: uma_si = 1.66057d-27 ! Kg
double precision, parameter :: debye_si = 3.33564d-30 ! coulomb meter
double precision, parameter :: debye_au = 0.393427228d0 ! e * Bohr
double precision, parameter :: angstrom_to_au = 1.889727d0 ! au
double precision, parameter :: au_to_ohmcmm1 = 46000.0d0 ! (ohm cm)^-1
double precision, parameter :: au_to_kb = 294210.0d0 ! kbar
double precision, parameter :: au_to_eV = 27.211652d0
double precision, parameter :: uma_to_au = 1822.89d0
double precision, parameter :: au_to_terahertz = 2.4189d-5
double precision, parameter :: au_to_sec = 2.4189d-17
double precision, parameter :: au_to_fsec = 2.4189d-2
double precision, parameter :: Wcm2 = 3.5d16
double precision, parameter :: amconv = 1.66042d-24/9.1095d-28*0.5d0 ! mass conversion: a.m.u to a.u. (ry)
double precision, parameter :: uakbar = 147105.d0 ! pressure conversion from ry/(a.u)^3 to k

15
src/utils/linear_algebra.irp.f
View File

@ -1920,8 +1920,12 @@ subroutine exp_matrix(X,n,exp_X)
call get_A_squared(X,n,r2_mat)
call lapack_diagd(eigvalues,eigvectors,r2_mat,n,n)
eigvalues=-eigvalues
do i = 1,n
! t = dsqrt(t^2) where t^2 are eigenvalues of X^2
eigvalues(i) = dsqrt(eigvalues(i))
enddo
if(.False.)then
if(.false.)then
!!! For debugging and following the book intermediate
! rebuilding the matrix : X^2 = -W t^2 W^T as in 3.1.30
! matrix_tmp1 = W t^2
@ -1932,14 +1936,16 @@ subroutine exp_matrix(X,n,exp_X)
enddo
eigvalues_mat=0.d0
do i = 1,n
! t = dsqrt(t^2) where t^2 are eigenvalues of X^2
eigvalues(i) = dsqrt(eigvalues(i))
eigvalues_mat(i,i) = eigvalues(i)*eigvalues(i)
enddo
call dgemm('N','N',n,n,n,1.d0,eigvectors,size(eigvectors,1), &
eigvalues_mat,size(eigvalues_mat,1),0.d0,matrix_tmp1,size(matrix_tmp1,1))
call dgemm('N','T',n,n,n,-1.d0,matrix_tmp1,size(matrix_tmp1,1), &
eigvectors,size(eigvectors,1),0.d0,matrix_tmp2,size(matrix_tmp2,1))
print*,'r2_mat = '
do i = 1, n
write(*,'(100(F16.10,X))')r2_mat(:,i)
enddo
print*,'r2_mat new = '
do i = 1, n
write(*,'(100(F16.10,X))')matrix_tmp2(:,i)
@ -1964,7 +1970,8 @@ subroutine exp_matrix(X,n,exp_X)
if(dabs(eigvalues(i)).gt.1.d-4)then
eigvalues_mat(i,i) = dsin(eigvalues(i))/eigvalues(i)
else ! Taylor development of sin(x)/x near x=0 = 1 - x^2/6
eigvalues_mat(i,i) = 1.d0 - eigvalues(i)*eigvalues(i)*c_1_3*0.5d0
eigvalues_mat(i,i) = 1.d0 - eigvalues(i)*eigvalues(i)*c_1_3*0.5d0 &
+ eigvalues(i)*eigvalues(i)*eigvalues(i)*eigvalues(i)*c_1_3*0.025d0
endif
enddo
! matrix_tmp1 = W t^-1 sin(t)