From 42d8b4c4043d9b8b013303f403472ffc215bdf99 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 12 Jun 2014 22:08:53 +0200 Subject: [PATCH] Improved Hartree-Fock --- src/AOs/aos.irp.f | 40 ---------------- src/Hartree_Fock/HF_density_matrix_ao.irp.f | 32 ++++++------- src/Hartree_Fock/mo_SCF_iterations.irp.f | 53 +++++++++++++++++---- src/MOs/cholesky_mo.irp.f | 15 +++--- 4 files changed, 66 insertions(+), 74 deletions(-) diff --git a/src/AOs/aos.irp.f b/src/AOs/aos.irp.f index 4b3a74c0..396d9aab 100644 --- a/src/AOs/aos.irp.f +++ b/src/AOs/aos.irp.f @@ -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 diff --git a/src/Hartree_Fock/HF_density_matrix_ao.irp.f b/src/Hartree_Fock/HF_density_matrix_ao.irp.f index b0d58344..3183025b 100644 --- a/src/Hartree_Fock/HF_density_matrix_ao.irp.f +++ b/src/Hartree_Fock/HF_density_matrix_ao.irp.f @@ -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 diff --git a/src/Hartree_Fock/mo_SCF_iterations.irp.f b/src/Hartree_Fock/mo_SCF_iterations.irp.f index 01cab413..cf38dbd1 100644 --- a/src/Hartree_Fock/mo_SCF_iterations.irp.f +++ b/src/Hartree_Fock/mo_SCF_iterations.irp.f @@ -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 - mo_coef = eigenvectors_Fock_matrix_mo - TOUCH mo_coef mo_label + 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 diff --git a/src/MOs/cholesky_mo.irp.f b/src/MOs/cholesky_mo.irp.f index d3e04dfb..9afea3e6 100644 --- a/src/MOs/cholesky_mo.irp.f +++ b/src/MOs/cholesky_mo.irp.f @@ -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