Better Hartree-Fock

src/Hartree_Fock/README.rst

@@ -83,6 +83,27 @@ Documentation
 `hf_density_matrix_ao_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L/BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta,  (ao_num_align,ao_num) ]/;">`_
   Beta density matrix in the AO basis
 
+`fock_mo_to_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO)/;">`_
+  Undocumented
+
+`insert_new_scf_density_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/subroutine insert_new_SCF_density_matrix/;">`_
+  Undocumented
+
+`it_scf <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/BEGIN_PROVIDER [ integer, it_scf ]/;">`_
+  Number of the current SCF iteration
+
+`scf_density_matrices <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,n_it_scf_max) ]/;">`_
+  Density matrices at every SCF iteration
+
+`scf_energies <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/&BEGIN_PROVIDER [ double precision, SCF_energies, (n_it_scf_max) ]/;">`_
+  Density matrices at every SCF iteration
+
+`scf_interpolation_step <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/subroutine SCF_interpolation_step/;">`_
+  Undocumented
+
+`scf_iterations <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/subroutine scf_iterations/;">`_
+  Undocumented
+
 `diagonal_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L/BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (mo_tot_num) ]/;">`_
   Diagonal Fock matrix in the MO basis
 

src/Hartree_Fock/SCF.irp.f

@@ -6,12 +6,14 @@ BEGIN_PROVIDER [ integer, it_scf ]
   it_scf = 0
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,0:n_it_scf_max) ]
+ BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,n_it_scf_max) ]
+&BEGIN_PROVIDER [ double precision, SCF_energies, (n_it_scf_max) ]
   implicit none
   BEGIN_DOC
   ! Density matrices at every SCF iteration
   END_DOC
   SCF_density_matrices = 0.d0
+  SCF_energies = 0.d0
 END_PROVIDER
 
 subroutine insert_new_SCF_density_matrix
@@ -20,11 +22,10 @@ subroutine insert_new_SCF_density_matrix
   do j=1,ao_num
     do i=1,ao_num
       SCF_density_matrices(i,j,1,it_scf) = HF_density_matrix_ao_alpha(i,j)
-    SCF_density_matrices(i,j,1,0) += HF_density_matrix_ao_alpha(i,j)
       SCF_density_matrices(i,j,2,it_scf) = HF_density_matrix_ao_beta(i,j)
-    SCF_density_matrices(i,j,2,0) += HF_density_matrix_ao_beta(i,j)
     enddo
   enddo
+  SCF_energies(it_scf) = HF_energy
 end
 
 subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO)
@@ -60,27 +61,31 @@ subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO)
   deallocate(T,M)
 end
 
-subroutine DIIS_step
+subroutine SCF_interpolation_step
   implicit none
   integer                        :: i,j
   double precision               :: c
- c = 1.d0/dble(it_scf)
+  
+  if (it_scf == 1) then
+    return
+  endif
+  call random_number(c)
   do j=1,ao_num
     do i=1,ao_num
-   HF_density_matrix_ao_alpha(i,j) = SCF_density_matrices(i,j,1,0) * c
-   HF_density_matrix_ao_beta (i,j) = SCF_density_matrices(i,j,2,0) * c
+      HF_density_matrix_ao_alpha(i,j) = c*SCF_density_matrices(i,j,1,it_scf)+SCF_density_matrices(i,j,1,it_scf-1) * (1.d0 - c)
+      HF_density_matrix_ao_beta (i,j) = c*SCF_density_matrices(i,j,2,it_scf)+SCF_density_matrices(i,j,2,it_scf-1) * (1.d0 - c)
     enddo
   enddo
   TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta
   
-! call Fock_mo_to_ao(Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1), &
-!      Fock_matrix_alpha_ao, size(Fock_matrix_alpha_ao,1) )
-! call Fock_mo_to_ao(Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1), &
-!      Fock_matrix_beta_ao, size(Fock_matrix_beta_ao,1) )
-! SOFT_TOUCH Fock_matrix_alpha_ao Fock_matrix_beta_ao Fock_matrix_mo_alpha Fock_matrix_mo_beta
+  ! call Fock_mo_to_ao(Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1),&
+      !      Fock_matrix_alpha_ao, size(Fock_matrix_alpha_ao,1) )
+  ! call Fock_mo_to_ao(Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1),&
+      !      Fock_matrix_beta_ao, size(Fock_matrix_beta_ao,1) )
+  ! SOFT_TOUCH Fock_matrix_alpha_ao Fock_matrix_beta_ao Fock_matrix_mo_alpha Fock_matrix_mo_beta
 end
 
-subroutine scf_iteration
+subroutine scf_iterations
   implicit none
   integer                        :: i,j
   do i=1,n_it_scf_max
@@ -89,7 +94,14 @@ subroutine scf_iteration
     mo_coef = eigenvectors_Fock_matrix_mo
     TOUCH mo_coef
     call insert_new_SCF_density_matrix
-   call DIIS_step
     print *,  HF_energy
+    if (SCF_energies(it_scf)>SCF_energies(it_scf-1)) then
+      call SCF_interpolation_step
+    endif
+    if (it_scf>1 ) then
+      if (dabs(SCF_energies(it_scf)-SCF_energies(it_scf-1)) < thresh_SCF) then
+        exit
+      endif
+    endif
   enddo
 end

src/Hartree_Fock/mo_SCF_iterations.irp.f

@@ -6,32 +6,10 @@ program xcf_iteration
   integer                        :: i_it
   
   E0 = HF_energy 
-  i_it = 0
-  n_it_scf_max = 10
-  SCF_energy_before = huge(1.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
-    call clear_mo_map
-    SCF_energy_after = HF_energy
-    print*,SCF_energy_after, dabs(SCF_energy_before - SCF_energy_after) 
-    i_it +=1
-    if(i_it > n_it_scf_max)exit
-  enddo
-  
-  if (i_it >= n_it_scf_max) then
-    stop 'Failed'
-  endif
-  if (SCF_energy_after - E0 > thresh_SCF) then
-    stop 'Failed'
-  endif
+  call scf_iterations
   mo_label = "Canonical"
   TOUCH mo_label mo_coef
-!  call save_mos
+  call save_mos
   
 end