quantum_package/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f

173 lines
5.6 KiB
Fortran

BEGIN_PROVIDER [ double precision, ao_ortho_mono_elec_integral_dressing, (ao_num,ao_num) ]
implicit none
BEGIN_DOC
! Dressing of the core hamiltonian in the orthogonal AO basis set
END_DOC
integer :: i,j,k
integer :: mu, nu, lambda, A
double precision :: tmp
ao_ortho_mono_elec_integral_dressing = 0.d0
i = idx_dressing
do mu=1,ao_num
if (dabs(mo_coef_in_ao_ortho_basis(mu,i)) > 1.d-5) then
do A=1,nucl_num
tmp = 0.d0
do nu=1,ao_num
tmp += AO_orthoSlaOverlap_matrix(nu,A) * ao_ortho_mono_elec_integral(mu,nu)
enddo
ao_ortho_mono_elec_integral_dressing(mu,mu) += cusp_C(A,i) * (AO_orthoSlaH_matrix(mu,A) - tmp)
enddo
ao_ortho_mono_elec_integral_dressing(mu,mu) *= 1.d0/mo_coef_in_ao_ortho_basis(mu,i)
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_ortho_mono_elec_integral, (ao_num, ao_num) ]
implicit none
BEGIN_DOC
! h core in orthogonal AO basis
END_DOC
double precision, allocatable :: T(:,:)
allocate(T(ao_num,ao_num))
call dgemm('T','N',ao_num,ao_num,ao_num,1.d0, &
ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), &
ao_mono_elec_integral, size(ao_mono_elec_integral,1), &
0.d0, T, size(T,1))
call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, &
T, size(T,1), &
ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), &
0.d0, ao_ortho_mono_elec_integral, size(ao_ortho_mono_elec_integral,1))
deallocate(T)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_mono_elec_integral_dressing, (ao_num,ao_num) ]
implicit none
BEGIN_DOC
! Dressing of the core hamiltonian in the AO basis set
END_DOC
call ao_ortho_cano_to_ao(ao_ortho_mono_elec_integral_dressing,size(ao_ortho_mono_elec_integral_dressing,1),&
ao_mono_elec_integral_dressing,size(ao_mono_elec_integral_dressing,1))
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_mono_elec_integral_dressing, (mo_tot_num,mo_tot_num) ]
implicit none
BEGIN_DOC
! Dressing of the core hamiltonian in the MO basis set
END_DOC
call ao_to_mo(ao_mono_elec_integral_dressing,size(ao_mono_elec_integral_dressing,1),&
mo_mono_elec_integral_dressing,size(mo_mono_elec_integral_dressing,1))
END_PROVIDER
BEGIN_PROVIDER [ integer, idx_dressing ]
implicit none
BEGIN_DOC
! Index of the MO which is being dressed
END_DOC
idx_dressing = 1
END_PROVIDER
BEGIN_PROVIDER [ double precision, cusp_corrected_mos, (ao_num,mo_tot_num) ]
implicit none
BEGIN_DOC
! Dressing core hamiltonian in the AO basis set
END_DOC
integer :: i,j
double precision, allocatable :: F(:,:), M(:,:)
allocate(F(mo_tot_num,mo_tot_num),M(ao_num,mo_tot_num))
logical :: oneshot
! oneshot = .True.
oneshot = .False.
if (oneshot) then
cusp_corrected_mos(1:ao_num,1:mo_tot_num) = mo_coef(1:ao_num,1:mo_tot_num)
slater_coef(1:nucl_num,1:mo_tot_num) = cusp_C(1:nucl_num,1:mo_tot_num)
return
else
do idx_dressing=1,mo_tot_num
if (idx_dressing>1) then
TOUCH idx_dressing
endif
do j=1,mo_tot_num
do i=1,mo_tot_num
F(i,j) = Fock_matrix_mo(i,j)
enddo
enddo
do j=1,mo_tot_num
do i=1,ao_num
M(i,j) = mo_coef(i,j)
enddo
enddo
integer :: it
do it=1,128
! print *, 'X', ao_ortho_canonical_coef(1:ao_num,1:ao_num)
! print *, 'C', mo_coef(1:ao_num,1:mo_tot_num)
! print *, 'Cp', mo_coef_in_ao_ortho_basis(1:ao_num,1:mo_tot_num)
! print *, 'cAi', cusp_C(1:nucl_num,1:mo_tot_num)
! print *, 'FmuA', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num)
! print *, 'Fock:', Fock_matrix_ao(1:ao_num,1:ao_num)
! print *, 'Diag Dressing:', ao_ortho_mono_elec_integral_dressing(1:ao_num,1:ao_num)
! print *, 'Dressing:', ao_mono_elec_integral_dressing(1:ao_num,1:ao_num)
! print *, 'Dressed Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) + ao_mono_elec_integral_dressing(1:ao_num,1:ao_num)
! print *, 'AO_orthoSlaOverlap_matrix', AO_orthoSlaOverlap_matrix(1:ao_num,1:nucl_num)
! print *, 'AO_orthoSlaH_matrix', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num)
! print *, 'ao_ortho_mono_elec_integral', ao_ortho_mono_elec_integral(1:ao_num,1:ao_num)
! print *, 'Fock MO:', Fock_matrix_mo(1:mo_tot_num,1:mo_tot_num)
do j=1,mo_tot_num
do i=1,mo_tot_num
Fock_matrix_mo(i,j) += mo_mono_elec_integral_dressing(i,j)
enddo
enddo
do i=1,mo_tot_num
Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i)
enddo
! print *, 'Dressed Fock MO:', Fock_matrix_mo(1:mo_tot_num,1:mo_tot_num)
double precision :: conv
conv = 0.d0
do j=1,mo_tot_num
do i=1,mo_tot_num
if (i==j) cycle
conv = max(conv,Fock_matrix_mo(i,j))
enddo
enddo
TOUCH Fock_matrix_mo Fock_matrix_diag_mo
mo_coef(1:ao_num,1:mo_tot_num) = eigenvectors_fock_matrix_mo(1:ao_num,1:mo_tot_num)
TOUCH mo_coef
!print *, 'C', mo_coef(1:ao_num,1:mo_tot_num)
!print *, '-----'
print *, idx_dressing, it, real(mo_coef(1,idx_dressing)), real(conv)
if (conv < 1.d-5) exit
!stop
enddo
cusp_corrected_mos(1:ao_num,idx_dressing) = mo_coef(1:ao_num,idx_dressing)
slater_coef(1:nucl_num,idx_dressing) = cusp_C(1:nucl_num,idx_dressing)
enddo
idx_dressing = 1
mo_coef(1:ao_num,1:mo_tot_num) = M(1:ao_num,1:mo_tot_num)
soft_TOUCH mo_coef idx_dressing slater_coef
endif
END_PROVIDER