From 7b7db6786f89d5179bc6df4c963c2ea4a77a52c6 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 17 Oct 2023 00:29:34 +0200
Subject: [PATCH] SR correction OK

---
 devel/sr_correction/energy_mu.irp.f     | 48 ++++++++++++++++++++++---
 devel/sr_correction/sr_correction.irp.f | 16 +++++++--
 2 files changed, 58 insertions(+), 6 deletions(-)

diff --git a/devel/sr_correction/energy_mu.irp.f b/devel/sr_correction/energy_mu.irp.f
index 0ee1533..7cf8b20 100644
--- a/devel/sr_correction/energy_mu.irp.f
+++ b/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
diff --git a/devel/sr_correction/sr_correction.irp.f b/devel/sr_correction/sr_correction.irp.f
index 76c9866..dda293c 100644
--- a/devel/sr_correction/sr_correction.irp.f
+++ b/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 *, 'E = ', energy_mu + correction_alpha_0
+ print *, '---'
+ 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