From 31ec3ace0540177c2476e478d763a068d73bd41b Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 11 Jul 2024 14:22:27 +0200 Subject: [PATCH] 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