mirror of
https://gitlab.com/scemama/qp_plugins_scemama.git
synced 2024-06-02 11:25:23 +02:00
116 lines
3.2 KiB
Fortran
116 lines
3.2 KiB
Fortran
BEGIN_PROVIDER [ double precision, alpha_coef, (0:1) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! SavCar-JCP-23
|
|
END_DOC
|
|
|
|
double precision :: num, den
|
|
|
|
num = 0.319820d0 + mu_erf * (1.063846d0 + mu_erf)
|
|
den = 0.487806d0 + mu_erf * (1.375439d0 + mu_erf)
|
|
alpha_coef(0) = num/den
|
|
|
|
num = 0.113074d0 + mu_erf * (0.638308d0 + mu_erf)
|
|
den = 0.122652d0 + mu_erf * (0.674813d0 + mu_erf)
|
|
alpha_coef(1) = num/den
|
|
END_PROVIDER
|
|
|
|
|
|
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
|
|
|
|
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)
|
|
|
|
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)
|
|
|
|
|
|
! 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, (N_states) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! E(mu=0)
|
|
END_DOC
|
|
|
|
integer :: l, istate
|
|
double precision, external :: ddot
|
|
|
|
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
|
|
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, correction_alpha_0_r, (N_states) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! alpha_0r(mu) * <Psi(mu)|W_bar_mu|Psi(mu)>
|
|
END_DOC
|
|
|
|
integer :: k,l, istate
|
|
double precision :: c0
|
|
double precision, external :: ddot
|
|
|
|
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, (N_states) ]
|
|
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, istate
|
|
double precision, external :: ddot
|
|
|
|
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(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
|
|
|