LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-11-03 20:54:00 +01:00
Code Issues Releases Wiki Activity

Improved Hartree-Fock

Browse Source
This commit is contained in:
Anthony Scemama 2014-06-12 22:08:53 +02:00
parent 3beea8d230
commit 42d8b4c404
4 changed files with 66 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
src/AOs/aos.irp.f
View File

@ -93,46 +93,6 @@ END_PROVIDER
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_overlap, (ao_num_align,ao_num) ]
implicit none
BEGIN_DOC
! matrix of the overlap for tha AOs
! S(i,j) = overlap between the ith and the jth atomic basis function
END_DOC
integer :: i,j,k,l,nz,num_i,num_j,powA(3),powB(3)
double precision :: accu,overlap_x,overlap_y,overlap_z,A_center(3),B_center(3),norm
nz=100
do i = 1, ao_num
num_i = ao_nucl(i)
powA(1) = ao_power(i,1)
powA(2) = ao_power(i,2)
powA(3) = ao_power(i,3)
A_center(1)=nucl_coord(num_i,1)
A_center(2)=nucl_coord(num_i,2)
A_center(3)=nucl_coord(num_i,3)
do j = i, ao_num
num_j = ao_nucl(j)
powB(1) = ao_power(j,1)
powB(2) = ao_power(j,2)
powB(3) = ao_power(j,3)
B_center(1)=nucl_coord(num_j,1)
B_center(2)=nucl_coord(num_j,2)
B_center(3)=nucl_coord(num_j,3)
accu = 0.d0
do k = 1, ao_prim_num(i)
do l = 1, ao_prim_num(j)
call overlap_gaussian_xyz(A_center,B_center,ao_expo(i,k),ao_expo(j,l),powA,powB,overlap_x,overlap_y,overlap_z,norm,nz)
accu = accu + norm * ao_coef(i,k) * ao_coef(j,l)
enddo
enddo
ao_overlap(i,j) = accu
ao_overlap(j,i) = accu
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_coef_transp, (ao_prim_num_max_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_expo_transp, (ao_prim_num_max_align,ao_num) ]
implicit none

32
src/Hartree_Fock/HF_density_matrix_ao.irp.f
View File

@ -5,13 +5,13 @@
! Alpha and Beta density matrix in the AO basis
END_DOC
integer :: i,j,k,l1,l2
integer, allocatable :: mo_occ(:,:)
integer, allocatable :: occ(:,:)
allocate ( mo_occ(elec_alpha_num,2) )
call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int)
allocate ( occ(elec_alpha_num,2) )
call bitstring_to_list( HF_bitmask(1,1), occ(1,1), j, N_int)
ASSERT ( j==elec_alpha_num )
call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int)
call bitstring_to_list( HF_bitmask(1,2), occ(1,2), j, N_int)
ASSERT ( j==elec_beta_num )
do j=1,ao_num
@ -21,8 +21,8 @@
HF_density_matrix_ao_beta (i,j) = 0.d0
enddo
do k=1,elec_beta_num
l1 = mo_occ(k,1)
l2 = mo_occ(k,2)
l1 = occ(k,1)
l2 = occ(k,2)
!DIR$ VECTOR ALIGNED
do i=1,ao_num
HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +&
@ -32,7 +32,7 @@
enddo
enddo
do k=elec_beta_num+1,elec_alpha_num
l1 = mo_occ(k,1)
l1 = occ(k,1)
!DIR$ VECTOR ALIGNED
do i=1,ao_num
HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +&
@ -40,7 +40,7 @@
enddo
enddo
enddo
deallocate(mo_occ)
deallocate(occ)
END_PROVIDER
BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]
@ -49,13 +49,13 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]
! Density matrix in the AO basis
END_DOC
integer :: i,j,k,l1,l2
integer, allocatable :: mo_occ(:,:)
integer, allocatable :: occ(:,:)
allocate ( mo_occ(elec_alpha_num,2) )
call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int)
allocate ( occ(elec_alpha_num,2) )
call bitstring_to_list( HF_bitmask(1,1), occ(1,1), j, N_int)
ASSERT ( j==elec_alpha_num )
call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int)
call bitstring_to_list( HF_bitmask(1,2), occ(1,2), j, N_int)
ASSERT ( j==elec_beta_num )
do j=1,ao_num
@ -64,8 +64,8 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]
HF_density_matrix_ao(i,j) = 0.d0
enddo
do k=1,elec_beta_num
l1 = mo_occ(k,1)
l2 = mo_occ(k,2)
l1 = occ(k,1)
l2 = occ(k,2)
!DIR$ VECTOR ALIGNED
do i=1,ao_num
HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) + &
@ -74,7 +74,7 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]
enddo
enddo
do k=elec_beta_num+1,elec_alpha_num
l1 = mo_occ(k,1)
l1 = occ(k,1)
!DIR$ VECTOR ALIGNED
do i=1,ao_num
HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) + &
@ -82,6 +82,6 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]
enddo
enddo
enddo
deallocate(mo_occ)
deallocate(occ)
END_PROVIDER

47
src/Hartree_Fock/mo_SCF_iterations.irp.f
View File

@ -3,21 +3,55 @@ program scf_iteration
implicit none
double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral
double precision :: E0
integer :: i_it
integer :: i_it, i, j, k
integer, allocatable :: iorder(:)
double precision, allocatable :: DM_occ(:,:), E_new(:), R(:,:)
E0 = HF_energy
i_it = 0
n_it_scf_max = 10
SCF_energy_before = huge(1.d0)
i_it = 1
n_it_scf_max = 100
SCF_energy_before = 0.d0
SCF_energy_after = E0
print *, E0
mo_label = "Canonical"
thresh_SCF = 1.d-10
do while (i_it < 40 .and. dabs(SCF_energy_before - SCF_energy_after) > thresh_SCF)
SCF_energy_before = SCF_energy_after
DM_occ = mo_density_matrix
allocate (DM_occ(size(mo_density_matrix,1),mo_tot_num), &
E_new(mo_tot_num), R(mo_tot_num,mo_tot_num), iorder(mo_tot_num))
do while (i_it < n_it_scf_max .and. dabs(SCF_energy_before - SCF_energy_after) > thresh_SCF)
if (SCF_energy_after <= SCF_energy_before+1.d-4) then
mo_coef = eigenvectors_Fock_matrix_mo
TOUCH mo_coef mo_label
DM_occ = mo_density_matrix
else
DM_occ = mo_density_matrix
mo_coef = eigenvectors_Fock_matrix_mo
TOUCH mo_coef mo_label mo_integrals_map
DM_occ = DM_occ + 0.0d0*mo_density_matrix
integer :: rank
call cholesky_mo(ao_num,mo_tot_num,DM_occ,size(DM_occ,1),mo_coef,size(mo_coef,1),-1.d0,rank)
print *, rank
TOUCH mo_coef mo_label
call orthonormalize_mos
call find_rotation(eigenvectors_Fock_matrix_mo,mo_tot_num_align,mo_coef,ao_num,R, mo_tot_num)
do i=1,mo_tot_num
iorder(i) = i
E_new(i) = 0.d0
do k=1,mo_tot_num
E_new(i) += R(k,i)*R(k,i)*diagonal_fock_matrix_mo(k)
enddo
enddo
call dsort(E_new(1),iorder(1),mo_tot_num)
eigenvectors_Fock_matrix_mo = mo_coef
do j=1,mo_tot_num
do i=1,ao_num
mo_coef(i,j) = eigenvectors_Fock_matrix_mo(i,iorder(j))
enddo
enddo
TOUCH mo_coef mo_label mo_integrals_map
endif
call clear_mo_map
SCF_energy_before = SCF_energy_after
SCF_energy_after = HF_energy
print*,SCF_energy_after, dabs(SCF_energy_before - SCF_energy_after)
i_it +=1
@ -31,6 +65,7 @@ program scf_iteration
stop 'Failed'
endif
mo_label = "Canonical"
deallocate (DM_occ)
TOUCH mo_label mo_coef
call save_mos

15
src/MOs/cholesky_mo.irp.f
View File

@ -1,11 +1,11 @@
subroutine cholesky_mo(n,m,P,C,LDC,tol_in,rank)
subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank)
implicit none
BEGIN_DOC
! Cholesky decomposition of MO Density matrix to
! generate MOs
END_DOC
integer, intent(in) :: n,m, LDC
double precision, intent(in) :: P(LDC,n)
integer, intent(in) :: n,m, LDC, LDP
double precision, intent(in) :: P(LDP,n)
double precision, intent(out) :: C(LDC,m)
double precision, intent(in) :: tol_in
integer, intent(out) :: rank
@ -49,7 +49,7 @@ BEGIN_PROVIDER [ double precision, mo_density_matrix, (mo_tot_num_align, mo_tot_
integer :: i,j,k
mo_density_matrix = 0.d0
do k=1,mo_tot_num
if (mo_occ(k) == 0) then
if (mo_occ(k) == 0.d0) then
cycle
endif
do j=1,ao_num
@ -67,14 +67,11 @@ BEGIN_PROVIDER [ double precision, mo_density_matrix_virtual, (mo_tot_num_align,
! Density matrix in MO basis (virtual MOs)
END_DOC
integer :: i,j,k
mo_density_matrix = 0.d0
mo_density_matrix_virtual = 0.d0
do k=1,mo_tot_num
if (mo_occ(k) == 0) then
cycle
endif
do j=1,ao_num
do i=1,ao_num
mo_density_matrix(i,j) = mo_density_matrix(i,j) + &
mo_density_matrix_virtual(i,j) = mo_density_matrix_virtual(i,j) + &
(2.d0-mo_occ(k)) * mo_coef(i,k) * mo_coef(j,k)
enddo
enddo