devel/sr_correction/density_matrices.irp.f

@@ -1,36 +0,0 @@
-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,57 +64,17 @@ 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)
+     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)
    enddo
  enddo
+ correction_alpha_0 = correction_alpha_0*alpha_coef(0)
 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
-
- 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_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
-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,19 +8,7 @@ end
 
 subroutine run
  implicit none
- print *, '---'
- print *, 'mu', mu_erf
  print *, 'E(mu)', energy_mu
- print *, '---'
+ print *, 'correction', correction_alpha_0
- print *, 'W_bar(mu)', correction_mu
+ print *, 'E = ', energy_mu + correction_alpha_0
- 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