LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-26 07:02:14 +02:00
Code Issues Releases Wiki Activity

Working on Slater dressing

Browse Source
This commit is contained in:
Anthony Scemama 2017-06-08 22:15:42 +02:00
parent 9242ca4584
commit 4cb4d5e416
10 changed files with 364 additions and 138 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
config/gfortran_debug.cfg
View File

@ -13,7 +13,7 @@
FC : gfortran -g -ffree-line-length-none -I .
LAPACK_LIB : -llapack -lblas
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32
IRPF90_FLAGS : --ninja --align=32 --assert
# Global options
################

19
plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f
View File

@ -216,8 +216,25 @@ END_DOC
double precision, allocatable :: AF(:,:)
allocate (AF(dim_DIIS+1,dim_DIIS+1))
double precision :: rcond, ferr, berr
integer :: iwork(dim_DIIS+1)
integer :: iwork(dim_DIIS+1), lwork
call dsysvx('N','U',dim_DIIS+1,1, &
B_matrix_DIIS,size(B_matrix_DIIS,1), &
AF, size(AF,1), &
ipiv, &
C_vector_DIIS,size(C_vector_DIIS,1), &
X_vector_DIIS,size(X_vector_DIIS,1), &
rcond, &
ferr, &
berr, &
scratch,-1, &
iwork, &
info &
)
lwork = int(scratch(1,1))
deallocate(scratch)
allocate(scratch(lwork,1))
call dsysvx('N','U',dim_DIIS+1,1, &
B_matrix_DIIS,size(B_matrix_DIIS,1), &
AF, size(AF,1), &

172
plugins/Hartree_Fock/diagonalize_fock.irp.f
View File

@ -1,115 +1,129 @@
BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (mo_tot_num) ]
BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (ao_num) ]
&BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Diagonal Fock matrix in the MO basis
END_DOC
integer :: i,j, m
integer :: i,j
integer :: liwork, lwork, n, info
integer, allocatable :: iwork(:), isuppz(:)
double precision, allocatable :: work(:), F(:,:), F2(:,:)
integer :: iorb,jorb
integer, allocatable :: iwork(:)
double precision, allocatable :: work(:), F(:,:), S(:,:)
allocate( F(mo_tot_num,mo_tot_num),F2(mo_tot_num,mo_tot_num), isuppz(2*mo_tot_num) )
allocate( F(mo_tot_num,mo_tot_num) )
do j=1,mo_tot_num
do i=1,mo_tot_num
F(i,j) = Fock_matrix_mo(i,j)
enddo
enddo
if(no_oa_or_av_opt)then
integer :: iorb,jorb
do i = 1, n_act_orb
iorb = list_act(i)
ASSERT (iorb > 0)
ASSERT (iorb <= mo_tot_num)
do j = 1, n_inact_orb
jorb = list_inact(j)
ASSERT (jorb > 0)
ASSERT (jorb <= mo_tot_num)
F(iorb,jorb) = 0.d0
F(jorb,iorb) = 0.d0
enddo
do j = 1, n_virt_orb
jorb = list_virt(j)
ASSERT (jorb > 0)
ASSERT (jorb <= mo_tot_num)
F(iorb,jorb) = 0.d0
F(jorb,iorb) = 0.d0
enddo
do j = 1, n_core_orb
jorb = list_core(j)
ASSERT (jorb > 0)
ASSERT (jorb <= mo_tot_num)
F(iorb,jorb) = 0.d0
F(jorb,iorb) = 0.d0
enddo
enddo
endif
! Insert level shift here
do i = elec_beta_num+1, elec_alpha_num
F(i,i) += 0.5d0*level_shift
enddo
do i = elec_alpha_num+1, mo_tot_num
F(i,i) += level_shift
enddo
n = mo_tot_num
lwork = 1+6*n + 2*n*n
liwork = 10*n
allocate(work(lwork))
allocate(iwork(liwork) )
call dsyevr('V', 'A', 'U', mo_tot_num, F, size(F,1), &
-100.d0, 100.d0, 1, mo_tot_num, 0.d0, &
m, diagonal_Fock_matrix_mo, &
F2, size(F2,1), &
isuppz, work, lwork, iwork, liwork, info)
if (info /= 0) then
print *, irp_here//' DSYEV failed : ', info
stop 1
endif
call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, &
mo_coef, size(mo_coef,1), F2, size(F2,1), &
0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1))
deallocate(work, F2, F)
deallocate(iwork, isuppz)
END_PROVIDER
endif
BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)]
implicit none
BEGIN_DOC
! diagonal element of the fock matrix calculated as the sum over all the interactions
! with all the electrons in the RHF determinant
! diagonal_Fock_matrix_mo_sum(i) = sum_{j=1, N_elec} 2 J_ij -K_ij
END_DOC
integer :: i,j
double precision :: accu
do j = 1,elec_alpha_num
accu = 0.d0
do i = 1, elec_alpha_num
accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j)
enddo
diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j)
enddo
do j = elec_alpha_num+1,mo_tot_num
accu = 0.d0
do i = 1, elec_alpha_num
accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j)
enddo
diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j)
enddo
! Insert level shift here
do i = elec_beta_num+1, elec_alpha_num
F(i,i) += 0.5d0*level_shift
enddo
do i = elec_alpha_num+1, mo_tot_num
F(i,i) += level_shift
enddo
n = mo_tot_num
lwork = 1+6*n + 2*n*n
liwork = 3 + 5*n
allocate(work(lwork))
allocate(iwork(liwork) )
lwork = -1
liwork = -1
call dsyevd( 'V', 'U', mo_tot_num, F, &
size(F,1), diagonal_Fock_matrix_mo, &
work, lwork, iwork, liwork, info)
if (info /= 0) then
print *, irp_here//' DSYEVD failed : ', info
stop 1
endif
lwork = int(work(1))
liwork = iwork(1)
deallocate(iwork)
deallocate(work)
allocate(work(lwork))
allocate(iwork(liwork) )
call dsyevd( 'V', 'U', mo_tot_num, F, &
size(F,1), diagonal_Fock_matrix_mo, &
work, lwork, iwork, liwork, info)
deallocate(iwork)
if (info /= 0) then
call dsyev( 'V', 'L', mo_tot_num, F, &
size(F,1), diagonal_Fock_matrix_mo, &
work, lwork, info)
if (info /= 0) then
print *, irp_here//' DSYEV failed : ', info
stop 1
endif
endif
call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, &
mo_coef, size(mo_coef,1), F, size(F,1), &
0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1))
deallocate(work, F)
END_PROVIDER
BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)]
implicit none
BEGIN_DOC
! diagonal element of the fock matrix calculated as the sum over all the interactions
! with all the electrons in the RHF determinant
! diagonal_Fock_matrix_mo_sum(i) = sum_{j=1, N_elec} 2 J_ij -K_ij
END_DOC
integer :: i,j
double precision :: accu
do j = 1,elec_alpha_num
accu = 0.d0
do i = 1, elec_alpha_num
accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j)
enddo
diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j)
enddo
do j = elec_alpha_num+1,mo_tot_num
accu = 0.d0
do i = 1, elec_alpha_num
accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j)
enddo
diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j)
enddo
END_PROVIDER

14
plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg Normal file
View File

@ -0,0 +1,14 @@
[slater_expo_ezfio]
type: double precision
doc: Exponents of the additional Slater functions
size: (nuclei.nucl_num)
interface: ezfio, provider
[slater_coef_ezfio]
type: double precision
doc: Exponents of the additional Slater functions
size: (mo_basis.mo_tot_num,nuclei.nucl_num)
interface: ezfio, provider

53
plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f
View File

@ -6,21 +6,19 @@ BEGIN_PROVIDER [ double precision, cusp_A, (nucl_num, nucl_num) ]
integer :: mu, A, B
do B=1,nucl_num
do A=1,nucl_num
cusp_A(A,B) = 0.d0
if (A/=B) then
cusp_A(A,B) -= slater_value_at_nucl(A,B)
endif
do mu=1,ao_num
cusp_A(A,B) += slater_overlap(mu,B) * ao_value_at_nucl(mu,A)
enddo
enddo
cusp_A = 0.d0
do A=1,nucl_num
cusp_A(A,A) = slater_expo(A)/nucl_charge(A) * slater_value_at_nucl(A,A)
do B=1,nucl_num
cusp_A(A,B) -= slater_value_at_nucl(B,A)
do mu=1,ao_num
cusp_A(A,B) += GauSlaOverlap_matrix(mu,B) * ao_value_at_nucl(mu,A)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, cusp_C, (nucl_num, mo_tot_num) ]
BEGIN_PROVIDER [ double precision, cusp_B, (nucl_num, mo_tot_num) ]
implicit none
BEGIN_DOC
! Equations to solve : A.C = B
@ -30,20 +28,25 @@ BEGIN_PROVIDER [ double precision, cusp_C, (nucl_num, mo_tot_num) ]
do i=1,mo_tot_num
do A=1,nucl_num
cusp_C(A,i) = mo_value_at_nucl(i,A)
cusp_B(A,i) = mo_value_at_nucl(i,A)
enddo
enddo
integer, allocatable :: ipiv(:)
allocate ( ipiv(nucl_num) )
call dgegv(nucl_num, mo_tot_num, cusp_A, size(cusp_A,1), &
ipiv, cusp_C, size(cusp_C,1), info)
deallocate (ipiv)
if (info /= 0) then
print *, 'Cusp : linear solve failed'
stop -1
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, cusp_C, (nucl_num, mo_tot_num) ]
implicit none
BEGIN_DOC
! Equations to solve : A.C = B
END_DOC
double precision, allocatable :: AF(:,:)
integer :: info
allocate ( AF(nucl_num,nucl_num) )
call get_pseudo_inverse(cusp_A,nucl_num,nucl_num,AF,size(AF,1))
call dgemm('N','N',nucl_num,mo_tot_num,nucl_num,1.d0, &
AF,size(AF,1), cusp_B, size(cusp_B,1), 0.d0, cusp_C, size(cusp_C,1))
END_PROVIDER

69
plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f Normal file
View File

@ -0,0 +1,69 @@
program scf
BEGIN_DOC
! Produce `Hartree_Fock` MO orbital with Slater cusp dressing
! output: mo_basis.mo_tot_num mo_basis.mo_label mo_basis.ao_md5 mo_basis.mo_coef mo_basis.mo_occ
! output: hartree_fock.energy
! optional: mo_basis.mo_coef
END_DOC
call check_mos
call debug
call run
end
subroutine check_mos
implicit none
BEGIN_DOC
! Create a MO guess if no MOs are present in the EZFIO directory
END_DOC
logical :: exists
PROVIDE ezfio_filename
call ezfio_has_mo_basis_mo_coef(exists)
if (.not.exists) then
print *, 'Please run SCF first'
stop
endif
end
subroutine debug
implicit none
integer :: i
print *, 'A'
do i=1,nucl_num
print *, i, cusp_A(1:nucl_num, i)
enddo
print *, 'B'
do i=1,mo_tot_num
print *, i, cusp_B(1:nucl_num, i)
enddo
print *, 'C'
do i=1,mo_tot_num
print *, i, cusp_C(1:nucl_num, i)
enddo
end
subroutine run
BEGIN_DOC
! Run SCF calculation
END_DOC
use bitmasks
implicit none
double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem
double precision :: EHF
integer :: i_it, i, j, k
EHF = HF_energy
mo_label = "CuspDressed"
call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(cusp_B)
! Choose SCF algorithm
! call Roothaan_Hall_SCF
end

2
plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f
View File

@ -58,7 +58,7 @@ BEGIN_PROVIDER [ double precision , slater_value_at_nucl, (nucl_num,nucl_num) ]
expo = slater_expo(i)*slater_expo(i)*((x*x) + (y*y) + (z*z))
if (expo > 160.d0) cycle
expo = dsqrt(expo)
slater_value_at_nucl(i,k) = dexp(-expo)
slater_value_at_nucl(i,k) = dexp(-expo) * slater_normalization(i)
enddo
enddo
END_PROVIDER

101
plugins/Hartree_Fock_SlaterDressed/integrals.irp.f
View File

@ -1,15 +1,17 @@
!*****************************************************************************
subroutine GauSlaOverlap(expGau,cGau,aGau,expSla,cSla)
! Compute the overlap integral between a Gaussian function
! with arbitrary angular momemtum and a s-type Slater function
subroutine GauSlaOverlap(expGau,cGau,aGau,expSla,cSla,result)
implicit none
BEGIN_DOC
! Compute the overlap integral between a Gaussian function
! with arbitrary angular momemtum and a s-type Slater function
END_DOC
! Input variables
double precision,intent(in) :: expGau,expSla
double precision,intent(in) :: cGau(3),cSla(3)
integer,intent(in) :: aGau(3)
double precision,intent(out) :: result
! Final value of the integrals
double precision :: ss,ps,ds
@ -82,13 +84,38 @@ subroutine GauSlaOverlap(expGau,cGau,aGau,expSla,cSla)
dxzs = AxBx*AzBz*ds
dyzs = AyBy*AzBz*ds
! Print result
write(*,'(A10,F16.10)') &
'(s|s) = ',ss
write(*,'(A10,F16.10,3X,A10,F16.10,3X,A10,F16.10)') &
'(px|s) = ',pxs,'(py|s) = ',pys,'(pz|s) = ',pzs
write(*,'(A10,F16.10,3X,A10,F16.10,3X,A10,F16.10,3X,A10,F16.10,3X,A10,F16.10,3X,A10,F16.10)') &
'(dx2|s) = ',dxxs,'(dy2|s) = ',dyys,'(dz2|s) = ',dzzs,'(dxy|s) = ',dxys,'(dxz|s) = ',dxzs,'(dyz|s) = ',dyzs
select case (sum(aGau))
case (0)
result = ss
case (1)
if (aGau(1) == 1) then
result = pxs
else if (aGau(2) == 1) then
result = pys
else if (aGau(3) == 1) then
result = pzs
endif
case (2)
if (aGau(1) == 2) then
result = dxxs
else if (aGau(2) == 2) then
result = dyys
else if (aGau(3) == 2) then
result = dzzs
else if (aGau(1)+aGau(2) == 2) then
result = dxys
else if (aGau(1)+aGau(3) == 2) then
result = dxzs
else if (aGau(2)+aGau(3) == 2) then
result = dyzs
endif
case default
stop 'GauSlaOverlap not implemented'
end select
end
!*****************************************************************************
@ -97,11 +124,13 @@ end
!*****************************************************************************
subroutine GauSlaKinetic(expGau,cGau,aGau,expSla,cSla)
! Compute the kinetic energy integral between a Gaussian function
! with arbitrary angular momemtum and a s-type Slater function
implicit none
BEGIN_DOC
! Compute the kinetic energy integral between a Gaussian function
! with arbitrary angular momemtum and a s-type Slater function
END_DOC
! Input variables
double precision,intent(in) :: expGau,expSla
double precision,intent(in) :: cGau(3),cSla(3)
@ -195,11 +224,13 @@ end
!*****************************************************************************
subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc)
! Compute the nuclear attraction integral between a Gaussian function
! with arbitrary angular momemtum and a s-type Slater function
implicit none
BEGIN_DOC
! Compute the nuclear attraction integral between a Gaussian function
! with arbitrary angular momemtum and a s-type Slater function
END_DOC
! Input variables
double precision,intent(in) :: expGau,expSla
double precision,intent(in) :: cGau(3),cSla(3)
@ -242,7 +273,8 @@ subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc)
end
!*****************************************************************************
double precision function BoysF0(t)
implicit none
double precision, intent(in) :: t
double precision :: pi
pi = 4d0*atan(1d0)
@ -257,4 +289,35 @@ end
!*****************************************************************************
BEGIN_PROVIDER [ double precision, GauSlaOverlap_matrix, (ao_num, nucl_num) ]
implicit none
BEGIN_DOC
! <Gaussian | Slater> overlap matrix
END_DOC
integer :: i,j,k
double precision :: cGau(3)
double precision :: cSla(3)
double precision :: expSla, res, expGau
integer :: aGau(3)
do k=1,nucl_num
cSla(1:3) = nucl_coord_transp(1:3,k)
expSla = slater_expo(k)
do i=1,ao_num
cGau(1:3) = nucl_coord_transp(1:3, ao_nucl(i))
aGau(1:3) = ao_power(i,1:3)
GauSlaOverlap_matrix(i,k) = 0.d0
do j=1,ao_prim_num(i)
expGau = ao_expo_ordered_transp(j,i)
call GauSlaOverlap(expGau,cGau,aGau,expSla,cSla,res)
GauSlaOverlap_matrix(i,k) += ao_coef_normalized_ordered_transp(j,i) * res
enddo
enddo
enddo
END_PROVIDER

43
plugins/Hartree_Fock_SlaterDressed/slater.irp.f Normal file
View File

@ -0,0 +1,43 @@
BEGIN_PROVIDER [ double precision, slater_expo, (nucl_num) ]
implicit none
BEGIN_DOC
! Exponents of the Slater functions
END_DOC
logical :: exists
call ezfio_has_Hartree_Fock_SlaterDressed_slater_expo_ezfio(exists)
if (exists) then
slater_expo(1:nucl_num) = slater_expo_ezfio(1:nucl_num)
else
slater_expo(1:nucl_num) = nucl_charge(1:nucl_num)
call ezfio_set_Hartree_Fock_SlaterDressed_slater_expo_ezfio(slater_expo)
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, slater_coef, (nucl_num,mo_tot_num) ]
implicit none
BEGIN_DOC
! Exponents of the Slater functions
END_DOC
logical :: exists
slater_coef = 0.d0
call ezfio_has_Hartree_Fock_SlaterDressed_slater_coef_ezfio(exists)
if (exists) then
slater_coef = slater_coef_ezfio
else
call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(slater_coef)
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, slater_normalization, (nucl_num) ]
implicit none
BEGIN_DOC
! Normalization of Slater functions : sqrt(expo^3/pi)
END_DOC
integer :: i
do i=1,nucl_num
slater_normalization(i) = dsqrt( slater_expo(i)**3/dacos(-1.d0) )
enddo
END_PROVIDER

27
plugins/mrcepa0/dressing.irp.f
View File

@ -252,25 +252,28 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
else if (perturbative_triples) then
! Linked
call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv)
hka = hij_cache(idx_alpha(k_sd))
do i_state=1,N_states
ASSERT (Delta_E_inv(i_state) < 0.d0)
dka(i_state) = hka / Delta_E_inv(i_state)
enddo
if (dabs(hka) > 1.d-12) then
call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv)
do i_state=1,N_states
ASSERT (Delta_E_inv(i_state) < 0.d0)
dka(i_state) = hka / Delta_E_inv(i_state)
enddo
endif
endif
if (perturbative_triples.and. (degree2 == 1) ) then
call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv)
call i_h_j(psi_ref(1,1,i_I),tmp_det,Nint,hka)
hka = hij_cache(idx_alpha(k_sd)) - hka
do i_state=1,N_states
ASSERT (Delta_E_inv(i_state) < 0.d0)
dka(i_state) = hka / Delta_E_inv(i_state)
enddo
if (dabs(hka) > 1.d-12) then
call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv)
do i_state=1,N_states
ASSERT (Delta_E_inv(i_state) < 0.d0)
dka(i_state) = hka / Delta_E_inv(i_state)
enddo
endif
endif