diff --git a/src/bi_ort_ints/semi_num_ints_mo.irp.f b/src/bi_ort_ints/semi_num_ints_mo.irp.f
index 0d727785..771d3274 100644
--- a/src/bi_ort_ints/semi_num_ints_mo.irp.f
+++ b/src/bi_ort_ints/semi_num_ints_mo.irp.f
@@ -138,10 +138,13 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3,
enddo
enddo
+ FREE int2_grad1_u12_ao
+
endif
call wall_time(wall1)
print *, ' wall time for int2_grad1_u12_ao_transp ', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
@@ -200,6 +203,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid,
enddo
enddo
+ FREE int2_grad1_u12_bimo_transp
+
END_PROVIDER
! ---
diff --git a/src/bi_ort_ints/three_body_ijm.irp.f b/src/bi_ort_ints/three_body_ijm.irp.f
index 4d21cb93..b34638b8 100644
--- a/src/bi_ort_ints/three_body_ijm.irp.f
+++ b/src/bi_ort_ints/three_body_ijm.irp.f
@@ -49,6 +49,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_direct_bi_ort, (mo_num, mo_num,
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_direct_bi_ort', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
@@ -102,6 +103,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_1_bi_ort, (mo_num, mo_num
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_cycle_1_bi_ort', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
@@ -155,6 +157,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_2_bi_ort, (mo_num, mo_num
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_cycle_2_bi_ort', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
@@ -208,6 +211,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch23_bi_ort, (mo_num, mo_num,
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_exch23_bi_ort', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
@@ -261,6 +265,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch13_bi_ort, (mo_num, mo_num,
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_exch13_bi_ort', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
@@ -306,6 +311,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch12_bi_ort, (mo_num, mo_num,
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_exch12_bi_ort', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
@@ -359,6 +365,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch12_bi_ort_new, (mo_num, mo_
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_exch12_bi_ort_new', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
diff --git a/src/bi_ort_ints/three_body_ijmk.irp.f b/src/bi_ort_ints/three_body_ijmk.irp.f
index 5afd49ab..95b57e37 100644
--- a/src/bi_ort_ints/three_body_ijmk.irp.f
+++ b/src/bi_ort_ints/three_body_ijmk.irp.f
@@ -43,6 +43,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_direct_bi_ort, (mo_num, mo_num,
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_direct_bi_ort', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
@@ -90,6 +91,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_1_bi_ort, (mo_num, mo_num
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_cycle_1_bi_ort', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
@@ -137,6 +139,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_2_bi_ort, (mo_num, mo_num
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_cycle_2_bi_ort', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
@@ -184,6 +187,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch23_bi_ort, (mo_num, mo_num,
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_exch23_bi_ort', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
@@ -230,6 +234,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch13_bi_ort, (mo_num, mo_num,
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_exch13_bi_ort', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
@@ -277,6 +282,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch12_bi_ort, (mo_num, mo_num,
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_exch12_bi_ort', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
diff --git a/src/bi_ort_ints/three_body_ijmkl.irp.f b/src/bi_ort_ints/three_body_ijmkl.irp.f
index bd669163..1b2e777e 100644
--- a/src/bi_ort_ints/three_body_ijmkl.irp.f
+++ b/src/bi_ort_ints/three_body_ijmkl.irp.f
@@ -58,9 +58,11 @@
enddo
enddo
enddo
+
!$OMP END DO
!$OMP END PARALLEL
+
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, n_points_final_grid, 1.d0, &
orb_mat, n_points_final_grid, &
grad_mli, n_points_final_grid, 0.d0, &
@@ -81,6 +83,7 @@
deallocate(orb_mat,grad_mli)
+
allocate(lm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
allocate(rm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
allocate(rk_grad_im(n_points_final_grid,3,mo_num,mo_num))
@@ -132,6 +135,7 @@
enddo
!$OMP END PARALLEL DO
+
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
lm_grad_ik, 3*n_points_final_grid, &
rk_grad_im, 3*n_points_final_grid, 0.d0, &
@@ -152,6 +156,7 @@
enddo
!$OMP END PARALLEL DO
+
deallocate(lm_grad_ik)
allocate(lk_grad_mi(n_points_final_grid,3,mo_num,mo_num))
@@ -198,6 +203,7 @@
enddo
!$OMP END PARALLEL DO
+
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
lk_grad_mi, 3*n_points_final_grid, &
rk_grad_im, 3*n_points_final_grid, 0.d0, &
@@ -223,8 +229,10 @@
deallocate(rk_grad_im)
enddo
+
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_bi_ort', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
diff --git a/src/bi_ort_ints/three_body_ints_bi_ort.irp.f b/src/bi_ort_ints/three_body_ints_bi_ort.irp.f
index 1962c8d6..d8145c3e 100644
--- a/src/bi_ort_ints/three_body_ints_bi_ort.irp.f
+++ b/src/bi_ort_ints/three_body_ints_bi_ort.irp.f
@@ -57,6 +57,7 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n
call wall_time(wall1)
print *, ' wall time for three_body_ints_bi_ort', wall1 - wall0
+ call print_memory_usage()
! if(write_three_body_ints_bi_ort)then
! print*,'Writing three_body_ints_bi_ort on disk ...'
! call write_array_6_index_tensor(mo_num,three_body_ints_bi_ort,name_file)
diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f
index 3f1a9bf5..44a6ae65 100644
--- a/src/non_h_ints_mu/grad_squared.irp.f
+++ b/src/non_h_ints_mu/grad_squared.irp.f
@@ -231,6 +231,7 @@ BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_g
call wall_time(time0)
PROVIDE j1b_type
+ PROVIDE int2_grad1u2_grad2u2_j1b2
do ipoint = 1, n_points_final_grid
tmp1 = v_1b(ipoint)
@@ -242,6 +243,8 @@ BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_g
enddo
enddo
+ FREE int2_grad1u2_grad2u2_j1b2
+
!if(j1b_type .eq. 0) then
! grad12_j12 = 0.d0
! do ipoint = 1, n_points_final_grid
@@ -262,6 +265,7 @@ BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_g
call wall_time(time1)
print*, ' Wall time for grad12_j12 = ', time1 - time0
+ call print_memory_usage()
END_PROVIDER
@@ -278,6 +282,9 @@ BEGIN_PROVIDER [double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_g
print*, ' providing u12sq_j1bsq ...'
call wall_time(time0)
+ ! do not free here
+ PROVIDE int2_u2_j1b2
+
do ipoint = 1, n_points_final_grid
tmp_x = v_1b_grad(1,ipoint)
tmp_y = v_1b_grad(2,ipoint)
@@ -292,6 +299,7 @@ BEGIN_PROVIDER [double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_g
call wall_time(time1)
print*, ' Wall time for u12sq_j1bsq = ', time1 - time0
+ call print_memory_usage()
END_PROVIDER
@@ -310,6 +318,9 @@ BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num,
print*, ' providing u12_grad1_u12_j1b_grad1_j1b ...'
call wall_time(time0)
+ PROVIDE int2_u_grad1u_j1b2
+ PROVIDE int2_u_grad1u_x_j1b2
+
do ipoint = 1, n_points_final_grid
x = final_grid_points(1,ipoint)
@@ -340,14 +351,17 @@ BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num,
enddo
enddo
+ FREE int2_u_grad1u_j1b2
+ FREE int2_u_grad1u_x_j1b2
+
call wall_time(time1)
print*, ' Wall time for u12_grad1_u12_j1b_grad1_j1b = ', time1 - time0
+ call print_memory_usage()
END_PROVIDER
! ---
-
BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao_num)]
BEGIN_DOC
@@ -401,6 +415,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
, int2_grad1_u12_square_ao(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
, 0.d0, tc_grad_square_ao, ao_num*ao_num)
+ FREE int2_grad1_u12_square_ao
+
! ---
if(((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) .and. use_ipp) then
@@ -442,6 +458,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
, int2_u2_j1b2(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
, 1.d0, tc_grad_square_ao, ao_num*ao_num)
+
+ FREE int2_u2_j1b2
endif
! ---
@@ -478,6 +496,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
call wall_time(time1)
print*, ' Wall time for tc_grad_square_ao = ', time1 - time0
+ call print_memory_usage()
END_PROVIDER
diff --git a/src/non_h_ints_mu/new_grad_tc.irp.f b/src/non_h_ints_mu/new_grad_tc.irp.f
index 24e7e743..499ffe9d 100644
--- a/src/non_h_ints_mu/new_grad_tc.irp.f
+++ b/src/non_h_ints_mu/new_grad_tc.irp.f
@@ -284,6 +284,7 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num,
call wall_time(time1)
print*, ' Wall time for tc_grad_and_lapl_ao = ', time1 - time0
+ call print_memory_usage()
END_PROVIDER
diff --git a/src/non_h_ints_mu/tc_integ.irp.f b/src/non_h_ints_mu/tc_integ.irp.f
index d5995ae5..b2c0df31 100644
--- a/src/non_h_ints_mu/tc_integ.irp.f
+++ b/src/non_h_ints_mu/tc_integ.irp.f
@@ -100,6 +100,8 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
!$OMP END DO
!$OMP END PARALLEL
+ FREE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b
+
elseif(j1b_type .ge. 100) then
PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
@@ -176,6 +178,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
call wall_time(time1)
print*, ' wall time for int2_grad1_u12_ao =', time1-time0
+ call print_memory_usage()
END_PROVIDER
@@ -242,6 +245,8 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
!$OMP END DO
!$OMP END PARALLEL
+ FREE u12sq_j1bsq grad12_j12
+
else
PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
@@ -262,6 +267,8 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
!$OMP END DO
!$OMP END PARALLEL
+ FREE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
+
endif
elseif(j1b_type .ge. 100) then
@@ -324,6 +331,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
call wall_time(time1)
print*, ' wall time for int2_grad1_u12_square_ao =', time1-time0
+ call print_memory_usage()
END_PROVIDER
diff --git a/src/non_h_ints_mu/total_tc_int.irp.f b/src/non_h_ints_mu/total_tc_int.irp.f
index 450bbef0..afa10305 100644
--- a/src/non_h_ints_mu/total_tc_int.irp.f
+++ b/src/non_h_ints_mu/total_tc_int.irp.f
@@ -84,8 +84,11 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao
enddo
endif
+ FREE tc_grad_square_ao tc_grad_and_lapl_ao ao_two_e_coul
+
call wall_time(wall1)
print *, ' wall time for ao_tc_int_chemist ', wall1 - wall0
+ call print_memory_usage()
END_PROVIDER
diff --git a/src/tc_bi_ortho/normal_ordered.irp.f b/src/tc_bi_ortho/normal_ordered.irp.f
index 8adc7a63..cc01d144 100644
--- a/src/tc_bi_ortho/normal_ordered.irp.f
+++ b/src/tc_bi_ortho/normal_ordered.irp.f
@@ -1,3 +1,6 @@
+
+! ---
+
BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
@@ -8,90 +11,116 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_
implicit none
- integer :: i,h1,p1,h2,p2
- integer :: hh1,hh2,pp1,pp2
+ integer :: i, h1, p1, h2, p2
+ integer :: hh1, hh2, pp1, pp2
integer :: Ne(2)
+ double precision :: hthree_aba, hthree_aaa, hthree_aab
+ double precision :: wall0, wall1
integer, allocatable :: occ(:,:)
integer(bit_kind), allocatable :: key_i_core(:,:)
- double precision :: hthree_aba,hthree_aaa,hthree_aab
- double precision :: wall0,wall1
+
+ print*,' Providing normal_two_body_bi_orth ...'
+ call wall_time(wall0)
PROVIDE N_int
- allocate( occ(N_int*bit_kind_size,2) )
- allocate( key_i_core(N_int,2) )
-
- if(core_tc_op) then
- do i = 1, N_int
- key_i_core(i,1) = xor(ref_bitmask(i,1),core_bitmask(i,1))
- key_i_core(i,2) = xor(ref_bitmask(i,2),core_bitmask(i,2))
- enddo
- call bitstring_to_list_ab(key_i_core,occ,Ne,N_int)
- else
- call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int)
- endif
-
- normal_two_body_bi_orth = 0.d0
- print*,'Providing normal_two_body_bi_orth ...'
+ print*,' Providing normal_two_body_bi_orth ...'
call wall_time(wall0)
- !$OMP PARALLEL &
- !$OMP DEFAULT (NONE) &
- !$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, hthree_aba, hthree_aab, hthree_aaa) &
- !$OMP SHARED (N_int, n_act_orb, list_act, Ne, occ, normal_two_body_bi_orth)
- !$OMP DO SCHEDULE (static)
- do hh1 = 1, n_act_orb
- h1 = list_act(hh1)
- do pp1 = 1, n_act_orb
- p1 = list_act(pp1)
- do hh2 = 1, n_act_orb
- h2 = list_act(hh2)
- do pp2 = 1, n_act_orb
- p2 = list_act(pp2)
- ! all contributions from the 3-e terms to the double excitations
- ! s1:(h1-->p1), s2:(h2-->p2) from the HF reference determinant
-
+ if(read_tc_norm_ord) then
- ! opposite spin double excitations : s1 /= s2
- call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aba)
+ open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/normal_two_body_bi_orth', action="read")
+ read(11) normal_two_body_bi_orth
+ close(11)
- ! same spin double excitations : s1 == s2
- if(h1
h2
- ! same spin double excitations with same spin contributions
- if(Ne(2).ge.3)then
- call give_aaa_contraction(N_int, h2, h1, p1, p2, Ne, occ, hthree_aaa) ! exchange h1<->h2
- else
- hthree_aaa = 0.d0
- endif
- else
- ! with opposite spin contributions
- call give_aab_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aab)
- if(Ne(2).ge.3)then
- ! same spin double excitations with same spin contributions
- call give_aaa_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aaa)
- else
- hthree_aaa = 0.d0
- endif
- endif
- normal_two_body_bi_orth(p2,h2,p1,h1) = 0.5d0*(hthree_aba + hthree_aab + hthree_aaa)
+ else
+
+ PROVIDE N_int
+
+ allocate( occ(N_int*bit_kind_size,2) )
+ allocate( key_i_core(N_int,2) )
+
+ if(core_tc_op) then
+ do i = 1, N_int
+ key_i_core(i,1) = xor(ref_bitmask(i,1),core_bitmask(i,1))
+ key_i_core(i,2) = xor(ref_bitmask(i,2),core_bitmask(i,2))
+ enddo
+ call bitstring_to_list_ab(key_i_core,occ,Ne,N_int)
+ else
+ call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int)
+ endif
+
+ normal_two_body_bi_orth = 0.d0
+
+ !$OMP PARALLEL &
+ !$OMP DEFAULT (NONE) &
+ !$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, hthree_aba, hthree_aab, hthree_aaa) &
+ !$OMP SHARED (N_int, n_act_orb, list_act, Ne, occ, normal_two_body_bi_orth)
+ !$OMP DO SCHEDULE (static)
+ do hh1 = 1, n_act_orb
+ h1 = list_act(hh1)
+ do pp1 = 1, n_act_orb
+ p1 = list_act(pp1)
+ do hh2 = 1, n_act_orb
+ h2 = list_act(hh2)
+ do pp2 = 1, n_act_orb
+ p2 = list_act(pp2)
+ ! all contributions from the 3-e terms to the double excitations
+ ! s1:(h1-->p1), s2:(h2-->p2) from the HF reference determinant
+
+
+ ! opposite spin double excitations : s1 /= s2
+ call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aba)
+
+ ! same spin double excitations : s1 == s2
+ if(h1h2
+ ! same spin double excitations with same spin contributions
+ if(Ne(2).ge.3)then
+ call give_aaa_contraction(N_int, h2, h1, p1, p2, Ne, occ, hthree_aaa) ! exchange h1<->h2
+ else
+ hthree_aaa = 0.d0
+ endif
+ else
+ ! with opposite spin contributions
+ call give_aab_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aab)
+ if(Ne(2).ge.3)then
+ ! same spin double excitations with same spin contributions
+ call give_aaa_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aaa)
+ else
+ hthree_aaa = 0.d0
+ endif
+ endif
+ normal_two_body_bi_orth(p2,h2,p1,h1) = 0.5d0*(hthree_aba + hthree_aab + hthree_aaa)
+ enddo
enddo
enddo
enddo
- enddo
- !$OMP END DO
- !$OMP END PARALLEL
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+ deallocate( occ )
+ deallocate( key_i_core )
+ endif
+
+ if(write_tc_norm_ord.and.mpi_master) then
+ open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/normal_two_body_bi_orth', action="write")
+ call ezfio_set_work_empty(.False.)
+ write(11) normal_two_body_bi_orth
+ close(11)
+ call ezfio_set_tc_keywords_io_tc_integ('Read')
+ endif
call wall_time(wall1)
- print*,'Wall time for normal_two_body_bi_orth ',wall1-wall0
+ print*,' Wall time for normal_two_body_bi_orth ', wall1-wall0
- deallocate( occ )
- deallocate( key_i_core )
+ call wall_time(wall1)
+ print*,' Wall time for normal_two_body_bi_orth ', wall1-wall0
END_PROVIDER
-
+! ---
subroutine give_aba_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
@@ -106,30 +135,41 @@ subroutine give_aba_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
!!!! double alpha/beta
hthree = 0.d0
+
do ii = 1, Ne(2) ! purely closed shell part
i = occ(ii,2)
- call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral)
+
+ call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
int_direct = -1.d0 * integral
- call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral)
+
+ call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
int_exc_13 = -1.d0 * integral
- call give_integrals_3_body_bi_ort(p2, i,p1,i,h2,h1,integral)
+
+ call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
int_exc_12 = -1.d0 * integral
- hthree += 2.d0 * int_direct - 1.d0 * ( int_exc_13 + int_exc_12)
+
+ hthree += 2.d0 * int_direct - 1.d0 * (int_exc_13 + int_exc_12)
enddo
+
do ii = Ne(2) + 1, Ne(1) ! purely open-shell part
- i = occ(ii,1)
- call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral)
+ i = occ(ii,1)
+
+ call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
int_direct = -1.d0 * integral
- call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral)
+
+ call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
int_exc_13 = -1.d0 * integral
- call give_integrals_3_body_bi_ort(p2, i,p1,i,h2,h1,integral)
+
+ call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
int_exc_12 = -1.d0 * integral
- hthree += 1.d0 * int_direct - 0.5d0* ( int_exc_13 + int_exc_12)
+
+ hthree += 1.d0 * int_direct - 0.5d0 * (int_exc_13 + int_exc_12)
enddo
-end subroutine give_aba_contraction
-
+ return
+end
+! ---
BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_ab, (mo_num, mo_num, mo_num, mo_num)]
@@ -152,29 +192,31 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_ab, (mo_num, mo_num,
allocate( key_i_core(N_int,2) )
allocate( occ(N_int*bit_kind_size,2) )
- if(core_tc_op)then
- do i = 1, N_int
- key_i_core(i,1) = xor(ref_bitmask(i,1),core_bitmask(i,1))
- key_i_core(i,2) = xor(ref_bitmask(i,2),core_bitmask(i,2))
- enddo
- call bitstring_to_list_ab(key_i_core,occ,Ne,N_int)
+ if(core_tc_op) then
+ do i = 1, N_int
+ key_i_core(i,1) = xor(ref_bitmask(i,1),core_bitmask(i,1))
+ key_i_core(i,2) = xor(ref_bitmask(i,2),core_bitmask(i,2))
+ enddo
+ call bitstring_to_list_ab(key_i_core,occ,Ne,N_int)
else
- call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int)
+ call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int)
endif
+
normal_two_body_bi_orth_ab = 0.d0
do hh1 = 1, n_act_orb
- h1 = list_act(hh1)
- do pp1 = 1, n_act_orb
- p1 = list_act(pp1)
- do hh2 = 1, n_act_orb
- h2 = list_act(hh2)
- do pp2 = 1, n_act_orb
- p2 = list_act(pp2)
- call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree)
- normal_two_body_bi_orth_ab(p2,h2,p1,h1) = hthree
- enddo
+ h1 = list_act(hh1)
+ do pp1 = 1, n_act_orb
+ p1 = list_act(pp1)
+ do hh2 = 1, n_act_orb
+ h2 = list_act(hh2)
+ do pp2 = 1, n_act_orb
+ p2 = list_act(pp2)
+ call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree)
+
+ normal_two_body_bi_orth_ab(p2,h2,p1,h1) = hthree
+ enddo
+ enddo
enddo
- enddo
enddo
deallocate( key_i_core )
@@ -182,7 +224,7 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_ab, (mo_num, mo_num,
END_PROVIDER
-
+! ---
BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_aa_bb, (n_act_orb, n_act_orb, n_act_orb, n_act_orb)]
@@ -250,13 +292,14 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_aa_bb, (n_act_orb, n_
END_PROVIDER
-
+! ---
subroutine give_aaa_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
BEGIN_DOC
-! pure same spin contribution to same spin double excitation s1=h1,p1, s2=h2,p2, with s1==s2
+ ! pure same spin contribution to same spin double excitation s1=h1,p1, s2=h2,p2, with s1==s2
END_DOC
+
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
@@ -270,48 +313,64 @@ subroutine give_aaa_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
hthree = 0.d0
do ii = 1, Ne(2) ! purely closed shell part
i = occ(ii,2)
- call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral)
+
+ call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
int_direct = -1.d0 * integral
- call give_integrals_3_body_bi_ort(p2,p1,i ,i,h2,h1,integral)
+
+ call give_integrals_3_body_bi_ort(p2, p1, i, i, h2, h1, integral)
int_exc_l = -1.d0 * integral
- call give_integrals_3_body_bi_ort(p1,i ,p2,i,h2,h1,integral)
+
+ call give_integrals_3_body_bi_ort(p1, i, p2, i, h2, h1, integral)
int_exc_ll= -1.d0 * integral
- call give_integrals_3_body_bi_ort(p2,i ,p1,i,h2,h1,integral)
+
+ call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
int_exc_12= -1.d0 * integral
- call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral)
+
+ call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
int_exc_13= -1.d0 * integral
- call give_integrals_3_body_bi_ort(i ,p1,p2,i,h2,h1,integral)
+
+ call give_integrals_3_body_bi_ort(i, p1, p2, i, h2, h1, integral)
int_exc_23= -1.d0 * integral
- hthree += 1.d0 * int_direct + int_exc_l + int_exc_ll -( int_exc_12+ int_exc_13+ int_exc_23 )
+ hthree += 1.d0 * int_direct + int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23)
enddo
+
do ii = Ne(2)+1,Ne(1) ! purely open-shell part
i = occ(ii,1)
- call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral)
- int_direct = -1.d0 * integral
- call give_integrals_3_body_bi_ort(p2,p1,i ,i,h2,h1,integral)
- int_exc_l = -1.d0 * integral
- call give_integrals_3_body_bi_ort(p1,i ,p2,i,h2,h1,integral)
- int_exc_ll= -1.d0 * integral
- call give_integrals_3_body_bi_ort(p2,i ,p1,i,h2,h1,integral)
- int_exc_12= -1.d0 * integral
- call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral)
- int_exc_13= -1.d0 * integral
- call give_integrals_3_body_bi_ort(i ,p1,p2,i,h2,h1,integral)
- int_exc_23= -1.d0 * integral
- hthree += 1.d0 * int_direct + 0.5d0 * (int_exc_l + int_exc_ll -( int_exc_12+ int_exc_13+ int_exc_23 ))
+ call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
+ int_direct = -1.d0 * integral
+
+ call give_integrals_3_body_bi_ort(p2, p1, i , i, h2, h1, integral)
+ int_exc_l = -1.d0 * integral
+
+ call give_integrals_3_body_bi_ort(p1, i, p2, i, h2, h1, integral)
+ int_exc_ll = -1.d0 * integral
+
+ call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
+ int_exc_12 = -1.d0 * integral
+
+ call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
+ int_exc_13 = -1.d0 * integral
+
+ call give_integrals_3_body_bi_ort(i, p1, p2, i, h2, h1, integral)
+ int_exc_23 = -1.d0 * integral
+
+ hthree += 1.d0 * int_direct + 0.5d0 * (int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23))
enddo
-end subroutine give_aaa_contraction
-
+ return
+end
+! ---
subroutine give_aab_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
- implicit none
+
use bitmasks ! you need to include the bitmasks_module.f90 features
- integer, intent(in) :: Nint, h1, h2, p1, p2
- integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2)
+
+ implicit none
+ integer, intent(in) :: Nint, h1, h2, p1, p2
+ integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2)
double precision, intent(out) :: hthree
integer :: ii, i
double precision :: int_direct, int_exc_12, int_exc_13, int_exc_23
@@ -320,11 +379,18 @@ subroutine give_aab_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
hthree = 0.d0
do ii = 1, Ne(2) ! purely closed shell part
i = occ(ii,2)
- call give_integrals_3_body_bi_ort(p2,p1,i,h2,h1,i,integral)
+
+ call give_integrals_3_body_bi_ort(p2, p1, i, h2, h1, i, integral)
int_direct = -1.d0 * integral
- call give_integrals_3_body_bi_ort(p1,p2,i,h2,h1,i,integral)
+
+ call give_integrals_3_body_bi_ort(p1, p2, i, h2, h1, i, integral)
int_exc_23= -1.d0 * integral
- hthree += 1.d0 * int_direct - int_exc_23
+
+ hthree += 1.d0 * int_direct - int_exc_23
enddo
-end subroutine give_aab_contraction
+ return
+end
+
+! ---
+
diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg
index 62adb068..a69f5bac 100644
--- a/src/tc_keywords/EZFIO.cfg
+++ b/src/tc_keywords/EZFIO.cfg
@@ -226,6 +226,12 @@ doc: Read/Write integrals int2_grad1_u12_ao, tc_grad_square_ao and tc_grad_and_l
interface: ezfio,provider,ocaml
default: None
+[io_tc_norm_ord]
+type: Disk_access
+doc: Read/Write normal_two_body_bi_orth from/to disk [ Write | Read | None ]
+interface: ezfio,provider,ocaml
+default: None
+
[debug_tc_pt2]
type: integer
doc: If :: 1 then you compute the TC-PT2 the old way, :: 2 then you check with the new version but without three-body
diff --git a/src/tc_scf/rh_tcscf_diis.irp.f b/src/tc_scf/rh_tcscf_diis.irp.f
index 20260a95..0504373c 100644
--- a/src/tc_scf/rh_tcscf_diis.irp.f
+++ b/src/tc_scf/rh_tcscf_diis.irp.f
@@ -11,6 +11,7 @@ subroutine rh_tcscf_diis()
integer :: i, j, it
integer :: dim_DIIS, index_dim_DIIS
+ logical :: converged
double precision :: etc_tot, etc_1e, etc_2e, etc_3e, e_save, e_delta
double precision :: tc_grad, g_save, g_delta, g_delta_th
double precision :: level_shift_save, rate_th
@@ -92,8 +93,9 @@ subroutine rh_tcscf_diis()
PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot
+ converged = .false.
!do while((tc_grad .gt. dsqrt(thresh_tcscf)) .and. (er_DIIS .gt. dsqrt(thresh_tcscf)))
- do while(er_DIIS .gt. dsqrt(thresh_tcscf))
+ do while(.not. converged)
call wall_time(t0)
@@ -218,21 +220,56 @@ subroutine rh_tcscf_diis()
!g_delta_th = dabs(tc_grad) ! g_delta)
er_delta_th = dabs(er_DIIS) !er_delta)
+ converged = er_DIIS .lt. dsqrt(thresh_tcscf)
+
call wall_time(t1)
!write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') &
! it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') &
it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
+
+! Write data in JSON file
+
+ call lock_io
+ if (it == 1) then
+ write(json_unit, json_dict_uopen_fmt)
+ else
+ write(json_unit, json_dict_close_uopen_fmt)
+ endif
+ write(json_unit, json_int_fmt) ' iteration ', it
+ write(json_unit, json_real_fmt) ' SCF TC Energy ', etc_tot
+ write(json_unit, json_real_fmt) ' E(1e) ', etc_1e
+ write(json_unit, json_real_fmt) ' E(2e) ', etc_2e
+ write(json_unit, json_real_fmt) ' E(3e) ', etc_3e
+ write(json_unit, json_real_fmt) ' delta Energy ', e_delta
+ write(json_unit, json_real_fmt) ' DIIS error ', er_DIIS
+ write(json_unit, json_real_fmt) ' level_shift ', level_shift_tcscf
+ write(json_unit, json_real_fmt) ' DIIS ', dim_DIIS
+ write(json_unit, json_real_fmt) ' Wall time (min)', (t1-t0)/60.d0
+ call unlock_io
+
if(er_delta .lt. 0.d0) then
call ezfio_set_tc_scf_bitc_energy(etc_tot)
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
+ write(json_unit, json_true_fmt) 'saved'
+ else
+ write(json_unit, json_false_fmt) 'saved'
endif
+ call lock_io
+ if (converged) then
+ write(json_unit, json_true_fmtx) 'converged'
+ else
+ write(json_unit, json_false_fmtx) 'converged'
+ endif
+ call unlock_io
if(qp_stop()) exit
enddo
+ write(json_unit, json_dict_close_fmtx)
+
! ---
print *, ' TCSCF DIIS converged !'
diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f
index 88ddd26c..04c4f92d 100644
--- a/src/tc_scf/tc_scf.irp.f
+++ b/src/tc_scf/tc_scf.irp.f
@@ -8,6 +8,8 @@ program tc_scf
implicit none
+ write(json_unit,json_array_open_fmt) 'tc-scf'
+
print *, ' starting ...'
my_grid_becke = .True.
@@ -57,6 +59,8 @@ program tc_scf
endif
+ write(json_unit,json_array_close_fmtx)
+ call json_close
end