Fix density matrices in SR correction

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

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

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

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
+  call dgemm('T','N', mo_num, (2*pseudo_lmax+1)*(pseudo_lmax+1)*nucl_num*pseudo_grid_size, &
-   do k=1,nucl_num
+    ao_num, 1.d0, mo_coef, size(mo_coef,1), ao_pseudo_grid, &
-     do l=0,pseudo_lmax
+    size(ao_pseudo_grid,1), 0.d0, mo_pseudo_grid, size(mo_pseudo_grid,1))
-       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)
 
 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