fixed two rdms

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

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,6 +93,17 @@ 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
@@ -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_bb_st_av_2  += wee_bb(istate)  * state_average_weight(istate)
-   wee_ab_st_av_2  += wee_aa(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

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

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(1,nkeys) = h2
-   keys(2,nkeys) = p2
+   keys(2,nkeys) = h1
-   keys(3,nkeys) = h1
+   keys(3,nkeys) = p2
-   keys(4,nkeys) = h2 
+   keys(4,nkeys) = p1 
  endif
  end
 

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(1,nkeys) = h2
-   keys(2,nkeys) = p2
+   keys(2,nkeys) = h1
-   keys(3,nkeys) = h1
+   keys(3,nkeys) = p2
-   keys(4,nkeys) = h2 
+   keys(4,nkeys) = p1 
  endif
  end
 

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)

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
+!  if (.not.restore_symm) then
-    return
+!    return
-  endif
+!  endif
 
   ! TODO:  Costs O(n^4), but can be improved to (2 n^2 * log(n)):
   ! - copy all values in a 1D array