From 36fda1d1bbff650e5bb828169c5fb4ce6795dfef Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 10 Nov 2023 14:23:26 +0100 Subject: [PATCH] Fix density matrices in SR correction --- devel/sr_correction/alpha.irp.f | 2 +- devel/sr_correction/density_matrices.irp.f | 4 +- devel/sr_correction/energy_mu.irp.f | 8 ++-- stable/qmcchem/pot_ao_pseudo_ints.irp.f | 54 +++------------------- 4 files changed, 14 insertions(+), 54 deletions(-) diff --git a/devel/sr_correction/alpha.irp.f b/devel/sr_correction/alpha.irp.f index 3ce05c3..8d198f5 100644 --- a/devel/sr_correction/alpha.irp.f +++ b/devel/sr_correction/alpha.irp.f @@ -82,7 +82,7 @@ BEGIN_PROVIDER [ double precision, correction_alpha_0_r, (N_states) ] correction_alpha_0_r(istate) = 0.d0 do l=1,mo_num do k=1,mo_num - correction_alpha_0_r += c0 * ddot (mo_num*mo_num, full_occ_2_rdm_ab_mo(1,1,k,l,istate), 1, W_bar_mu(1,1,k,l), 1) + correction_alpha_0_r += c0 * ddot (mo_num*mo_num, full_occ_2_rdm_spin_trace_mo(1,1,k,l,istate), 1, W_bar_mu(1,1,k,l), 1) enddo enddo enddo diff --git a/devel/sr_correction/density_matrices.irp.f b/devel/sr_correction/density_matrices.irp.f index 21672d7..6a0695b 100644 --- a/devel/sr_correction/density_matrices.irp.f +++ b/devel/sr_correction/density_matrices.irp.f @@ -10,7 +10,7 @@ BEGIN_PROVIDER [ double precision, two_e_dm_mo_singlet, (mo_num, mo_num, mo_num, do k=1,mo_num do j=1,mo_num do i=1,mo_num - two_e_dm_mo_singlet(i,j,k,l,istate) = 0.5d0 * (full_occ_2_rdm_ab_mo(i,j,k,l,istate) + full_occ_2_rdm_ab_mo(i,j,l,k,istate)) + two_e_dm_mo_singlet(i,j,k,l,istate) = 0.5d0 * (full_occ_2_rdm_spin_trace_mo(i,j,k,l,istate) + full_occ_2_rdm_spin_trace_mo(i,j,l,k,istate)) enddo enddo enddo @@ -30,7 +30,7 @@ BEGIN_PROVIDER [ double precision, two_e_dm_mo_triplet, (mo_num, mo_num, mo_num, do k=1,mo_num do j=1,mo_num do i=1,mo_num - two_e_dm_mo_triplet(i,j,k,l,istate) = 0.5d0 * (full_occ_2_rdm_ab_mo(i,j,k,l,istate) - full_occ_2_rdm_ab_mo(i,j,l,k,istate)) + two_e_dm_mo_triplet(i,j,k,l,istate) = 0.5d0 * (full_occ_2_rdm_spin_trace_mo(i,j,k,l,istate) - full_occ_2_rdm_spin_trace_mo(i,j,l,k,istate)) enddo enddo enddo diff --git a/devel/sr_correction/energy_mu.irp.f b/devel/sr_correction/energy_mu.irp.f index 878db6d..2f105d4 100644 --- a/devel/sr_correction/energy_mu.irp.f +++ b/devel/sr_correction/energy_mu.irp.f @@ -73,7 +73,7 @@ BEGIN_PROVIDER [ double precision, energy_mu, (N_states) ] one_e += ddot(mo_num, one_e_dm_mo_alpha(1,l,istate), 1, mo_one_e_integrals(1,l), 1) one_e += ddot(mo_num, one_e_dm_mo_beta (1,l,istate), 1, mo_one_e_integrals(1,l), 1) do k=1,mo_num - two_e += 0.5d0 * ddot (mo_num*mo_num, full_occ_2_rdm_ab_mo(1,1,k,l,istate), 1, W_mu(1,1,k,l), 1) + two_e += 0.5d0 * ddot (mo_num*mo_num, full_occ_2_rdm_spin_trace_mo(1,1,k,l,istate), 1, W_mu(1,1,k,l), 1) enddo enddo energy_mu(istate) = one_e + two_e + nuclear_repulsion @@ -95,7 +95,7 @@ BEGIN_PROVIDER [ double precision, correction_alpha_0, (N_states) ] correction_alpha_0(istate) = 0.d0 do l=1,mo_num do k=1,mo_num - correction_alpha_0(istate) += c0 * ddot (mo_num*mo_num, full_occ_2_rdm_ab_mo(1,1,k,l,istate), 1, W_bar_mu(1,1,k,l), 1) + correction_alpha_0(istate) += c0 * ddot (mo_num*mo_num, full_occ_2_rdm_spin_trace_mo(1,1,k,l,istate), 1, W_bar_mu(1,1,k,l), 1) enddo enddo enddo @@ -142,7 +142,7 @@ BEGIN_PROVIDER [ double precision, correction_mu, (N_states) ] correction_mu(istate) = 0.d0 do l=1,mo_num do k=1,mo_num - correction_mu(istate) += c0 * ddot (mo_num*mo_num, full_occ_2_rdm_ab_mo(1,1,k,l,istate), 1, W_bar_mu(1,1,k,l), 1) + correction_mu(istate) += c0 * ddot (mo_num*mo_num, full_occ_2_rdm_spin_trace_mo(1,1,k,l,istate), 1, W_bar_mu(1,1,k,l), 1) enddo enddo enddo @@ -163,7 +163,7 @@ BEGIN_PROVIDER [ double precision, correction_mu0, (N_states) ] correction_mu0(istate) = 0.d0 do l=1,mo_num do k=1,mo_num - correction_mu0(istate) += c0 * ddot (mo_num*mo_num, full_occ_2_rdm_ab_mo(1,1,k,l,istate), 1, W_bar_mu0(1,1,k,l), 1) + correction_mu0(istate) += c0 * ddot (mo_num*mo_num, full_occ_2_rdm_spin_trace_mo(1,1,k,l,istate), 1, W_bar_mu0(1,1,k,l), 1) enddo enddo enddo diff --git a/stable/qmcchem/pot_ao_pseudo_ints.irp.f b/stable/qmcchem/pot_ao_pseudo_ints.irp.f index eacb968..8c56b97 100644 --- a/stable/qmcchem/pot_ao_pseudo_ints.irp.f +++ b/stable/qmcchem/pot_ao_pseudo_ints.irp.f @@ -29,6 +29,8 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_l enddo ao_pseudo_grid = 0.d0 + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP PRIVATE(j,k,c,l,i,a,n_a,g_a,m,y) do j=1,pseudo_grid_size do k=1,nucl_num c(1:3) = nucl_coord(k,1:3) @@ -48,6 +50,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_l enddo enddo enddo + !$OMP END PARALLEL DO deallocate(r) END_PROVIDER @@ -75,54 +78,11 @@ BEGIN_PROVIDER [ double precision, mo_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_l double precision :: y double precision, allocatable :: r(:) - allocate (r(pseudo_grid_size)) - dr = pseudo_grid_rmax/dble(pseudo_grid_size) - r(1) = 0.d0 - do j=2,pseudo_grid_size - r(j) = r(j-1) + dr - enddo - - mo_pseudo_grid = 0.d0 - do n=1,pseudo_grid_size - do k=1,nucl_num - do l=0,pseudo_lmax - do m=-l,l - do i=1,ao_num - if (dabs(ao_pseudo_grid(i,m,l,k,n)) < 1.e-12) then - cycle - endif - do j=1,mo_num - mo_pseudo_grid(j,m,l,k,n) = mo_pseudo_grid(j,m,l,k,n) + & - ao_pseudo_grid(i,m,l,k,n) * mo_coef(i,j) - enddo - enddo - enddo - enddo - enddo - enddo - deallocate(r) + mo_pseudo_grid = 0.d0 + call dgemm('T','N', mo_num, (2*pseudo_lmax+1)*(pseudo_lmax+1)*nucl_num*pseudo_grid_size, & + ao_num, 1.d0, mo_coef, size(mo_coef,1), ao_pseudo_grid, & + size(ao_pseudo_grid,1), 0.d0, mo_pseudo_grid, size(mo_pseudo_grid,1)) END_PROVIDER -double precision function test_pseudo_grid_ao(i,j) - implicit none - integer, intent(in) :: i,j - integer :: k,l,m,n - double precision :: r, dr,u - dr = pseudo_grid_rmax/dble(pseudo_grid_size) - - test_pseudo_grid_ao = 0.d0 - r = 0.d0 - do k=1,pseudo_grid_size - do n=1,nucl_num - do l = 0,pseudo_lmax - u = pseudo_v_kl(n,l,1) * exp(-pseudo_dz_kl(n,l,1)*r*r)* r*r*dr - do m=-l,l - test_pseudo_grid_ao += ao_pseudo_grid(i,m,l,n,k) * ao_pseudo_grid(j,m,l,n,k) * u - enddo - enddo - enddo - r = r+dr - enddo -end