From ca2b58b495eaa783b9c6857857bced0f811f235d Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 3 Jun 2021 14:57:45 +0200 Subject: [PATCH] fixed two rdms --- src/dft_utils_in_r/mo_in_r.irp.f | 18 +++++ src/two_body_rdm/example.irp.f | 72 +++++++++++++++++-- src/two_body_rdm/state_av_act_2rdm.irp.f | 6 +- src/two_rdm_routines/update_rdm.irp.f | 30 ++++++-- .../update_state_av_rdm.irp.f | 25 +++++-- src/utils/constants.include.F | 1 + src/utils/linear_algebra.irp.f | 6 +- 7 files changed, 138 insertions(+), 20 deletions(-) diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f index d4afcf76..b901ea50 100644 --- a/src/dft_utils_in_r/mo_in_r.irp.f +++ b/src/dft_utils_in_r/mo_in_r.irp.f @@ -89,6 +89,24 @@ enddo END_PROVIDER + BEGIN_PROVIDER[double precision, mos_grad_in_r_array_transp_bis, (n_points_final_grid,mo_num,3)] + implicit none + BEGIN_DOC +! Transposed gradients +! + END_DOC + integer :: i,j,m + do m = 1, 3 + do j = 1, mo_num + do i = 1, n_points_final_grid + mos_grad_in_r_array_transp_bis(i,j,m) = mos_grad_in_r_array(j,i,m) + enddo + enddo + enddo + END_PROVIDER + + + BEGIN_PROVIDER [double precision, alpha_dens_kin_in_r, (n_points_final_grid)] &BEGIN_PROVIDER [double precision, beta_dens_kin_in_r, (n_points_final_grid)] implicit none diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index e9cbd1a2..6cf0209e 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -15,7 +15,7 @@ subroutine routine_active_only double precision :: wee_aa_st_av, rdm_aa_st_av double precision :: wee_bb_st_av, rdm_bb_st_av double precision :: wee_ab_st_av, rdm_ab_st_av - double precision :: wee_tot_st_av, rdm_tot_st_av + double precision :: wee_tot_st_av, rdm_tot_st_av,spin_trace double precision :: wee_aa_st_av_2,wee_ab_st_av_2,wee_bb_st_av_2,wee_tot_st_av_2,wee_tot_st_av_3 wee_ab = 0.d0 @@ -38,6 +38,27 @@ subroutine routine_active_only provide act_2_rdm_ab_mo act_2_rdm_aa_mo act_2_rdm_bb_mo act_2_rdm_spin_trace_mo provide state_av_act_2_rdm_ab_mo state_av_act_2_rdm_aa_mo provide state_av_act_2_rdm_bb_mo state_av_act_2_rdm_spin_trace_mo + i = 1 + j = 2 +! print*,'testing stuffs' +! istate = 1 +! print*,'alpha/beta' +! print*,'',j,i,j,i +! print*,act_2_rdm_ab_mo(j,i,j,i,istate) +! print*,'',i,j,i,j +! print*,act_2_rdm_ab_mo(i,j,i,j,istate) +! print*,'alpha/alpha' +! print*,'',j,i,j,i +! print*,act_2_rdm_aa_mo(j,i,j,i,istate) +! print*,'',i,j,i,j +! print*,act_2_rdm_aa_mo(i,j,i,j,istate) +! print*,'spin_trace' +! print*,'',j,i,j,i +! print*,act_2_rdm_spin_trace_mo(j,i,j,i,istate) +! print*,'',i,j,i,j +! print*,act_2_rdm_spin_trace_mo(i,j,i,j,istate) +! stop +! print*,'**************************' print*,'**************************' do istate = 1, N_states @@ -51,6 +72,19 @@ subroutine routine_active_only korb = list_act(k) do l = 1, n_act_orb lorb = list_act(l) + if(dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(j,i,l,k,istate)).gt.1.d-10)then + print*,'Error in act_2_rdm_spin_trace_mo' + print*,"dabs(act_2_rdm_spin_trace_mo(i,j,k,l) - act_2_rdm_spin_trace_mo(j,i,l,k)).gt.1.d-10" + print*,i,j,k,l + print*,act_2_rdm_spin_trace_mo(i,j,k,l,istate),act_2_rdm_spin_trace_mo(j,i,l,k,istate),dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(j,i,l,k,istate)) + endif + + if(dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(k,l,i,j,istate)).gt.1.d-10)then + print*,'Error in act_2_rdm_spin_trace_mo' + print*,"dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(k,l,i,j,istate),istate).gt.1.d-10" + print*,i,j,k,l + print*,act_2_rdm_spin_trace_mo(i,j,k,l,istate),act_2_rdm_spin_trace_mo(k,l,i,j,istate),dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(k,l,i,j,istate)) + endif vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) @@ -59,7 +93,18 @@ subroutine routine_active_only rdmaa = act_2_rdm_aa_mo(l,k,j,i,istate) rdmbb = act_2_rdm_bb_mo(l,k,j,i,istate) rdmtot = act_2_rdm_spin_trace_mo(l,k,j,i,istate) + spin_trace = rdmaa + rdmbb + rdmab + if(dabs(rdmtot- spin_trace).gt.1.d-10)then + print*,'Error in non state average !!!!' + print*,l,k,j,i + print*,lorb,korb,jorb,iorb + print*,spin_trace,rdmtot,dabs(spin_trace - rdmtot) + print*,'rdmab,rdmaa,rdmbb' + print*, rdmab,rdmaa,rdmbb + + endif + wee_ab(istate) += vijkl * rdmab wee_aa(istate) += vijkl * rdmaa @@ -71,8 +116,8 @@ subroutine routine_active_only enddo enddo wee_aa_st_av_2 += wee_aa(istate) * state_average_weight(istate) - wee_bb_st_av_2 += wee_aa(istate) * state_average_weight(istate) - wee_ab_st_av_2 += wee_aa(istate) * state_average_weight(istate) + wee_bb_st_av_2 += wee_bb(istate) * state_average_weight(istate) + wee_ab_st_av_2 += wee_ab(istate) * state_average_weight(istate) wee_tot_st_av_2 += wee_tot(istate) * state_average_weight(istate) wee_tot_st_av_3 += psi_energy_two_e(istate) * state_average_weight(istate) print*,'' @@ -87,7 +132,6 @@ subroutine routine_active_only print*,'Full energy ' print*,'psi_energy_two_e(istate)= ',psi_energy_two_e(istate) enddo - wee_aa_st_av = 0.d0 wee_bb_st_av = 0.d0 wee_ab_st_av = 0.d0 @@ -103,10 +147,30 @@ subroutine routine_active_only vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) + if(dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(j,i,l,k)).gt.1.d-10)then + print*,'Error in state_av_act_2_rdm_spin_trace_mo' + print*,"dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(j,i,l,k)).gt.1.d-10" + print*,i,j,k,l + print*,state_av_act_2_rdm_spin_trace_mo(i,j,k,l),state_av_act_2_rdm_spin_trace_mo(j,i,l,k),dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(j,i,l,k)) + endif + + if(dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(k,l,i,j)).gt.1.d-10)then + print*,'Error in state_av_act_2_rdm_spin_trace_mo' + print*,"dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(k,l,i,j)).gt.1.d-10" + print*,i,j,k,l + print*,state_av_act_2_rdm_spin_trace_mo(i,j,k,l),state_av_act_2_rdm_spin_trace_mo(k,l,i,j),dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(k,l,i,j)) + endif + rdm_aa_st_av = state_av_act_2_rdm_aa_mo(l,k,j,i) rdm_bb_st_av = state_av_act_2_rdm_bb_mo(l,k,j,i) rdm_ab_st_av = state_av_act_2_rdm_ab_mo(l,k,j,i) + spin_trace = rdm_aa_st_av + rdm_bb_st_av + rdm_ab_st_av rdm_tot_st_av = state_av_act_2_rdm_spin_trace_mo(l,k,j,i) + if(dabs(spin_trace - rdm_tot_st_av).gt.1.d-10)then + print*,'Error !!!!' + print*,l,k,j,i + print*,spin_trace,rdm_tot_st_av,dabs(spin_trace - rdm_tot_st_av) + endif wee_aa_st_av += vijkl * rdm_aa_st_av wee_bb_st_av += vijkl * rdm_bb_st_av diff --git a/src/two_body_rdm/state_av_act_2rdm.irp.f b/src/two_body_rdm/state_av_act_2rdm.irp.f index d85c3cdb..db793047 100644 --- a/src/two_body_rdm/state_av_act_2rdm.irp.f +++ b/src/two_body_rdm/state_av_act_2rdm.irp.f @@ -119,7 +119,11 @@ call wall_time(wall_1) double precision :: wall_1, wall_2 print*,'providing state_av_act_2_rdm_spin_trace_mo ' - call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + state_av_act_2_rdm_spin_trace_mo = state_av_act_2_rdm_ab_mo & + + state_av_act_2_rdm_aa_mo & + + state_av_act_2_rdm_bb_mo + +! call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call wall_time(wall_2) print*,'Time to provide state_av_act_2_rdm_spin_trace_mo',wall_2 - wall_1 diff --git a/src/two_rdm_routines/update_rdm.irp.f b/src/two_rdm_routines/update_rdm.irp.f index 662b3ee9..8aeb0699 100644 --- a/src/two_rdm_routines/update_rdm.irp.f +++ b/src/two_rdm_routines/update_rdm.irp.f @@ -61,14 +61,24 @@ ! Therefore you don't necessayr have symmetry between electron 1 and 2 nkeys += 1 do istate = 1, N_st - values(istate,nkeys) = c_1(istate) + values(istate,nkeys) = 0.5d0 * c_1(istate) enddo keys(1,nkeys) = h1 keys(2,nkeys) = h2 keys(3,nkeys) = h1 keys(4,nkeys) = h2 + + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 enddo enddo + else if (alpha_alpha)then do i = 1, n_occ_ab(1) i1 = occ(i,1) @@ -259,12 +269,20 @@ if(alpha_beta)then nkeys += 1 do istate = 1, N_st - values(istate,nkeys) = c_1(istate) * phase + values(istate,nkeys) = 0.5d0 * c_1(istate) * phase enddo keys(1,nkeys) = h1 keys(2,nkeys) = h2 keys(3,nkeys) = p1 keys(4,nkeys) = p2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 else if(spin_trace)then nkeys += 1 do istate = 1, N_st @@ -278,10 +296,10 @@ do istate = 1, N_st values(istate,nkeys) = 0.5d0 * c_1(istate) * phase enddo - keys(1,nkeys) = p1 - keys(2,nkeys) = p2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 endif end diff --git a/src/two_rdm_routines/update_state_av_rdm.irp.f b/src/two_rdm_routines/update_state_av_rdm.irp.f index bf7b6148..0038e94b 100644 --- a/src/two_rdm_routines/update_state_av_rdm.irp.f +++ b/src/two_rdm_routines/update_state_av_rdm.irp.f @@ -60,11 +60,18 @@ ! If alpha/beta, electron 1 is alpha, electron 2 is beta ! Therefore you don't necessayr have symmetry between electron 1 and 2 nkeys += 1 - values(nkeys) = c_1 + values(nkeys) = 0.5d0 * c_1 keys(1,nkeys) = h1 keys(2,nkeys) = h2 keys(3,nkeys) = h1 keys(4,nkeys) = h2 + + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 enddo enddo else if (alpha_alpha)then @@ -236,11 +243,17 @@ p2 = list_orb_reverse(p2) if(alpha_beta)then nkeys += 1 - values(nkeys) = c_1 * phase + values(nkeys) = 0.5d0 * c_1 * phase keys(1,nkeys) = h1 keys(2,nkeys) = h2 keys(3,nkeys) = p1 keys(4,nkeys) = p2 + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 else if(spin_trace)then nkeys += 1 values(nkeys) = 0.5d0 * c_1 * phase @@ -250,10 +263,10 @@ keys(4,nkeys) = p2 nkeys += 1 values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = p1 - keys(2,nkeys) = p2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 endif end diff --git a/src/utils/constants.include.F b/src/utils/constants.include.F index 7399b4a6..297a839e 100644 --- a/src/utils/constants.include.F +++ b/src/utils/constants.include.F @@ -3,6 +3,7 @@ integer, parameter :: SIMD_vector = 32 integer, parameter :: N_int_max = 32 double precision, parameter :: pi = dacos(-1.d0) +double precision, parameter :: inv_pi = 1.d0/dacos(-1.d0) double precision, parameter :: sqpi = dsqrt(dacos(-1.d0)) double precision, parameter :: pi_5_2 = 34.9868366552d0 double precision, parameter :: dfour_pi = 4.d0*dacos(-1.d0) diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 5b2f9067..405d2d20 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1589,9 +1589,9 @@ subroutine restore_symmetry(m,n,A,LDA,thresh) thresh2 = dsqrt(thresh) call nullify_small_elements(m,n,A,LDA,thresh) - if (.not.restore_symm) then - return - endif +! if (.not.restore_symm) then +! return +! endif ! TODO: Costs O(n^4), but can be improved to (2 n^2 * log(n)): ! - copy all values in a 1D array