diff --git a/devel/sr_correction/alpha.irp.f b/devel/sr_correction/alpha.irp.f index 2199697..3ce05c3 100644 --- a/devel/sr_correction/alpha.irp.f +++ b/devel/sr_correction/alpha.irp.f @@ -16,70 +16,79 @@ BEGIN_PROVIDER [ double precision, alpha_coef, (0:1) ] END_PROVIDER -BEGIN_PROVIDER [ double precision, alpha_coef_r, (0:1) ] +BEGIN_PROVIDER [ double precision, alpha_coef_r, (0:1,N_states) ] implicit none BEGIN_DOC ! SavCar-JCP-23, corrected for small mu END_DOC + integer :: istate double precision :: num, den double precision :: a0, b0, delta_E, E - E = energy_mu + correction_alpha_1 + do istate = 1,N_states + delta_E = (energy_mu(istate) + correction_alpha_1(istate)) - energy_mu0(istate) + a0 = delta_E * correction_mu0(istate) + 2.d0/dsqrt(dacos(-1.d0))*(delta_E - correction_mu0(istate)) * mu_erf + b0 = correction_mu0(istate)*correction_mu0(istate) - delta_E = E - energy_mu0 + num = 0.319820d0 + mu_erf * (1.063846d0 + mu_erf) + den = 0.487806d0 + mu_erf * (1.375439d0 + mu_erf) + alpha_coef_r(0,istate) = (a0 + mu_erf*num)/(b0 + mu_erf*den) - a0 = delta_E * correction_mu0 + 2.d0/dsqrt(dacos(-1.d0))*(delta_E - correction_mu0) * mu_erf - b0 = correction_mu0*correction_mu0 - num = 0.319820d0 + mu_erf * (1.063846d0 + mu_erf) - den = 0.487806d0 + mu_erf * (1.375439d0 + mu_erf) - alpha_coef_r(0) = (a0 + mu_erf*num)/(b0 + mu_erf*den) - - num = 0.113074d0 + mu_erf * (0.638308d0 + mu_erf) - den = 0.122652d0 + mu_erf * (0.674813d0 + mu_erf) - alpha_coef_r(1) = (a0 + mu_erf*num)/(b0 + mu_erf*den) + ! delta_E = (energy_mu + correction_alpha_1) - energy_mu0 + ! a0 = delta_E * correction_mu0 + 2.d0/dsqrt(dacos(-1.d0))*(delta_E - correction_mu0) * mu_erf + ! b0 = correction_mu0*correction_mu0 + num = 0.113074d0 + mu_erf * (0.638308d0 + mu_erf) + den = 0.122652d0 + mu_erf * (0.674813d0 + mu_erf) + alpha_coef_r(1,istate) = (a0 + mu_erf*num)/(b0 + mu_erf*den) + enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, energy_mu0 ] +BEGIN_PROVIDER [ double precision, energy_mu0, (N_states) ] implicit none BEGIN_DOC ! E(mu=0) END_DOC - integer :: l + integer :: l, istate double precision, external :: ddot - energy_mu0 = 0.d0 - do l=1,mo_num - energy_mu0 += ddot(mo_num, one_e_dm_mo(1,l), 1, mo_one_e_integrals(1,l), 1) + do istate = 1,N_states + energy_mu0(istate) = 0.d0 + do l=1,mo_num + energy_mu0(istate) += ddot(mo_num, one_e_dm_mo_alpha(1,l,istate), 1, mo_one_e_integrals(1,l), 1) + energy_mu0(istate) += ddot(mo_num, one_e_dm_mo_beta(1,l,istate), 1, mo_one_e_integrals(1,l), 1) + enddo + energy_mu0(istate) = energy_mu0(istate) + nuclear_repulsion enddo - energy_mu0 = energy_mu0 + nuclear_repulsion END_PROVIDER -BEGIN_PROVIDER [ double precision, correction_alpha_0_r ] +BEGIN_PROVIDER [ double precision, correction_alpha_0_r, (N_states) ] implicit none BEGIN_DOC ! alpha_0r(mu) * END_DOC - integer :: k,l + integer :: k,l, istate double precision :: c0 double precision, external :: ddot - c0 = 0.5d0 * alpha_coef_r(0) - correction_alpha_0_r = 0.d0 - do l=1,mo_num - do k=1,mo_num - correction_alpha_0_r += c0 * ddot (mo_num*mo_num, two_e_dm_mo(1,1,k,l), 1, W_bar_mu(1,1,k,l), 1) + do istate = 1,N_states + c0 = 0.5d0 * alpha_coef_r(0,istate) + 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) + enddo enddo enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, correction_alpha_1_r ] +BEGIN_PROVIDER [ double precision, correction_alpha_1_r, (N_states) ] implicit none BEGIN_DOC ! alpha_0(mu) * _s + alpha_1(mu) * _t @@ -87,17 +96,19 @@ BEGIN_PROVIDER [ double precision, correction_alpha_1_r ] double precision :: c1, c0 - integer :: k,l + integer :: k,l, istate double precision, external :: ddot - c0 = 0.5d0 * alpha_coef_r(0) - c1 = 0.5d0 * alpha_coef_r(1) + do istate = 1,N_states + c0 = 0.5d0 * alpha_coef_r(0,istate) + c1 = 0.5d0 * alpha_coef_r(1,istate) - correction_alpha_1_r = 0.d0 - do l=1,mo_num - do k=1,mo_num - correction_alpha_1_r += 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_r += c1 * ddot (mo_num*mo_num, two_e_dm_mo_triplet(1,1,k,l), 1, W_bar_mu(1,1,k,l), 1) + correction_alpha_1_r(istate) = 0.d0 + do l=1,mo_num + do k=1,mo_num + correction_alpha_1_r(istate) += c0 * ddot (mo_num*mo_num, two_e_dm_mo_singlet(1,1,k,l,istate), 1, W_bar_mu(1,1,k,l), 1) + correction_alpha_1_r(istate) += c1 * ddot (mo_num*mo_num, two_e_dm_mo_triplet(1,1,k,l,istate), 1, W_bar_mu(1,1,k,l), 1) + enddo enddo enddo END_PROVIDER diff --git a/devel/sr_correction/density_matrices.irp.f b/devel/sr_correction/density_matrices.irp.f index 82daa1f..21672d7 100644 --- a/devel/sr_correction/density_matrices.irp.f +++ b/devel/sr_correction/density_matrices.irp.f @@ -1,33 +1,37 @@ -BEGIN_PROVIDER [ double precision, two_e_dm_mo_singlet, (mo_num, mo_num, mo_num, mo_num) ] +BEGIN_PROVIDER [ double precision, two_e_dm_mo_singlet, (mo_num, mo_num, mo_num, mo_num, N_states) ] 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 + integer :: i,j,k,l,istate - 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)) + do istate = 1,N_states + 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,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)) + enddo enddo enddo enddo enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, two_e_dm_mo_triplet, (mo_num, mo_num, mo_num, mo_num) ] +BEGIN_PROVIDER [ double precision, two_e_dm_mo_triplet, (mo_num, mo_num, mo_num, mo_num, N_states) ] 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 + integer :: i,j,k,l, istate - 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)) + do istate = 1,N_states + 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,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)) + enddo enddo enddo enddo diff --git a/devel/sr_correction/energy_mu.irp.f b/devel/sr_correction/energy_mu.irp.f index da269f2..878db6d 100644 --- a/devel/sr_correction/energy_mu.irp.f +++ b/devel/sr_correction/energy_mu.irp.f @@ -55,7 +55,7 @@ BEGIN_PROVIDER [ double precision, W_bar_mu, (mo_num, mo_num, mo_num, mo_num) ] enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, energy_mu ] +BEGIN_PROVIDER [ double precision, energy_mu, (N_states) ] implicit none BEGIN_DOC ! E(mu) @@ -63,40 +63,45 @@ BEGIN_PROVIDER [ double precision, energy_mu ] double precision :: one_e, two_e - integer :: k,l + integer :: k,l, istate double precision, external :: ddot - one_e = 0.d0 - two_e = 0.d0 - do l=1,mo_num - one_e += ddot(mo_num, one_e_dm_mo(1,l), 1, mo_one_e_integrals(1,l), 1) - do k=1,mo_num - two_e += 0.5d0 * ddot (mo_num*mo_num, two_e_dm_mo(1,1,k,l), 1, W_mu(1,1,k,l), 1) + do istate = 1,N_states + one_e = 0.d0 + two_e = 0.d0 + do l=1,mo_num + 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) + enddo + enddo + energy_mu(istate) = one_e + two_e + nuclear_repulsion enddo - enddo - energy_mu = one_e + two_e + nuclear_repulsion END_PROVIDER -BEGIN_PROVIDER [ double precision, correction_alpha_0 ] +BEGIN_PROVIDER [ double precision, correction_alpha_0, (N_states) ] implicit none BEGIN_DOC ! alpha_0(mu) * END_DOC - integer :: k,l + integer :: k,l, istate 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) + do istate=1,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) + enddo enddo enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, correction_alpha_1 ] +BEGIN_PROVIDER [ double precision, correction_alpha_1, (N_states) ] implicit none BEGIN_DOC ! alpha_0(mu) * _s + alpha_1(mu) * _t @@ -104,55 +109,62 @@ BEGIN_PROVIDER [ double precision, correction_alpha_1 ] double precision :: c1, c0 - integer :: k,l + integer :: k,l, istate 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) + do istate=1,N_states + correction_alpha_1(istate) = 0.d0 + do l=1,mo_num + do k=1,mo_num + correction_alpha_1(istate) += c0 * ddot (mo_num*mo_num, two_e_dm_mo_singlet(1,1,k,l,istate), 1, W_bar_mu(1,1,k,l), 1) + correction_alpha_1(istate) += c1 * ddot (mo_num*mo_num, two_e_dm_mo_triplet(1,1,k,l,istate), 1, W_bar_mu(1,1,k,l), 1) + enddo enddo enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, correction_mu ] +BEGIN_PROVIDER [ double precision, correction_mu, (N_states) ] implicit none BEGIN_DOC ! END_DOC - integer :: k,l + integer :: k,l, istate 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) + + do istate=1,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) + enddo enddo enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, correction_mu0 ] +BEGIN_PROVIDER [ double precision, correction_mu0, (N_states) ] implicit none BEGIN_DOC ! : Should be END_DOC - integer :: k,l + integer :: k,l,istate double precision :: c0 double precision, external :: ddot c0 = 0.5d0 - correction_mu0 = 0.d0 - do l=1,mo_num - do k=1,mo_num - correction_mu0 += c0 * ddot (mo_num*mo_num, two_e_dm_mo(1,1,k,l), 1, W_bar_mu0(1,1,k,l), 1) + do istate=1,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) + enddo enddo enddo END_PROVIDER diff --git a/devel/sr_correction/sr_correction.irp.f b/devel/sr_correction/sr_correction.irp.f index 667d237..7b05222 100644 --- a/devel/sr_correction/sr_correction.irp.f +++ b/devel/sr_correction/sr_correction.irp.f @@ -8,27 +8,40 @@ end subroutine run implicit none + integer :: istate + do istate=1,N_states + print *, '', istate + print *, 'State ', istate print *, '---' print *, 'mu', mu_erf - print *, 'E(mu)', energy_mu + print *, 'E(mu)', energy_mu(istate) print *, '---' - print *, 'W_bar(mu)', correction_mu - print *, 'E(mu) + ', energy_mu + correction_mu + print *, 'W_bar(mu)', correction_mu(istate) + print *, 'E(mu) + ', energy_mu(istate) + correction_mu(istate) print *, '---' print *, 'alpha_0', alpha_coef(0) - print *, 'correction 0', correction_alpha_0 - print *, 'E(mu) + alpha_0 = ', energy_mu + correction_alpha_0 + print *, 'correction 0', correction_alpha_0(istate) + print *, 'E(mu) + alpha_0 = ', energy_mu(istate) + correction_alpha_0(istate) print *, '---' print *, 'alpha_1', alpha_coef(1) - print *, 'correction 1', correction_alpha_1 - print *, 'E(mu) + alpha_0 _s + alpha_1 _t = ', energy_mu + correction_alpha_1 + print *, 'correction 1', correction_alpha_1(istate) + print *, 'E(mu) + alpha_0 _s + alpha_1 _t = ', energy_mu(istate) + correction_alpha_1(istate) print *, '---' - print *, 'alpha_0_r', alpha_coef_r(0) - print *, 'correction 0', correction_alpha_0_r - print *, 'E(mu) + alpha_0_r = ', energy_mu + correction_alpha_0_r + print *, 'alpha_0_r', alpha_coef_r(0,istate) + print *, 'correction 0', correction_alpha_0_r(istate) + print *, 'E(mu) + alpha_0_r = ', energy_mu(istate) + correction_alpha_0_r(istate) print *, '---' - print *, 'alpha_1_r', alpha_coef_r(1) - print *, 'correction 1', correction_alpha_1_r - print *, 'E(mu) + alpha_0_r _s + alpha_1_r _t = ', energy_mu + correction_alpha_1_r + print *, 'alpha_1_r', alpha_coef_r(1,istate) + print *, 'correction 1', correction_alpha_1_r(istate) + print *, 'E(mu) + alpha_0_r _s + alpha_1_r _t = ', energy_mu(istate) + correction_alpha_1_r(istate) print *, '---' + + + print *,'' + print '(''|'',A6,''|'',5(A20,''|''))', 'E(mu)', '', 'E(mu) + ', 'E(mu) + \alpha_0', & + 'E(mu) + \alpha_0_s + \alpha_1_t', 'E(mu) + \alpha_0_r_s + \alpha_1_r_t' + print '(''|'',F6.1,''|'',5(F20.15,''|''))', mu_erf, energy_mu(istate), energy_mu(istate) + correction_mu(istate), energy_mu(istate) + & + correction_alpha_0(istate), energy_mu(istate) + correction_alpha_1(istate), energy_mu(istate) + & + correction_alpha_1_r(istate) + enddo end