From d219dc10267c0fe86fbae4683c00f5051229c8c0 Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 11 Jul 2024 13:57:28 +0200 Subject: [PATCH 1/4] beginning to put cholesky in CASSCF --- plugins/local/spher_harm/spher_harm.irp.f | 4 +- .../local/spher_harm/spher_harm_func.irp.f | 13 +++ src/casscf_cipsi/chol_bielec.irp.f | 93 ++++++++++++++++ src/casscf_cipsi/test_chol.irp.f | 23 ++++ src/mo_two_e_ints/cholesky.irp.f | 31 ++++++ src/mu_of_r/f_hf_cholesky.irp.f | 100 ++++++++++++++++++ src/mu_of_r/test_proj_op.irp.f | 20 +++- 7 files changed, 280 insertions(+), 4 deletions(-) create mode 100644 src/casscf_cipsi/chol_bielec.irp.f create mode 100644 src/casscf_cipsi/test_chol.irp.f diff --git a/plugins/local/spher_harm/spher_harm.irp.f b/plugins/local/spher_harm/spher_harm.irp.f index 7a2eea06..e8deafb9 100644 --- a/plugins/local/spher_harm/spher_harm.irp.f +++ b/plugins/local/spher_harm/spher_harm.irp.f @@ -1,7 +1,7 @@ program spher_harm implicit none -! call test_spher_harm + call test_spher_harm ! call test_cart - call test_brutal_spheric +! call test_brutal_spheric end diff --git a/plugins/local/spher_harm/spher_harm_func.irp.f b/plugins/local/spher_harm/spher_harm_func.irp.f index 825bd8ac..f12c8fb9 100644 --- a/plugins/local/spher_harm/spher_harm_func.irp.f +++ b/plugins/local/spher_harm/spher_harm_func.irp.f @@ -7,6 +7,7 @@ subroutine spher_harm_func_r3(r,l,m,re_ylm, im_ylm) double precision :: theta, phi,r_abs call cartesian_to_spherical(r,theta,phi,r_abs) call spher_harm_func(l,m,theta,phi,re_ylm, im_ylm) +! call spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm) end @@ -131,6 +132,10 @@ subroutine spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm) tmp = - inv_sq_pi * dsqrt(3.d0/8.d0) * dsin(theta) re_ylm = tmp * dcos(phi) im_ylm = tmp * dsin(phi) + else if (l==1.and.m==-1)then + tmp = - inv_sq_pi * dsqrt(3.d0/8.d0) * dsin(theta) + re_ylm = tmp * dcos(phi) + im_ylm = -tmp * dsin(phi) else if(l==1.and.m==0)then tmp = inv_sq_pi * dsqrt(3.d0/4.d0) * dcos(theta) re_ylm = tmp @@ -139,10 +144,18 @@ subroutine spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm) tmp = 0.25d0 * inv_sq_pi * dsqrt(0.5d0*15.d0) * dsin(theta)*dsin(theta) re_ylm = tmp * dcos(2.d0*phi) im_ylm = tmp * dsin(2.d0*phi) + else if(l==2.and.m==-2)then + tmp = 0.25d0 * inv_sq_pi * dsqrt(0.5d0*15.d0) * dsin(theta)*dsin(theta) + re_ylm = tmp * dcos(2.d0*phi) + im_ylm =-tmp * dsin(2.d0*phi) else if(l==2.and.m==1)then tmp = - inv_sq_pi * dsqrt(15.d0/8.d0) * dsin(theta) * dcos(theta) re_ylm = tmp * dcos(phi) im_ylm = tmp * dsin(phi) + else if(l==2.and.m==-1)then + tmp = - inv_sq_pi * dsqrt(15.d0/8.d0) * dsin(theta) * dcos(theta) + re_ylm = tmp * dcos(phi) + im_ylm =-tmp * dsin(phi) else if(l==2.and.m==0)then tmp = dsqrt(5.d0/4.d0) * inv_sq_pi* (1.5d0*dcos(theta)*dcos(theta)-0.5d0) re_ylm = tmp diff --git a/src/casscf_cipsi/chol_bielec.irp.f b/src/casscf_cipsi/chol_bielec.irp.f new file mode 100644 index 00000000..1fe985ad --- /dev/null +++ b/src/casscf_cipsi/chol_bielec.irp.f @@ -0,0 +1,93 @@ + +BEGIN_PROVIDER [double precision, cholesky_no_1_idx_transp, (cholesky_mo_num, n_act_orb, mo_num)] + BEGIN_DOC + ! Cholesky vectors with ONE orbital on the active natural orbital basis + END_DOC + implicit none + integer :: i_chol,i_act,i_mo,jj_act + double precision, allocatable :: chol_tmp(:,:) + allocate(chol_tmp(cholesky_mo_num,n_act_orb)) + cholesky_no_1_idx_transp = 0.D0 + do i_mo = 1, mo_num + ! Get all the integrals corresponding to the "i_mo" + do i_act = 1, n_act_orb + jj_act = list_act(i_act) + do i_chol = 1, cholesky_mo_num + chol_tmp(i_chol, i_act) = cholesky_mo_transp(i_chol, jj_act, i_mo) + enddo + enddo +! ! Do the matrix product +! do i_act = 1, n_act_orb +! do jj_act = 1, n_act_orb +! do i_chol = 1, cholesky_mo_num +! cholesky_no_1_idx_transp(i_chol, i_act, i_mo) += chol_tmp(i_chol, jj_act) * natorbsCI(jj_act,i_act) +! enddo +! enddo +! enddo + call dgemm('N','N',cholesky_mo_num,n_act_orb,n_act_orb,1.d0, & + chol_tmp, size(chol_tmp,1), & + natorbsCI, size(natorbsCI,1), & + 0.d0, & + cholesky_no_1_idx_transp(1,1,i_mo), size(cholesky_no_1_idx_transp,1)) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp, (cholesky_mo_num, n_act_orb, n_act_orb)] + BEGIN_DOC + ! Cholesky vectors with TWO orbital on the active natural orbital basis + END_DOC + implicit none + integer :: i_chol,i_act,j_act,jj_act + double precision, allocatable :: chol_tmp(:,:) + allocate(chol_tmp(cholesky_mo_num,n_act_orb)) + cholesky_no_2_idx_transp = 0.D0 + do j_act = 1, n_act_orb + do i_act = 1, n_act_orb + do jj_act = 1, n_act_orb + do i_chol = 1, cholesky_mo_num + cholesky_no_2_idx_transp(i_chol, i_act, j_act) += cholesky_no_1_idx_transp(i_chol, i_act,jj_act) * natorbsCI(jj_act,i_act) + enddo + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp_dgemm, (cholesky_mo_num, n_act_orb, n_act_orb)] + BEGIN_DOC + ! Cholesky vectors with TWO orbital on the active natural orbital basis + END_DOC + implicit none + integer :: i_chol,i_act,j_act,jj_act + double precision, allocatable :: chol_tmp(:,:) + allocate(chol_tmp(cholesky_mo_num,n_act_orb)) + cholesky_no_2_idx_transp_dgemm = 0.D0 + do j_act = 1, n_act_orb + ! Get all the integrals corresponding to the "j_act" + do i_act = 1, n_act_orb + jj_act = list_act(i_act) + do i_chol = 1, cholesky_mo_num + chol_tmp(i_chol, i_act) = cholesky_no_1_idx_transp(i_chol, j_act, jj_act) + enddo + enddo +! ! Do the matrix product +! do i_act = 1, n_act_orb +! do jj_act = 1, n_act_orb +! do i_chol = 1, cholesky_mo_num +! cholesky_no_1_idx_transp(i_chol, i_act, j_act) += chol_tmp(i_chol, jj_act) * natorbsCI(jj_act,i_act) +! enddo +! enddo +! enddo + call dgemm('N','N',cholesky_mo_num,n_act_orb,n_act_orb,1.d0, & + chol_tmp, size(chol_tmp,1), & + natorbsCI, size(natorbsCI,1), & + 0.d0, & + cholesky_no_2_idx_transp_dgemm(1,1,j_act), size(cholesky_no_2_idx_transp_dgemm,1)) + enddo + +END_PROVIDER + + diff --git a/src/casscf_cipsi/test_chol.irp.f b/src/casscf_cipsi/test_chol.irp.f new file mode 100644 index 00000000..b94851f9 --- /dev/null +++ b/src/casscf_cipsi/test_chol.irp.f @@ -0,0 +1,23 @@ +program test_chol + implicit none + read_wf= .True. + touch read_wf + call routine + +end + +subroutine routine + implicit none + integer :: i_chol, i_act, i_mo + double precision :: accu + accu = 0.d0 + do i_mo = 1, n_act_orb + do i_act = 1, n_act_orb + do i_chol = 1, cholesky_mo_num + accu += dabs(cholesky_no_2_idx_transp_dgemm(i_chol,i_act,i_mo) - cholesky_no_2_idx_transp(i_chol,i_act,i_mo)) + print*,cholesky_no_2_idx_transp_dgemm(i_chol,i_act,i_mo) , cholesky_no_2_idx_transp(i_chol,i_act,i_mo) + enddo + enddo + enddo + print*,'accu =', accu +end diff --git a/src/mo_two_e_ints/cholesky.irp.f b/src/mo_two_e_ints/cholesky.irp.f index 7e2c8b37..1fed949d 100644 --- a/src/mo_two_e_ints/cholesky.irp.f +++ b/src/mo_two_e_ints/cholesky.irp.f @@ -101,3 +101,34 @@ BEGIN_PROVIDER [ double precision, cholesky_mo_transp, (cholesky_mo_num, mo_num, END_PROVIDER + +BEGIN_PROVIDER [ double precision, cholesky_semi_mo_transp_simple, (cholesky_mo_num, ao_num, mo_num) ] + implicit none + BEGIN_DOC + ! Cholesky vectors in MO basis + END_DOC + + double precision, allocatable :: X(:,:,:) + double precision :: wall0, wall1 + integer :: ierr + print *, 'Semi AO->MO Transformation of Cholesky vectors' + call wall_time(wall0) + + allocate(X(mo_num,cholesky_mo_num,ao_num), stat=ierr) + if (ierr /= 0) then + print *, irp_here, ': Allocation failed' + endif + integer :: i_chol, i_mo, j_mo, i_ao + cholesky_semi_mo_transp_simple = 0.d0 + do i_mo = 1, mo_num + do i_ao = 1, ao_num + do j_mo = 1, mo_num + do i_chol = 1, cholesky_mo_num + cholesky_semi_mo_transp_simple(i_chol, i_ao,i_mo) += cholesky_mo_transp(i_chol,j_mo,i_mo) * mo_coef_transp(j_mo,i_ao) + enddo + enddo + enddo + enddo + +END_PROVIDER + diff --git a/src/mu_of_r/f_hf_cholesky.irp.f b/src/mu_of_r/f_hf_cholesky.irp.f index 5dd69eb6..179b80dd 100644 --- a/src/mu_of_r/f_hf_cholesky.irp.f +++ b/src/mu_of_r/f_hf_cholesky.irp.f @@ -289,6 +289,106 @@ BEGIN_PROVIDER [ double precision, f_hf_cholesky_sparse, (n_points_final_grid)] endif END_PROVIDER +BEGIN_PROVIDER [ double precision, f_hf_cholesky_sparse_bis, (n_points_final_grid)] + implicit none + integer :: ipoint,m,mm,i,ii,p + !!f(R) = \sum_{I} \sum_{J} Phi_I(R) Phi_J(R) V_IJ + !! = \sum_{I}\sum_{J}\sum_A Phi_I(R) Phi_J(R) V_AI V_AJ + !! = \sum_A \sum_{I}Phi_I(R)V_AI \sum_{J}V_AJ Phi_J(R) + !! = \sum_A V_AR G_AR + !! V_AR = \sum_{I}Phi_IR V_AI = \sum_{I}Phi^t_RI V_AI + double precision :: u_dot_v,wall0,wall1,accu_1, accu_2,mo_i_r1,mo_b_r1 + double precision :: thresh_1,thresh_2 + double precision, allocatable :: accu_vec(:),delta_vec(:) + thresh_2 = ao_cholesky_threshold * 100.d0 + thresh_1 = dsqrt(thresh_2) + provide cholesky_mo_transp + if(elec_alpha_num == elec_beta_num)then + call wall_time(wall0) + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (accu_vec,ipoint,p,ii,i,mm,m,mo_i_r1,mo_b_r1) & + !$OMP ShARED (n_occ_val_orb_for_hf,list_valence_orb_for_hf,mos_in_r_array_omp,aos_in_r_array,thresh_1,thresh_2) & + !$OMP ShARED (cholesky_mo_num,f_hf_cholesky_sparse_bis,n_points_final_grid,cholesky_semi_mo_transp_simple,ao_num) + allocate(accu_vec(cholesky_mo_num)) + !$OMP DO + do ipoint = 1, n_points_final_grid + f_hf_cholesky_sparse_bis(ipoint) = 0.d0 + accu_vec = 0.d0 + do ii = 1, n_occ_val_orb_for_hf(1) + i = list_valence_orb_for_hf(ii,1) + mo_i_r1 = mos_in_r_array_omp(i,ipoint) + if(dabs(mo_i_r1).lt.thresh_1)cycle + do mm = 1, ao_num ! electron 1 + mo_b_r1 = aos_in_r_array(mm,ipoint)*mo_i_r1 + if(dabs(mo_b_r1).lt.thresh_2)cycle + do p = 1, cholesky_mo_num + accu_vec(p) = accu_vec(p) + mo_b_r1 * cholesky_semi_mo_transp_simple(p,mm,i) + enddo + enddo + enddo + do p = 1, cholesky_mo_num + f_hf_cholesky_sparse_bis(ipoint) = f_hf_cholesky_sparse_bis(ipoint) + accu_vec(p) * accu_vec(p) + enddo + f_hf_cholesky_sparse_bis(ipoint) *= 2.D0 + enddo + !$OMP END DO + deallocate(accu_vec) + !$OMP END PARALLEL + + call wall_time(wall1) + print*,'Time to provide f_hf_cholesky_sparse_bis = ',wall1-wall0 + else + call wall_time(wall0) + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (accu_vec,delta_vec,ipoint,p,ii,i,mm,m,mo_i_r1,mo_b_r1) & + !$OMP ShARED (n_occ_val_orb_for_hf,list_valence_orb_for_hf,list_basis,mos_in_r_array_omp,thresh_1,thresh_2) & + !$OMP ShARED (cholesky_mo_num,f_hf_cholesky_sparse_bis,n_points_final_grid,cholesky_mo_transp,n_basis_orb) + allocate(accu_vec(cholesky_mo_num),delta_vec(cholesky_mo_num)) + !$OMP DO + do ipoint = 1, n_points_final_grid + f_hf_cholesky_sparse_bis(ipoint) = 0.d0 + accu_vec = 0.d0 + do ii = 1, n_occ_val_orb_for_hf(2) + i = list_valence_orb_for_hf(ii,2) + mo_i_r1 = mos_in_r_array_omp(i,ipoint) + if(dabs(mo_i_r1).lt.thresh_1)cycle + do mm = 1, n_basis_orb ! electron 1 + m = list_basis(mm) + mo_b_r1 = mos_in_r_array_omp(m,ipoint) + if(dabs(mo_i_r1*mo_b_r1).lt.thresh_2)cycle + do p = 1, cholesky_mo_num + accu_vec(p) = accu_vec(p) + mo_i_r1 * mo_b_r1 * cholesky_mo_transp(p,m,i) + enddo + enddo + enddo + delta_vec = 0.d0 + do ii = n_occ_val_orb_for_hf(2)+1,n_occ_val_orb_for_hf(1) + i = list_valence_orb_for_hf(ii,1) + mo_i_r1 = mos_in_r_array_omp(i,ipoint) + if(dabs(mo_i_r1).lt.thresh_1)cycle + do mm = 1, n_basis_orb ! electron 1 + m = list_basis(mm) + mo_b_r1 = mos_in_r_array_omp(m,ipoint) + if(dabs(mo_i_r1*mo_b_r1).lt.thresh_2)cycle + do p = 1, cholesky_mo_num + delta_vec(p) = delta_vec(p) + mo_i_r1 * mo_b_r1 * cholesky_mo_transp(p,m,i) + enddo + enddo + enddo + do p = 1, cholesky_mo_num + f_hf_cholesky_sparse_bis(ipoint) = f_hf_cholesky_sparse_bis(ipoint) + accu_vec(p) * accu_vec(p) + accu_vec(p) * delta_vec(p) + enddo + f_hf_cholesky_sparse_bis(ipoint) *= 2.D0 + enddo + !$OMP END DO + deallocate(accu_vec) + !$OMP END PARALLEL + call wall_time(wall1) + print*,'Time to provide f_hf_cholesky_sparse_bis = ',wall1-wall0 + endif +END_PROVIDER + + BEGIN_PROVIDER [ double precision, on_top_hf_grid, (n_points_final_grid)] implicit none integer :: ipoint,i,ii diff --git a/src/mu_of_r/test_proj_op.irp.f b/src/mu_of_r/test_proj_op.irp.f index f9aba094..fd5e976b 100644 --- a/src/mu_of_r/test_proj_op.irp.f +++ b/src/mu_of_r/test_proj_op.irp.f @@ -15,7 +15,23 @@ program projected_operators ! call test_f_HF_valence_ab ! call routine_full_mos ! call test_f_ii_valence_ab - call test_f_ia_valence_ab - call test_f_ii_ia_aa_valence_ab +! call test_f_ia_valence_ab +! call test_f_ii_ia_aa_valence_ab + call test end + +subroutine test + implicit none + integer :: i_point + double precision :: ref, new, accu, weight + accu = 0.d0 + do i_point = 1, n_points_final_grid + ref = f_hf_cholesky_sparse(i_point) + new = f_hf_cholesky_sparse_bis(i_point) + weight = final_weight_at_r_vector(i_point) + accu += dabs(ref - new) * weight + enddo + print*,'accu = ',accu + +end From 31ec3ace0540177c2476e478d763a068d73bd41b Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 11 Jul 2024 14:22:27 +0200 Subject: [PATCH 2/4] correct transformation of cholesky vectors on the NO basis --- src/casscf_cipsi/chol_bielec.irp.f | 41 +++++++++++++----------------- src/casscf_cipsi/test_chol.irp.f | 14 +++++++--- 2 files changed, 27 insertions(+), 28 deletions(-) diff --git a/src/casscf_cipsi/chol_bielec.irp.f b/src/casscf_cipsi/chol_bielec.irp.f index 1fe985ad..3104fe5f 100644 --- a/src/casscf_cipsi/chol_bielec.irp.f +++ b/src/casscf_cipsi/chol_bielec.irp.f @@ -34,20 +34,21 @@ BEGIN_PROVIDER [double precision, cholesky_no_1_idx_transp, (cholesky_mo_num, n_ END_PROVIDER -BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp, (cholesky_mo_num, n_act_orb, n_act_orb)] +BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp_old, (cholesky_mo_num, n_act_orb, n_act_orb)] BEGIN_DOC ! Cholesky vectors with TWO orbital on the active natural orbital basis END_DOC implicit none - integer :: i_chol,i_act,j_act,jj_act + integer :: i_chol,i_act,j_act,jj_act,jjj_act double precision, allocatable :: chol_tmp(:,:) allocate(chol_tmp(cholesky_mo_num,n_act_orb)) - cholesky_no_2_idx_transp = 0.D0 - do j_act = 1, n_act_orb - do i_act = 1, n_act_orb - do jj_act = 1, n_act_orb + cholesky_no_2_idx_transp_old = 0.D0 + do jj_act = 1, n_act_orb + jjj_act = list_act(jj_act) + do j_act = 1, n_act_orb + do i_act = 1, n_act_orb do i_chol = 1, cholesky_mo_num - cholesky_no_2_idx_transp(i_chol, i_act, j_act) += cholesky_no_1_idx_transp(i_chol, i_act,jj_act) * natorbsCI(jj_act,i_act) + cholesky_no_2_idx_transp_old(i_chol, i_act, j_act) += cholesky_no_1_idx_transp(i_chol, i_act,jjj_act) * natorbsCI(jj_act,j_act) enddo enddo enddo @@ -56,36 +57,28 @@ BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp, (cholesky_mo_num, n_ END_PROVIDER -BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp_dgemm, (cholesky_mo_num, n_act_orb, n_act_orb)] +BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp, (cholesky_mo_num, n_act_orb, n_act_orb)] BEGIN_DOC ! Cholesky vectors with TWO orbital on the active natural orbital basis END_DOC implicit none integer :: i_chol,i_act,j_act,jj_act - double precision, allocatable :: chol_tmp(:,:) - allocate(chol_tmp(cholesky_mo_num,n_act_orb)) - cholesky_no_2_idx_transp_dgemm = 0.D0 - do j_act = 1, n_act_orb + double precision, allocatable :: chol_tmp(:,:),chol_tmp_bis(:,:) + allocate(chol_tmp(cholesky_mo_num,n_act_orb),chol_tmp_bis(cholesky_mo_num,n_act_orb)) + cholesky_no_2_idx_transp = 0.D0 + do i_act = 1, n_act_orb ! Get all the integrals corresponding to the "j_act" - do i_act = 1, n_act_orb - jj_act = list_act(i_act) + do j_act = 1, n_act_orb + jj_act = list_act(j_act) do i_chol = 1, cholesky_mo_num - chol_tmp(i_chol, i_act) = cholesky_no_1_idx_transp(i_chol, j_act, jj_act) + chol_tmp(i_chol, j_act) = cholesky_no_1_idx_transp(i_chol, i_act, jj_act) enddo enddo -! ! Do the matrix product -! do i_act = 1, n_act_orb -! do jj_act = 1, n_act_orb -! do i_chol = 1, cholesky_mo_num -! cholesky_no_1_idx_transp(i_chol, i_act, j_act) += chol_tmp(i_chol, jj_act) * natorbsCI(jj_act,i_act) -! enddo -! enddo -! enddo call dgemm('N','N',cholesky_mo_num,n_act_orb,n_act_orb,1.d0, & chol_tmp, size(chol_tmp,1), & natorbsCI, size(natorbsCI,1), & 0.d0, & - cholesky_no_2_idx_transp_dgemm(1,1,j_act), size(cholesky_no_2_idx_transp_dgemm,1)) + cholesky_no_2_idx_transp(1,1,i_act), size(cholesky_no_2_idx_transp,1)) enddo END_PROVIDER diff --git a/src/casscf_cipsi/test_chol.irp.f b/src/casscf_cipsi/test_chol.irp.f index b94851f9..8d978817 100644 --- a/src/casscf_cipsi/test_chol.irp.f +++ b/src/casscf_cipsi/test_chol.irp.f @@ -9,15 +9,21 @@ end subroutine routine implicit none integer :: i_chol, i_act, i_mo - double precision :: accu + double precision :: accu,error,exact accu = 0.d0 do i_mo = 1, n_act_orb do i_act = 1, n_act_orb do i_chol = 1, cholesky_mo_num - accu += dabs(cholesky_no_2_idx_transp_dgemm(i_chol,i_act,i_mo) - cholesky_no_2_idx_transp(i_chol,i_act,i_mo)) - print*,cholesky_no_2_idx_transp_dgemm(i_chol,i_act,i_mo) , cholesky_no_2_idx_transp(i_chol,i_act,i_mo) + error = dabs(cholesky_no_2_idx_transp(i_chol,i_act,i_mo) - cholesky_no_2_idx_transp_old(i_chol,i_act,i_mo)) + exact = dabs(cholesky_no_2_idx_transp_old(i_chol,i_act,i_mo)) + accu += error + if(exact.gt.1.d-10)then + if(error/exact.gt.1.d-7)then + write(*,'(4(E16.10,X))')cholesky_no_2_idx_transp(i_chol,i_act,i_mo) , cholesky_no_2_idx_transp_old(i_chol,i_act,i_mo),error,error/exact + endif + endif enddo enddo enddo - print*,'accu =', accu + print*,'accu =', accu/(dble(n_act_orb*n_act_orb*cholesky_mo_num)) end From 56450ed0436c4c1be26f33757ea3a4ca35238b57 Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 11 Jul 2024 19:09:20 +0200 Subject: [PATCH 3/4] introduced functions mimicking the arrays --- src/casscf_cipsi/chol_bielec.irp.f | 158 ++++++++++++++++++++++++++++- src/casscf_cipsi/test_chol.irp.f | 92 ++++++++++++++--- 2 files changed, 234 insertions(+), 16 deletions(-) diff --git a/src/casscf_cipsi/chol_bielec.irp.f b/src/casscf_cipsi/chol_bielec.irp.f index 3104fe5f..94a76453 100644 --- a/src/casscf_cipsi/chol_bielec.irp.f +++ b/src/casscf_cipsi/chol_bielec.irp.f @@ -33,7 +33,6 @@ BEGIN_PROVIDER [double precision, cholesky_no_1_idx_transp, (cholesky_mo_num, n_ END_PROVIDER - BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp_old, (cholesky_mo_num, n_act_orb, n_act_orb)] BEGIN_DOC ! Cholesky vectors with TWO orbital on the active natural orbital basis @@ -83,4 +82,161 @@ BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp, (cholesky_mo_num, n_ END_PROVIDER +BEGIN_PROVIDER [ double precision, cholesky_no_total_transp, (cholesky_mo_num, mo_num, mo_num)] + implicit none + BEGIN_DOC + ! Cholesky vectors defined on all basis including the NO basis + END_DOC + integer :: i_chol, i_act, ii_act, j_act, jj_act, i_core_inact, j_core_inact, ii_core_inact, jj_core_inact + integer :: i_virt, ii_virt, j_virt, jj_virt + ! Block when two orbitals belong to the core/inact + do j_core_inact = 1, n_core_inact_orb + jj_core_inact = list_core_inact(j_core_inact) + do i_core_inact = 1, n_core_inact_orb + ii_core_inact = list_core_inact(i_core_inact) + do i_chol = 1, cholesky_mo_num + cholesky_no_total_transp(i_chol, ii_core_inact, jj_core_inact) = cholesky_mo_transp(i_chol,ii_core_inact,jj_core_inact) + enddo + enddo + enddo + ! Block when one orbitals belongs to the core/inact and one belongs to the active + do j_core_inact = 1, n_core_inact_orb + jj_core_inact = list_core_inact(j_core_inact) + do i_act = 1, n_act_orb + ii_act = list_act(i_act) + do i_chol = 1, cholesky_mo_num + cholesky_no_total_transp(i_chol,ii_act,j_core_inact) = cholesky_no_1_idx_transp(i_chol,i_act,jj_core_inact) + enddo + enddo + enddo + do j_core_inact = 1, n_core_inact_orb + jj_core_inact = list_core_inact(j_core_inact) + do i_act = 1, n_act_orb + ii_act = list_act(i_act) + do i_chol = 1, cholesky_mo_num + cholesky_no_total_transp(i_chol,j_core_inact,ii_act) = cholesky_no_1_idx_transp(i_chol,i_act,jj_core_inact) + enddo + enddo + enddo + + ! Block when two orbitals belong to the active + do j_act = 1, n_act_orb + jj_act = list_act(j_act) + do i_act = 1, n_act_orb + ii_act = list_act(i_act) + do i_chol = 1, cholesky_mo_num + cholesky_no_total_transp(i_chol,ii_act,jj_act) = cholesky_no_2_idx_transp(i_chol,i_act,j_act) + enddo + enddo + enddo + + ! Block when two orbitals belong to the virtuals + do i_virt = 1, n_virt_orb + ii_virt = list_virt(i_virt) + do j_virt = 1, n_virt_orb + jj_virt = list_virt(j_virt) + do i_chol = 1, cholesky_mo_num + cholesky_no_total_transp(i_chol,jj_virt,ii_virt) = cholesky_mo_transp(i_chol,jj_virt,ii_virt) + enddo + enddo + enddo + + ! Block when one orbital is in active and the other in the virtuals + do i_virt = 1, n_virt_orb + ii_virt = list_virt(i_virt) + do i_act = 1, n_act_orb + ii_act = list_act(i_act) + do i_chol = 1, cholesky_mo_num + cholesky_no_total_transp(i_chol,ii_act,ii_virt) = cholesky_no_1_idx_transp(i_chol, i_act,ii_virt) + enddo + enddo + enddo + do i_virt = 1, n_virt_orb + ii_virt = list_virt(i_virt) + do i_act = 1, n_act_orb + ii_act = list_act(i_act) + do i_chol = 1, cholesky_mo_num + cholesky_no_total_transp(i_chol,ii_virt,ii_act) = cholesky_no_1_idx_transp(i_chol, i_act,ii_virt) + enddo + enddo + enddo + ! Block when one orbital is in the virtual and one in the core-inact + do i_virt = 1, n_virt_orb + ii_virt = list_virt(i_virt) + do i_core_inact = 1, n_core_inact_orb + ii_core_inact = list_core_inact(i_core_inact) + do i_chol = 1, cholesky_mo_num + cholesky_no_total_transp(i_chol, ii_core_inact, ii_virt) = cholesky_mo_transp(i_chol, ii_core_inact, ii_virt) + enddo + enddo + enddo + do i_core_inact = 1, n_core_inact_orb + ii_core_inact = list_core_inact(i_core_inact) + do i_virt = 1, n_virt_orb + ii_virt = list_virt(i_virt) + do i_chol = 1, cholesky_mo_num + cholesky_no_total_transp(i_chol, ii_virt, ii_core_inact) = cholesky_mo_transp(i_chol, ii_virt, ii_core_inact) + enddo + enddo + enddo +END_PROVIDER + + +double precision function bielec_no_basis_chol(i_1,j_1,i_2,j_2) + implicit none + integer, intent(in) :: i_1,j_1,i_2,j_2 + BEGIN_DOC + ! integral (i_1 j_1|i_2 j_2) in the mixed basis of both MOs and natural MOs + ! + END_DOC + integer :: i_chol + bielec_no_basis_chol = 0.d0 + do i_chol = 1, cholesky_mo_num + bielec_no_basis_chol += cholesky_no_total_transp(i_chol,i_1, j_1) * cholesky_no_total_transp(i_chol,i_2,j_2) + enddo +end + +double precision function bielec_PQxx_no_chol(i_mo, j_mo, i_ca, j_ca) + implicit none + BEGIN_DOC + ! function that computes (i_mo j_mo| i_ca j_ca) with Cholesky decomposition + ! + ! indices are unshifted orbital numbers + END_DOC + integer, intent(in) :: i_ca, j_ca, i_mo, j_mo + integer :: ii_ca, jj_ca + double precision :: bielec_no_basis_chol + ii_ca = list_core_inact_act(i_ca) + jj_ca = list_core_inact_act(j_ca) + bielec_PQxx_no_chol = bielec_no_basis_chol(i_mo,j_mo,ii_ca,jj_ca) + +end + +double precision function bielec_PxxQ_no_chol(i_mo, j_ca, i_ca, j_mo) + implicit none + BEGIN_DOC + ! function that computes (i_mo j_ca |i_ca j_mo) with Cholesky decomposition + ! + ! indices are unshifted orbital numbers + END_DOC + integer, intent(in) :: i_ca, j_ca, i_mo, j_mo + integer :: ii_ca, jj_ca + double precision :: bielec_no_basis_chol + ii_ca = list_core_inact_act(i_ca) + jj_ca = list_core_inact_act(j_ca) + bielec_PxxQ_no_chol = bielec_no_basis_chol(i_mo, jj_ca, ii_ca, j_mo) + +end + +double precision function bielecCI_no_chol(i_ca, j_ca, k_ca, i_mo) + implicit none + integer, intent(in) :: i_ca, j_ca, k_ca, i_mo + integer :: ii_ca, jj_ca, kk_ca + double precision :: bielec_no_basis_chol + ii_ca = list_act(i_ca) + jj_ca = list_act(j_ca) + kk_ca = list_act(k_ca) + bielecCI_no_chol = bielec_no_basis_chol(ii_ca, jj_ca, kk_ca, i_mo) + +end diff --git a/src/casscf_cipsi/test_chol.irp.f b/src/casscf_cipsi/test_chol.irp.f index 8d978817..87c5c352 100644 --- a/src/casscf_cipsi/test_chol.irp.f +++ b/src/casscf_cipsi/test_chol.irp.f @@ -2,28 +2,90 @@ program test_chol implicit none read_wf= .True. touch read_wf - call routine +! call routine_bielec_PxxQ_no + call routine_bielecCI_no end -subroutine routine +subroutine routine_bielec_PQxx_no implicit none - integer :: i_chol, i_act, i_mo - double precision :: accu,error,exact + integer :: i_chol, i_act, ii_act, j_act, jj_act, i_core_inact, j_core_inact, ii_core_inact, jj_core_inact + integer :: i_virt, ii_virt, j_virt, jj_virt, i_mo, j_mo + double precision :: exact, new, error, accu, bielec_no_basis_chol + double precision :: bielec_PQxx_no_chol + accu = 0.d0 - do i_mo = 1, n_act_orb - do i_act = 1, n_act_orb - do i_chol = 1, cholesky_mo_num - error = dabs(cholesky_no_2_idx_transp(i_chol,i_act,i_mo) - cholesky_no_2_idx_transp_old(i_chol,i_act,i_mo)) - exact = dabs(cholesky_no_2_idx_transp_old(i_chol,i_act,i_mo)) - accu += error - if(exact.gt.1.d-10)then - if(error/exact.gt.1.d-7)then - write(*,'(4(E16.10,X))')cholesky_no_2_idx_transp(i_chol,i_act,i_mo) , cholesky_no_2_idx_transp_old(i_chol,i_act,i_mo),error,error/exact + do i_core_inact = 1, n_core_inact_act_orb + ii_core_inact = list_core_inact_act(i_core_inact) + do j_core_inact = 1, n_core_inact_act_orb + jj_core_inact = list_core_inact_act(j_core_inact) + do i_mo = 1, mo_num + do j_mo = 1, mo_num + exact = bielec_PQxx_no(j_mo,i_mo, j_core_inact, i_core_inact) +! new = bielec_no_basis_chol(j_mo,i_mo, jj_core_inact, ii_core_inact) + new = bielec_PQxx_no_chol(j_mo,i_mo, j_core_inact, i_core_inact) + error = dabs(exact-new) + if(dabs(exact).gt.1.d-10)then + print*,exact,new,error endif - endif + accu += error + enddo enddo enddo enddo - print*,'accu =', accu/(dble(n_act_orb*n_act_orb*cholesky_mo_num)) + print*,'accu = ',accu/(dble(mo_num*mo_num*n_core_inact_act_orb**2)) +end + +subroutine routine_bielec_PxxQ_no + implicit none + integer :: i_chol, i_act, ii_act, j_act, jj_act, i_core_inact, j_core_inact, ii_core_inact, jj_core_inact + integer :: i_virt, ii_virt, j_virt, jj_virt, i_mo, j_mo + double precision :: exact, new, error, accu, bielec_no_basis_chol + double precision :: bielec_PxxQ_no_chol + + accu = 0.d0 + do i_mo = 1, mo_num + do i_core_inact = 1, n_core_inact_act_orb + ii_core_inact = list_core_inact_act(i_core_inact) + do j_core_inact = 1, n_core_inact_act_orb + jj_core_inact = list_core_inact_act(j_core_inact) + do j_mo = 1, mo_num + exact = bielec_PxxQ_no(j_mo, j_core_inact, i_core_inact,i_mo) +! new = bielec_no_basis_chol(j_mo,i_mo, jj_core_inact, ii_core_inact) + new = bielec_PxxQ_no_chol(j_mo, j_core_inact, i_core_inact,i_mo) + error = dabs(exact-new) + accu += error + if(dabs(exact).gt.1.d-10)then + print*,exact,new,error + endif + enddo + enddo + enddo + enddo + print*,'accu = ',accu/(dble(mo_num*mo_num*n_core_inact_act_orb**2)) +end + +subroutine routine_bielecCI_no + implicit none + integer :: i_ca, j_ca, k_ca, i_mo + double precision :: exact, new, error, accu, bielec_no_basis_chol + double precision :: bielecCI_no_chol + + accu = 0.d0 + do i_mo = 1, mo_num + do i_ca = 1, n_act_orb + do j_ca = 1, n_act_orb + do k_ca = 1, n_act_orb + exact =bielecCI_no(k_ca, j_ca, i_ca, i_mo) + new = bielecCI_no_chol(k_ca, j_ca, i_ca, i_mo) + error = dabs(exact-new) + accu += error + if(dabs(exact).gt.1.d-10)then + print*,exact,new,error + endif + enddo + enddo + enddo + enddo + print*,'accu = ',accu/(dble(mo_num*mo_num*n_core_inact_act_orb**2)) end From 505d10084c8331c36ea6a2244bdc47ea8fb33a81 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 12 Jul 2024 16:19:53 +0200 Subject: [PATCH 4/4] Choleskization of the CASSCF --- src/casscf_cipsi/bielec.irp.f | 61 +++++++++----- src/casscf_cipsi/bielec_natorb.irp.f | 68 ++++++++++----- src/casscf_cipsi/casscf.irp.f | 2 +- src/casscf_cipsi/chol_bielec.irp.f | 118 ++++++++++++++------------- src/casscf_cipsi/chol_garb.irp.f | 34 ++++++++ src/casscf_cipsi/gradient.irp.f | 1 + src/casscf_cipsi/hessian.irp.f | 6 ++ src/casscf_cipsi/mcscf_fock.irp.f | 2 + src/casscf_cipsi/test_chol.irp.f | 63 +++++++++----- src/casscf_cipsi/tot_en.irp.f | 1 + 10 files changed, 238 insertions(+), 118 deletions(-) create mode 100644 src/casscf_cipsi/chol_garb.irp.f diff --git a/src/casscf_cipsi/bielec.irp.f b/src/casscf_cipsi/bielec.irp.f index 0a44f994..a4901985 100644 --- a/src/casscf_cipsi/bielec.irp.f +++ b/src/casscf_cipsi/bielec.irp.f @@ -1,18 +1,25 @@ -BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)] +BEGIN_PROVIDER [real*8, bielec_PQxx_array, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)] BEGIN_DOC - ! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active + ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! + ! + ! Replaced by the Cholesky-based function bielec_PQxx + ! + ! bielec_PQxx_array : integral (pq|xx) with p,q arbitrary, x core or active ! indices are unshifted orbital numbers END_DOC implicit none integer :: i,j,ii,jj,p,q,i3,j3,t3,v3 real*8 :: mo_two_e_integral + print*,'' + print*,'Providing bielec_PQxx_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!' + print*,'' - bielec_PQxx(:,:,:,:) = 0.d0 + bielec_PQxx_array(:,:,:,:) = 0.d0 PROVIDE mo_two_e_integrals_in_map !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,ii,j,jj,i3,j3) & - !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PQxx, & + !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PQxx_array, & !$OMP n_act_orb,mo_integrals_map,list_act) !$OMP DO @@ -20,14 +27,14 @@ BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_inact_act_orb,n_core ii=list_core_inact(i) do j=i,n_core_inact_orb jj=list_core_inact(j) - call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j),mo_integrals_map) - bielec_PQxx(:,:,j,i)=bielec_PQxx(:,:,i,j) + call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx_array(1,1,i,j),mo_integrals_map) + bielec_PQxx_array(:,:,j,i)=bielec_PQxx_array(:,:,i,j) end do do j=1,n_act_orb jj=list_act(j) j3=j+n_core_inact_orb - call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j3),mo_integrals_map) - bielec_PQxx(:,:,j3,i)=bielec_PQxx(:,:,i,j3) + call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx_array(1,1,i,j3),mo_integrals_map) + bielec_PQxx_array(:,:,j3,i)=bielec_PQxx_array(:,:,i,j3) end do end do !$OMP END DO @@ -40,8 +47,8 @@ BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_inact_act_orb,n_core do j=i,n_act_orb jj=list_act(j) j3=j+n_core_inact_orb - call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i3,j3),mo_integrals_map) - bielec_PQxx(:,:,j3,i3)=bielec_PQxx(:,:,i3,j3) + call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx_array(1,1,i3,j3),mo_integrals_map) + bielec_PQxx_array(:,:,j3,i3)=bielec_PQxx_array(:,:,i3,j3) end do end do !$OMP END DO @@ -52,9 +59,13 @@ END_PROVIDER -BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)] +BEGIN_PROVIDER [real*8, bielec_PxxQ_array, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)] BEGIN_DOC - ! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active + ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! + ! + ! Replaced by the Cholesky-based function bielec_PxxQ + ! + ! bielec_PxxQ_array : integral (px|xq) with p,q arbitrary, x core or active ! indices are unshifted orbital numbers END_DOC implicit none @@ -62,12 +73,15 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_a double precision, allocatable :: integrals_array(:,:) real*8 :: mo_two_e_integral + print*,'' + print*,'Providing bielec_PxxQ_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!' + print*,'' PROVIDE mo_two_e_integrals_in_map - bielec_PxxQ = 0.d0 + bielec_PxxQ_array = 0.d0 !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,ii,j,jj,i3,j3,integrals_array) & - !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PxxQ, & + !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PxxQ_array, & !$OMP n_act_orb,mo_integrals_map,list_act) allocate(integrals_array(mo_num,mo_num)) @@ -80,8 +94,8 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_a call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map) do q=1,mo_num do p=1,mo_num - bielec_PxxQ(p,i,j,q)=integrals_array(p,q) - bielec_PxxQ(p,j,i,q)=integrals_array(q,p) + bielec_PxxQ_array(p,i,j,q)=integrals_array(p,q) + bielec_PxxQ_array(p,j,i,q)=integrals_array(q,p) end do end do end do @@ -91,8 +105,8 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_a call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map) do q=1,mo_num do p=1,mo_num - bielec_PxxQ(p,i,j3,q)=integrals_array(p,q) - bielec_PxxQ(p,j3,i,q)=integrals_array(q,p) + bielec_PxxQ_array(p,i,j3,q)=integrals_array(p,q) + bielec_PxxQ_array(p,j3,i,q)=integrals_array(q,p) end do end do end do @@ -111,8 +125,8 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_a call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map) do q=1,mo_num do p=1,mo_num - bielec_PxxQ(p,i3,j3,q)=integrals_array(p,q) - bielec_PxxQ(p,j3,i3,q)=integrals_array(q,p) + bielec_PxxQ_array(p,i3,j3,q)=integrals_array(p,q) + bielec_PxxQ_array(p,j3,i3,q)=integrals_array(q,p) end do end do end do @@ -129,10 +143,15 @@ BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)] BEGIN_DOC ! bielecCI : integrals (tu|vp) with p arbitrary, tuv active ! index p runs over the whole basis, t,u,v only over the active orbitals + ! + ! This array can be stored anyway. Ex: 50 active orbitals, 1500 MOs ==> 8x50^3x1500 = 1.5 Gb END_DOC implicit none integer :: i,j,k,p,t,u,v double precision, external :: mo_two_e_integral + double precision :: wall0, wall1 + call wall_time(wall0) + print*,'Providing bielecCI' PROVIDE mo_two_e_integrals_in_map !$OMP PARALLEL DO DEFAULT(NONE) & @@ -151,5 +170,7 @@ BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)] end do end do !$OMP END PARALLEL DO + call wall_time(wall1) + print*,'Time to provide bielecCI = ',wall1 - wall0 END_PROVIDER diff --git a/src/casscf_cipsi/bielec_natorb.irp.f b/src/casscf_cipsi/bielec_natorb.irp.f index 9968530c..99734a0b 100644 --- a/src/casscf_cipsi/bielec_natorb.irp.f +++ b/src/casscf_cipsi/bielec_natorb.irp.f @@ -1,30 +1,38 @@ - BEGIN_PROVIDER [real*8, bielec_PQxx_no, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)] + BEGIN_PROVIDER [real*8, bielec_PQxx_no_array, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)] BEGIN_DOC + ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! + ! + ! Replaced by the Cholesky-based function bielec_PQxx_no + ! ! integral (pq|xx) in the basis of natural MOs ! indices are unshifted orbital numbers + ! END_DOC implicit none integer :: i,j,k,l,t,u,p,q double precision, allocatable :: f(:,:,:), d(:,:,:) + print*,'' + print*,'Providing bielec_PQxx_no_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!' + print*,'' !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(j,k,l,p,d,f) & !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, & - !$OMP bielec_PQxx_no,bielec_PQxx,list_act,natorbsCI) + !$OMP bielec_PQxx_no_array,bielec_PQxx_array,list_act,natorbsCI) allocate (f(n_act_orb,mo_num,n_core_inact_act_orb), & d(n_act_orb,mo_num,n_core_inact_act_orb)) !$OMP DO do l=1,n_core_inact_act_orb - bielec_PQxx_no(:,:,:,l) = bielec_PQxx(:,:,:,l) + bielec_PQxx_no_array(:,:,:,l) = bielec_PQxx_array(:,:,:,l) do k=1,n_core_inact_act_orb do j=1,mo_num do p=1,n_act_orb - f(p,j,k)=bielec_PQxx_no(list_act(p),j,k,l) + f(p,j,k)=bielec_PQxx_no_array(list_act(p),j,k,l) end do end do end do @@ -36,13 +44,13 @@ do k=1,n_core_inact_act_orb do j=1,mo_num do p=1,n_act_orb - bielec_PQxx_no(list_act(p),j,k,l)=d(p,j,k) + bielec_PQxx_no_array(list_act(p),j,k,l)=d(p,j,k) end do end do do j=1,mo_num do p=1,n_act_orb - f(p,j,k)=bielec_PQxx_no(j,list_act(p),k,l) + f(p,j,k)=bielec_PQxx_no_array(j,list_act(p),k,l) end do end do end do @@ -54,7 +62,7 @@ do k=1,n_core_inact_act_orb do p=1,n_act_orb do j=1,mo_num - bielec_PQxx_no(j,list_act(p),k,l)=d(p,j,k) + bielec_PQxx_no_array(j,list_act(p),k,l)=d(p,j,k) end do end do end do @@ -71,7 +79,7 @@ do p=1,n_act_orb do k=1,mo_num do j=1,mo_num - f(j,k,p) = bielec_PQxx_no(j,k,n_core_inact_orb+p,l) + f(j,k,p) = bielec_PQxx_no_array(j,k,n_core_inact_orb+p,l) end do end do end do @@ -83,7 +91,7 @@ do p=1,n_act_orb do k=1,mo_num do j=1,mo_num - bielec_PQxx_no(j,k,n_core_inact_orb+p,l)=d(j,k,p) + bielec_PQxx_no_array(j,k,n_core_inact_orb+p,l)=d(j,k,p) end do end do end do @@ -97,7 +105,7 @@ do p=1,n_act_orb do k=1,mo_num do j=1,mo_num - f(j,k,p) = bielec_PQxx_no(j,k,l,n_core_inact_orb+p) + f(j,k,p) = bielec_PQxx_no_array(j,k,l,n_core_inact_orb+p) end do end do end do @@ -109,7 +117,7 @@ do p=1,n_act_orb do k=1,mo_num do j=1,mo_num - bielec_PQxx_no(j,k,l,n_core_inact_orb+p)=d(j,k,p) + bielec_PQxx_no_array(j,k,l,n_core_inact_orb+p)=d(j,k,p) end do end do end do @@ -123,8 +131,12 @@ END_PROVIDER -BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)] +BEGIN_PROVIDER [real*8, bielec_PxxQ_no_array, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)] BEGIN_DOC + ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! + ! + ! Replaced by the Cholesky-based function bielec_PxxQ_no + ! ! integral (px|xq) in the basis of natural MOs ! indices are unshifted orbital numbers END_DOC @@ -132,10 +144,14 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac integer :: i,j,k,l,t,u,p,q double precision, allocatable :: f(:,:,:), d(:,:,:) + print*,'' + print*,'Providing bielec_PxxQ_no_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!' + print*,'' + !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(j,k,l,p,d,f) & !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, & - !$OMP bielec_PxxQ_no,bielec_PxxQ,list_act,natorbsCI) + !$OMP bielec_PxxQ_no_array,bielec_PxxQ_array,list_act,natorbsCI) allocate (f(n_act_orb,n_core_inact_act_orb,n_core_inact_act_orb), & @@ -143,11 +159,11 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac !$OMP DO do j=1,mo_num - bielec_PxxQ_no(:,:,:,j) = bielec_PxxQ(:,:,:,j) + bielec_PxxQ_no_array(:,:,:,j) = bielec_PxxQ_array(:,:,:,j) do l=1,n_core_inact_act_orb do k=1,n_core_inact_act_orb do p=1,n_act_orb - f(p,k,l) = bielec_PxxQ_no(list_act(p),k,l,j) + f(p,k,l) = bielec_PxxQ_no_array(list_act(p),k,l,j) end do end do end do @@ -159,7 +175,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac do l=1,n_core_inact_act_orb do k=1,n_core_inact_act_orb do p=1,n_act_orb - bielec_PxxQ_no(list_act(p),k,l,j)=d(p,k,l) + bielec_PxxQ_no_array(list_act(p),k,l,j)=d(p,k,l) end do end do end do @@ -176,7 +192,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac do l=1,n_core_inact_act_orb do j=1,mo_num do p=1,n_act_orb - f(p,j,l) = bielec_PxxQ_no(j,n_core_inact_orb+p,l,k) + f(p,j,l) = bielec_PxxQ_no_array(j,n_core_inact_orb+p,l,k) end do end do end do @@ -188,7 +204,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac do l=1,n_core_inact_act_orb do j=1,mo_num do p=1,n_act_orb - bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)=d(p,j,l) + bielec_PxxQ_no_array(j,n_core_inact_orb+p,l,k)=d(p,j,l) end do end do end do @@ -205,7 +221,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac do p=1,n_act_orb do l=1,n_core_inact_act_orb do j=1,mo_num - f(j,l,p) = bielec_PxxQ_no(j,l,n_core_inact_orb+p,k) + f(j,l,p) = bielec_PxxQ_no_array(j,l,n_core_inact_orb+p,k) end do end do end do @@ -217,7 +233,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac do p=1,n_act_orb do l=1,n_core_inact_act_orb do j=1,mo_num - bielec_PxxQ_no(j,l,n_core_inact_orb+p,k)=d(j,l,p) + bielec_PxxQ_no_array(j,l,n_core_inact_orb+p,k)=d(j,l,p) end do end do end do @@ -231,7 +247,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac do p=1,n_act_orb do k=1,n_core_inact_act_orb do j=1,mo_num - f(j,k,p) = bielec_PxxQ_no(j,k,l,n_core_inact_orb+p) + f(j,k,p) = bielec_PxxQ_no_array(j,k,l,n_core_inact_orb+p) end do end do end do @@ -243,7 +259,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac do p=1,n_act_orb do k=1,n_core_inact_act_orb do j=1,mo_num - bielec_PxxQ_no(j,k,l,n_core_inact_orb+p)=d(j,k,p) + bielec_PxxQ_no_array(j,k,l,n_core_inact_orb+p)=d(j,k,p) end do end do end do @@ -259,10 +275,16 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)] BEGIN_DOC ! integrals (tu|vp) in the basis of natural MOs ! index p runs over the whole basis, t,u,v only over the active orbitals + ! + ! This array can be stored anyway. Ex: 50 active orbitals, 1500 MOs ==> 8x50^3x1500 = 1.5 Gb END_DOC implicit none integer :: i,j,k,l,t,u,p,q double precision, allocatable :: f(:,:,:), d(:,:,:) + + double precision :: wall0, wall1 + call wall_time(wall0) + print*,'Providing bielecCI_no' !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(j,k,l,p,d,f) & @@ -363,6 +385,8 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)] deallocate(d,f) !$OMP END PARALLEL + call wall_time(wall1) + print*,'Time to provide bielecCI_no = ',wall1-wall0 END_PROVIDER diff --git a/src/casscf_cipsi/casscf.irp.f b/src/casscf_cipsi/casscf.irp.f index d0a26d36..dc3e2245 100644 --- a/src/casscf_cipsi/casscf.irp.f +++ b/src/casscf_cipsi/casscf.irp.f @@ -11,7 +11,7 @@ program casscf if(small_active_space)then pt2_relative_error = 0.00001 else - thresh_scf = 1.d-4 + thresh_scf = max(1.d-4,thresh_scf) pt2_relative_error = 0.04 endif touch pt2_relative_error diff --git a/src/casscf_cipsi/chol_bielec.irp.f b/src/casscf_cipsi/chol_bielec.irp.f index 94a76453..f69832c1 100644 --- a/src/casscf_cipsi/chol_bielec.irp.f +++ b/src/casscf_cipsi/chol_bielec.irp.f @@ -6,6 +6,9 @@ BEGIN_PROVIDER [double precision, cholesky_no_1_idx_transp, (cholesky_mo_num, n_ implicit none integer :: i_chol,i_act,i_mo,jj_act double precision, allocatable :: chol_tmp(:,:) + double precision :: wall0,wall1 + call wall_time(wall0) + print*,'Providing cholesky_no_1_idx_transp' allocate(chol_tmp(cholesky_mo_num,n_act_orb)) cholesky_no_1_idx_transp = 0.D0 do i_mo = 1, mo_num @@ -16,46 +19,17 @@ BEGIN_PROVIDER [double precision, cholesky_no_1_idx_transp, (cholesky_mo_num, n_ chol_tmp(i_chol, i_act) = cholesky_mo_transp(i_chol, jj_act, i_mo) enddo enddo -! ! Do the matrix product -! do i_act = 1, n_act_orb -! do jj_act = 1, n_act_orb -! do i_chol = 1, cholesky_mo_num -! cholesky_no_1_idx_transp(i_chol, i_act, i_mo) += chol_tmp(i_chol, jj_act) * natorbsCI(jj_act,i_act) -! enddo -! enddo -! enddo call dgemm('N','N',cholesky_mo_num,n_act_orb,n_act_orb,1.d0, & chol_tmp, size(chol_tmp,1), & natorbsCI, size(natorbsCI,1), & 0.d0, & cholesky_no_1_idx_transp(1,1,i_mo), size(cholesky_no_1_idx_transp,1)) enddo + call wall_time(wall1) + print*,'Time to provide cholesky_no_1_idx_transp = ', wall1 - wall0 END_PROVIDER -BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp_old, (cholesky_mo_num, n_act_orb, n_act_orb)] - BEGIN_DOC - ! Cholesky vectors with TWO orbital on the active natural orbital basis - END_DOC - implicit none - integer :: i_chol,i_act,j_act,jj_act,jjj_act - double precision, allocatable :: chol_tmp(:,:) - allocate(chol_tmp(cholesky_mo_num,n_act_orb)) - cholesky_no_2_idx_transp_old = 0.D0 - do jj_act = 1, n_act_orb - jjj_act = list_act(jj_act) - do j_act = 1, n_act_orb - do i_act = 1, n_act_orb - do i_chol = 1, cholesky_mo_num - cholesky_no_2_idx_transp_old(i_chol, i_act, j_act) += cholesky_no_1_idx_transp(i_chol, i_act,jjj_act) * natorbsCI(jj_act,j_act) - enddo - enddo - enddo - enddo - -END_PROVIDER - - BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp, (cholesky_mo_num, n_act_orb, n_act_orb)] BEGIN_DOC ! Cholesky vectors with TWO orbital on the active natural orbital basis @@ -64,6 +38,9 @@ BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp, (cholesky_mo_num, n_ integer :: i_chol,i_act,j_act,jj_act double precision, allocatable :: chol_tmp(:,:),chol_tmp_bis(:,:) allocate(chol_tmp(cholesky_mo_num,n_act_orb),chol_tmp_bis(cholesky_mo_num,n_act_orb)) + double precision :: wall0,wall1 + call wall_time(wall0) + print*,'Providing cholesky_no_2_idx_transp' cholesky_no_2_idx_transp = 0.D0 do i_act = 1, n_act_orb ! Get all the integrals corresponding to the "j_act" @@ -79,6 +56,8 @@ BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp, (cholesky_mo_num, n_ 0.d0, & cholesky_no_2_idx_transp(1,1,i_act), size(cholesky_no_2_idx_transp,1)) enddo + call wall_time(wall1) + print*,'Time to provide cholesky_no_2_idx_transp = ', wall1 - wall0 END_PROVIDER @@ -89,6 +68,9 @@ BEGIN_PROVIDER [ double precision, cholesky_no_total_transp, (cholesky_mo_num, m END_DOC integer :: i_chol, i_act, ii_act, j_act, jj_act, i_core_inact, j_core_inact, ii_core_inact, jj_core_inact integer :: i_virt, ii_virt, j_virt, jj_virt + double precision :: wall0,wall1 + call wall_time(wall0) + print*,'Providing cholesky_no_total_transp ' ! Block when two orbitals belong to the core/inact do j_core_inact = 1, n_core_inact_orb jj_core_inact = list_core_inact(j_core_inact) @@ -180,63 +162,87 @@ BEGIN_PROVIDER [ double precision, cholesky_no_total_transp, (cholesky_mo_num, m enddo enddo enddo + + call wall_time(wall1) + print*,'Time to provide cholesky_no_total_transp = ', wall1 - wall0 END_PROVIDER -double precision function bielec_no_basis_chol(i_1,j_1,i_2,j_2) +double precision function bielec_no_basis(i_1,j_1,i_2,j_2) implicit none integer, intent(in) :: i_1,j_1,i_2,j_2 BEGIN_DOC ! integral (i_1 j_1|i_2 j_2) in the mixed basis of both MOs and natural MOs ! END_DOC - integer :: i_chol - bielec_no_basis_chol = 0.d0 - do i_chol = 1, cholesky_mo_num - bielec_no_basis_chol += cholesky_no_total_transp(i_chol,i_1, j_1) * cholesky_no_total_transp(i_chol,i_2,j_2) + integer :: i + bielec_no_basis = 0.d0 + do i = 1, cholesky_mo_num + bielec_no_basis += cholesky_no_total_transp(i,i_1, j_1) * cholesky_no_total_transp(i,i_2,j_2) enddo end -double precision function bielec_PQxx_no_chol(i_mo, j_mo, i_ca, j_ca) +double precision function bielec_PQxx_no(i_mo, j_mo, i_ca, j_ca) implicit none BEGIN_DOC - ! function that computes (i_mo j_mo| i_ca j_ca) with Cholesky decomposition + ! function that computes (i_mo j_mo| i_ca j_ca) with Cholesky decomposition on the NO basis for active orbitals ! - ! indices are unshifted orbital numbers + ! where i_ca, j_ca are in [1:n_core_inact_act_orb] END_DOC integer, intent(in) :: i_ca, j_ca, i_mo, j_mo integer :: ii_ca, jj_ca - double precision :: bielec_no_basis_chol + double precision :: bielec_no_basis ii_ca = list_core_inact_act(i_ca) jj_ca = list_core_inact_act(j_ca) - bielec_PQxx_no_chol = bielec_no_basis_chol(i_mo,j_mo,ii_ca,jj_ca) - + bielec_PQxx_no = bielec_no_basis(i_mo,j_mo,ii_ca,jj_ca) end -double precision function bielec_PxxQ_no_chol(i_mo, j_ca, i_ca, j_mo) +double precision function bielec_PxxQ_no(i_mo, j_ca, i_ca, j_mo) implicit none BEGIN_DOC - ! function that computes (i_mo j_ca |i_ca j_mo) with Cholesky decomposition + ! function that computes (i_mo j_ca |i_ca j_mo) with Cholesky decomposition on the NO basis for active orbitals ! - ! indices are unshifted orbital numbers + ! where i_ca, j_ca are in [1:n_core_inact_act_orb] END_DOC integer, intent(in) :: i_ca, j_ca, i_mo, j_mo integer :: ii_ca, jj_ca - double precision :: bielec_no_basis_chol + double precision :: bielec_no_basis ii_ca = list_core_inact_act(i_ca) jj_ca = list_core_inact_act(j_ca) - bielec_PxxQ_no_chol = bielec_no_basis_chol(i_mo, jj_ca, ii_ca, j_mo) + bielec_PxxQ_no = bielec_no_basis(i_mo, jj_ca, ii_ca, j_mo) end -double precision function bielecCI_no_chol(i_ca, j_ca, k_ca, i_mo) + +double precision function bielec_PQxx(i_mo, j_mo, i_ca, j_ca) + BEGIN_DOC + ! function that computes (i_mo j_mo |i_ca j_ca) with Cholesky decomposition + ! + ! indices are unshifted orbital numbers + ! + ! where i_ca, j_ca are in [1:n_core_inact_act_orb] + END_DOC implicit none - integer, intent(in) :: i_ca, j_ca, k_ca, i_mo - integer :: ii_ca, jj_ca, kk_ca - double precision :: bielec_no_basis_chol - ii_ca = list_act(i_ca) - jj_ca = list_act(j_ca) - kk_ca = list_act(k_ca) - bielecCI_no_chol = bielec_no_basis_chol(ii_ca, jj_ca, kk_ca, i_mo) - + integer, intent(in) :: i_ca, j_ca, j_mo, i_mo + double precision :: mo_two_e_integral + integer :: ii_ca, jj_ca + ii_ca = list_core_inact_act(i_ca) + jj_ca = list_core_inact_act(j_ca) + bielec_PQxx = mo_two_e_integral(i_mo,ii_ca,j_mo,jj_ca) end + +double precision function bielec_PxxQ(i_mo, i_ca, j_ca, j_mo) + BEGIN_DOC + ! function that computes (i_mo j_mo |i_ca j_ca) with Cholesky decomposition + ! + ! where i_ca, j_ca are in [1:n_core_inact_act_orb] + END_DOC + implicit none + integer, intent(in) :: i_ca, j_ca, j_mo, i_mo + double precision :: mo_two_e_integral + integer :: ii_ca, jj_ca + ii_ca = list_core_inact_act(i_ca) + jj_ca = list_core_inact_act(j_ca) + bielec_PxxQ = mo_two_e_integral(i_mo,jj_ca,ii_ca,j_mo) +end + diff --git a/src/casscf_cipsi/chol_garb.irp.f b/src/casscf_cipsi/chol_garb.irp.f new file mode 100644 index 00000000..c4a8fa59 --- /dev/null +++ b/src/casscf_cipsi/chol_garb.irp.f @@ -0,0 +1,34 @@ + +!!!!! FUNCTIONS THAT WORK BUT WHICH ARE USELESS AS THE ARRAYS CAN ALWAYS BE STORED +!double precision function bielecCI_chol(i_a, j_a, k_a, i_mo) +! BEGIN_DOC +! ! function that computes (i_a j_a |k_a j_mo) with Cholesky decomposition +! ! +! ! where i_a, j_a, k_a are in [1:n_act_orb] !!! ONLY ON ACTIVE +! END_DOC +! implicit none +! integer, intent(in) :: i_a, j_a, k_a, i_mo +! integer :: ii_a, jj_a, kk_a +! double precision :: mo_two_e_integral +! ii_a = list_act(i_a) +! jj_a = list_act(j_a) +! kk_a = list_act(k_a) +! bielecCI_chol = mo_two_e_integral(ii_a,kk_a,jj_a,i_mo) +!end + +!double precision function bielecCI_no_chol(i_ca, j_ca, k_ca, i_mo) +! BEGIN_DOC +! ! function that computes (i_ca j_ca |k_ca j_mo) with Cholesky decomposition on the NO basis for active orbitals +! ! +! ! where i_ca, j_ca, k_ca are in [1:n_core_inact_act_orb] +! END_DOC +! implicit none +! integer, intent(in) :: i_ca, j_ca, k_ca, i_mo +! integer :: ii_ca, jj_ca, kk_ca +! double precision :: bielec_no_basis_chol +! ii_ca = list_act(i_ca) +! jj_ca = list_act(j_ca) +! kk_ca = list_act(k_ca) +! bielecCI_no_chol = bielec_no_basis_chol(ii_ca, jj_ca, kk_ca, i_mo) +! +!end diff --git a/src/casscf_cipsi/gradient.irp.f b/src/casscf_cipsi/gradient.irp.f index a1c5e947..961d260d 100644 --- a/src/casscf_cipsi/gradient.irp.f +++ b/src/casscf_cipsi/gradient.irp.f @@ -157,6 +157,7 @@ real*8 function gradvec_it(i,t) integer :: ii,tt,v,vv,x,y integer :: x3,y3 + double precision :: bielec_PQxx_no ii=list_core_inact(i) tt=list_act(t) diff --git a/src/casscf_cipsi/hessian.irp.f b/src/casscf_cipsi/hessian.irp.f index 458c6aa6..9a7a9031 100644 --- a/src/casscf_cipsi/hessian.irp.f +++ b/src/casscf_cipsi/hessian.irp.f @@ -10,6 +10,7 @@ real*8 function hessmat_itju(i,t,j,u) implicit none integer :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj real*8 :: term,t2 + double precision :: bielec_pqxx_no,bielec_pxxq_no ii=list_core_inact(i) tt=list_act(t) @@ -95,6 +96,7 @@ real*8 function hessmat_itja(i,t,j,a) implicit none integer :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y real*8 :: term + double precision :: bielec_pqxx_no,bielec_pxxq_no ! it/ja ii=list_core_inact(i) @@ -128,6 +130,7 @@ real*8 function hessmat_itua(i,t,u,a) implicit none integer :: i,t,u,a,ii,tt,uu,aa,v,vv,x,xx,u3,t3,v3 real*8 :: term + double precision :: bielec_pqxx_no,bielec_pxxq_no ii=list_core_inact(i) tt=list_act(t) @@ -169,6 +172,7 @@ real*8 function hessmat_iajb(i,a,j,b) implicit none integer :: i,a,j,b,ii,aa,jj,bb real*8 :: term + double precision :: bielec_pqxx_no,bielec_pxxq_no ii=list_core_inact(i) aa=list_virt(a) @@ -205,6 +209,7 @@ real*8 function hessmat_iatb(i,a,t,b) implicit none integer :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3 real*8 :: term + double precision :: bielec_pqxx_no,bielec_pxxq_no ii=list_core_inact(i) aa=list_virt(a) @@ -237,6 +242,7 @@ real*8 function hessmat_taub(t,a,u,b) integer :: t,a,u,b,tt,aa,uu,bb,v,vv,x,xx,y integer :: v3,x3 real*8 :: term,t1,t2,t3 + double precision :: bielec_pqxx_no,bielec_pxxq_no tt=list_act(t) aa=list_virt(a) diff --git a/src/casscf_cipsi/mcscf_fock.irp.f b/src/casscf_cipsi/mcscf_fock.irp.f index 0f4b7a99..82b710a7 100644 --- a/src/casscf_cipsi/mcscf_fock.irp.f +++ b/src/casscf_cipsi/mcscf_fock.irp.f @@ -4,6 +4,7 @@ BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ] END_DOC implicit none integer :: p,q,k,kk,t,tt,u,uu + double precision :: bielec_pxxq_no, bielec_pqxx_no do q=1,mo_num do p=1,mo_num @@ -44,6 +45,7 @@ BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ] END_DOC implicit none integer :: p,q,k,kk,t,tt,u,uu + double precision :: bielec_pxxq_no, bielec_pqxx_no Fapq = 0.d0 diff --git a/src/casscf_cipsi/test_chol.irp.f b/src/casscf_cipsi/test_chol.irp.f index 87c5c352..bcce7cf7 100644 --- a/src/casscf_cipsi/test_chol.irp.f +++ b/src/casscf_cipsi/test_chol.irp.f @@ -3,7 +3,9 @@ program test_chol read_wf= .True. touch read_wf ! call routine_bielec_PxxQ_no - call routine_bielecCI_no +! call routine_bielecCI_no +! call test_bielec_PxxQ_chol +! call test_bielecCI end @@ -12,7 +14,7 @@ subroutine routine_bielec_PQxx_no integer :: i_chol, i_act, ii_act, j_act, jj_act, i_core_inact, j_core_inact, ii_core_inact, jj_core_inact integer :: i_virt, ii_virt, j_virt, jj_virt, i_mo, j_mo double precision :: exact, new, error, accu, bielec_no_basis_chol - double precision :: bielec_PQxx_no_chol + double precision :: bielec_PQxx_no accu = 0.d0 do i_core_inact = 1, n_core_inact_act_orb @@ -21,9 +23,8 @@ subroutine routine_bielec_PQxx_no jj_core_inact = list_core_inact_act(j_core_inact) do i_mo = 1, mo_num do j_mo = 1, mo_num - exact = bielec_PQxx_no(j_mo,i_mo, j_core_inact, i_core_inact) -! new = bielec_no_basis_chol(j_mo,i_mo, jj_core_inact, ii_core_inact) - new = bielec_PQxx_no_chol(j_mo,i_mo, j_core_inact, i_core_inact) + exact = bielec_PQxx_no_array(j_mo,i_mo, j_core_inact, i_core_inact) + new = bielec_PQxx_no(j_mo,i_mo, j_core_inact, i_core_inact) error = dabs(exact-new) if(dabs(exact).gt.1.d-10)then print*,exact,new,error @@ -36,12 +37,12 @@ subroutine routine_bielec_PQxx_no print*,'accu = ',accu/(dble(mo_num*mo_num*n_core_inact_act_orb**2)) end -subroutine routine_bielec_PxxQ_no +subroutine routine_bielec_PxxQ_no_array implicit none integer :: i_chol, i_act, ii_act, j_act, jj_act, i_core_inact, j_core_inact, ii_core_inact, jj_core_inact integer :: i_virt, ii_virt, j_virt, jj_virt, i_mo, j_mo double precision :: exact, new, error, accu, bielec_no_basis_chol - double precision :: bielec_PxxQ_no_chol + double precision :: bielec_PxxQ_no accu = 0.d0 do i_mo = 1, mo_num @@ -50,9 +51,9 @@ subroutine routine_bielec_PxxQ_no do j_core_inact = 1, n_core_inact_act_orb jj_core_inact = list_core_inact_act(j_core_inact) do j_mo = 1, mo_num - exact = bielec_PxxQ_no(j_mo, j_core_inact, i_core_inact,i_mo) + exact = bielec_PxxQ_no_array(j_mo, j_core_inact, i_core_inact,i_mo) ! new = bielec_no_basis_chol(j_mo,i_mo, jj_core_inact, ii_core_inact) - new = bielec_PxxQ_no_chol(j_mo, j_core_inact, i_core_inact,i_mo) + new = bielec_PxxQ_no(j_mo, j_core_inact, i_core_inact,i_mo) error = dabs(exact-new) accu += error if(dabs(exact).gt.1.d-10)then @@ -65,19 +66,43 @@ subroutine routine_bielec_PxxQ_no print*,'accu = ',accu/(dble(mo_num*mo_num*n_core_inact_act_orb**2)) end -subroutine routine_bielecCI_no +subroutine test_bielec_PQxx(i_mo, j_mo, i_ca, j_ca) implicit none - integer :: i_ca, j_ca, k_ca, i_mo - double precision :: exact, new, error, accu, bielec_no_basis_chol - double precision :: bielecCI_no_chol + integer :: i_mo, j_mo, i_ca, j_ca + double precision :: exact, new, error, accu + double precision :: bielec_PQxx accu = 0.d0 - do i_mo = 1, mo_num - do i_ca = 1, n_act_orb - do j_ca = 1, n_act_orb - do k_ca = 1, n_act_orb - exact =bielecCI_no(k_ca, j_ca, i_ca, i_mo) - new = bielecCI_no_chol(k_ca, j_ca, i_ca, i_mo) + do j_ca = 1, n_core_inact_act_orb + do i_ca = 1, n_core_inact_act_orb + do j_mo = 1, mo_num + do i_mo = 1, mo_num + exact = bielec_PQxx_array(i_mo, j_mo, i_ca, j_ca) + new = bielec_PQxx(i_mo, j_mo, i_ca, j_ca) + error = dabs(exact-new) + accu += error + if(dabs(exact).gt.1.d-10)then + print*,exact,new,error + endif + enddo + enddo + enddo + enddo + print*,'accu = ',accu/(dble(mo_num*mo_num*n_core_inact_act_orb**2)) +end + +subroutine test_bielec_PxxQ_chol(i_mo, i_ca, j_ca, j_mo) + implicit none + integer :: i_mo, i_ca, j_ca, j_mo + double precision :: exact, new, error, accu + double precision :: bielec_PxxQ + accu = 0.d0 + do j_mo = 1, mo_num + do j_ca = 1, n_core_inact_act_orb + do i_ca =1, n_core_inact_act_orb + do i_mo = 1, mo_num + exact = bielec_PxxQ_array(i_mo, i_ca, j_ca, j_mo) + new = bielec_PxxQ(i_mo, i_ca, j_ca, j_mo) error = dabs(exact-new) accu += error if(dabs(exact).gt.1.d-10)then diff --git a/src/casscf_cipsi/tot_en.irp.f b/src/casscf_cipsi/tot_en.irp.f index 1d70e087..37ceac05 100644 --- a/src/casscf_cipsi/tot_en.irp.f +++ b/src/casscf_cipsi/tot_en.irp.f @@ -8,6 +8,7 @@ implicit none integer :: t,u,v,x,i,ii,tt,uu,vv,xx,j,jj,t3,u3,v3,x3 real*8 :: e_one_all,e_two_all + double precision :: bielec_PQxx,bielec_PxxQ e_one_all=0.D0 e_two_all=0.D0 do i=1,n_core_inact_orb