Forgot file

devel/sr_correction/density_matrices.irp.f

@@ -0,0 +1,36 @@
+BEGIN_PROVIDER [ double precision, two_e_dm_mo_singlet, (mo_num, mo_num, mo_num, mo_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Ps(r1,r2,r1',r2') = 1/2 (P(r1,r2,r1',r2') + P(r1,r2,r2',r1'))
+ END_DOC
+ integer :: i,j,k,l
+
+ do l=1,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) = 0.5d0 * (two_e_dm_mo(i,j,k,l) + two_e_dm_mo(i,j,l,k))
+    enddo
+   enddo
+  enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, two_e_dm_mo_triplet, (mo_num, mo_num, mo_num, mo_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Ps(r1,r2,r1',r2') = 1/2 (P(r1,r2,r1',r2') - P(r1,r2,r2',r1'))
+ END_DOC
+ integer :: i,j,k,l
+
+ do l=1,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) = 0.5d0 * (two_e_dm_mo(i,j,k,l) - two_e_dm_mo(i,j,l,k))
+    enddo
+   enddo
+  enddo
+ enddo
+END_PROVIDER
+

devel/sr_correction/energy_mu.irp.f

@@ -64,17 +64,57 @@ BEGIN_PROVIDER [ double precision, correction_alpha_0 ]
  ! alpha_0(mu) * <Psi(mu)|W_bar_mu|Psi(mu)>
  END_DOC
 
- double precision :: one_e, two_e
+ integer :: k,l
+ double precision :: c0
+ double precision, external :: ddot
+
+ c0 = 0.5d0 * alpha_coef(0)
+ correction_alpha_0 = 0.d0
+ do l=1,mo_num
+   do k=1,mo_num
+     correction_alpha_0 += c0 * ddot (mo_num*mo_num, two_e_dm_mo(1,1,k,l), 1, W_bar_mu(1,1,k,l), 1)
+   enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, correction_alpha_1 ]
+ implicit none
+ BEGIN_DOC
+ ! alpha_0(mu) * <Psi(mu)|W_bar_mu|Psi(mu)>_s + alpha_1(mu) * <Psi(mu)|W_bar_mu|Psi(mu)>_t
+ END_DOC
+
+ double precision :: c1, c0
 
  integer :: k,l
  double precision, external :: ddot
 
- correction_alpha_0 = 0.d0
+ c0 = 0.5d0 * alpha_coef(0)
+ c1 = 0.5d0 * alpha_coef(1)
+
+ correction_alpha_1 = 0.d0
  do l=1,mo_num
    do k=1,mo_num
-     correction_alpha_0 += 0.5d0 * ddot (mo_num*mo_num, two_e_dm_mo(1,1,k,l), 1, W_bar_mu(1,1,k,l), 1)
+     correction_alpha_1 += c0 * ddot (mo_num*mo_num, two_e_dm_mo_singlet(1,1,k,l), 1, W_bar_mu(1,1,k,l), 1)
+     correction_alpha_1 += c1 * ddot (mo_num*mo_num, two_e_dm_mo_triplet(1,1,k,l), 1, W_bar_mu(1,1,k,l), 1)
    enddo
  enddo
- correction_alpha_0 = correction_alpha_0*alpha_coef(0)
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, correction_mu ]
+ implicit none
+ BEGIN_DOC
+ ! alpha_0(mu) * <Psi(mu)|W_bar_mu|Psi(mu)>
+ END_DOC
+
+ integer :: k,l
+ double precision :: c0
+ double precision, external :: ddot
+
+ c0 = 0.5d0
+ correction_mu = 0.d0
+ do l=1,mo_num
+   do k=1,mo_num
+     correction_mu += c0 * ddot (mo_num*mo_num, two_e_dm_mo(1,1,k,l), 1, W_bar_mu(1,1,k,l), 1)
+   enddo
+ enddo
+END_PROVIDER

devel/sr_correction/sr_correction.irp.f

@@ -8,7 +8,19 @@ end
 
 subroutine run
  implicit none
+ print *, '---'
+ print *, 'mu', mu_erf
  print *, 'E(mu)', energy_mu
- print *, 'correction', correction_alpha_0
+ print *, '---'
- print *, 'E = ', energy_mu + correction_alpha_0
+ print *, 'W_bar(mu)', correction_mu
+ print *, 'E(mu) + <W_bar(mu)>', energy_mu + correction_mu
+ print *, '---'
+ print *, 'alpha_0', alpha_coef(0)
+ print *, 'correction 0', correction_alpha_0
+ print *, 'E(mu) + alpha_0 <W_bar(mu)> = ', energy_mu + correction_alpha_0
+ print *, '---'
+ print *, 'alpha_1', alpha_coef(1)
+ print *, 'correction 1', correction_alpha_1
+ print *, 'E(mu) + alpha_0 <W_bar(mu)>_s + alpha_1 <W_bar(mu)>_t = ', energy_mu + correction_alpha_1
+ print *, '---'
 end