From 6c930b70f589dc5335f9f31f9909e22cce911c20 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 11 Nov 2020 01:12:52 +0100 Subject: [PATCH] Added restoration of symmetry --- src/determinants/density_matrix.irp.f | 2 + src/hartree_fock/scf.irp.f | 19 +++--- src/mo_basis/EZFIO.cfg | 1 - src/mo_guess/h_core_guess_routine.irp.f | 1 + src/mo_one_e_ints/orthonormalize.irp.f | 13 ++-- src/scf_utils/EZFIO.cfg | 5 ++ src/scf_utils/huckel.irp.f | 1 + src/scf_utils/roothaan_hall_scf.irp.f | 4 +- src/tools/save_ortho_mos.irp.f | 2 + src/utils/linear_algebra.irp.f | 89 +++++++++++++++++++++++++ 10 files changed, 122 insertions(+), 15 deletions(-) diff --git a/src/determinants/density_matrix.irp.f b/src/determinants/density_matrix.irp.f index e69a1803..7c4a7fec 100644 --- a/src/determinants/density_matrix.irp.f +++ b/src/determinants/density_matrix.irp.f @@ -280,6 +280,8 @@ subroutine save_natural_mos ! the |MO| basis END_DOC call set_natural_mos + call nullify_small_elements(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10) + call orthonormalize_mos call save_mos end diff --git a/src/hartree_fock/scf.irp.f b/src/hartree_fock/scf.irp.f index 6ebb1b80..ace4e354 100644 --- a/src/hartree_fock/scf.irp.f +++ b/src/hartree_fock/scf.irp.f @@ -1,34 +1,34 @@ program scf BEGIN_DOC -! +! ! The :ref:`scf` program performs *Restricted* Hartree-Fock ! calculations (the spatial part of the |MOs| is common for alpha and beta ! spinorbitals). -! +! ! It performs the following actions: -! +! ! #. Compute/Read all the one- and two-electron integrals, and store them ! in memory ! #. Check in the |EZFIO| database if there is a set of |MOs|. ! If there is, it will read them as initial guess. Otherwise, it will ! create a guess. ! #. Perform the |SCF| iterations -! +! ! For the keywords related to the |SCF| procedure, see the ``scf_utils`` ! directory where you will find all options. -! +! ! At each iteration, the |MOs| are saved in the |EZFIO| database. Hence, ! if the calculation crashes for any unexpected reason, the calculation ! can be restarted by running again the |SCF| with the same |EZFIO| ! database. -! +! ! To start again a fresh |SCF| calculation, the |MOs| can be reset by ! running the :ref:`qp_reset` command. -! -! The `DIIS`_ algorithm is implemented, as well as the `level-shifting`_ +! +! The `DIIS`_ algorithm is implemented, as well as the `level-shifting`_ ! method. If the |SCF| does not converge, try again with a higher value of ! :option:`level_shift`. -! +! ! .. _DIIS: https://en.wikipedia.org/w/index.php?title=DIIS ! .. _level-shifting: https://doi.org/10.1002/qua.560070407 ! @@ -55,6 +55,7 @@ subroutine create_guess size(mo_one_e_integrals,1), & size(mo_one_e_integrals,2), & mo_label,1,.false.) + call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-10) SOFT_TOUCH mo_coef mo_label else if (mo_guess_type == "Huckel") then call huckel_guess diff --git a/src/mo_basis/EZFIO.cfg b/src/mo_basis/EZFIO.cfg index 874af46a..81ffba5c 100644 --- a/src/mo_basis/EZFIO.cfg +++ b/src/mo_basis/EZFIO.cfg @@ -36,4 +36,3 @@ size: (mo_basis.mo_num) type: character*(32) doc: MD5 checksum characterizing the |AO| basis set. interface: ezfio - diff --git a/src/mo_guess/h_core_guess_routine.irp.f b/src/mo_guess/h_core_guess_routine.irp.f index 8fc3f6f2..1073f052 100644 --- a/src/mo_guess/h_core_guess_routine.irp.f +++ b/src/mo_guess/h_core_guess_routine.irp.f @@ -8,6 +8,7 @@ subroutine hcore_guess call mo_as_eigvectors_of_mo_matrix(mo_one_e_integrals, & size(mo_one_e_integrals,1), & size(mo_one_e_integrals,2),label,1,.false.) + call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-12 ) call save_mos SOFT_TOUCH mo_coef mo_label end diff --git a/src/mo_one_e_ints/orthonormalize.irp.f b/src/mo_one_e_ints/orthonormalize.irp.f index 3a5d5488..c7cfbd90 100644 --- a/src/mo_one_e_ints/orthonormalize.irp.f +++ b/src/mo_one_e_ints/orthonormalize.irp.f @@ -1,11 +1,16 @@ subroutine orthonormalize_mos implicit none - integer :: m,p,s + integer :: m,p,s,i m = size(mo_coef,1) p = size(mo_overlap,1) - call ortho_lowdin(mo_overlap,p,mo_num,mo_coef,m,ao_num,lin_dep_cutoff) - mo_label = 'Orthonormalized' - SOFT_TOUCH mo_coef mo_label + do i=1,4 + call ortho_lowdin(mo_overlap,p,mo_num,mo_coef,m,ao_num,lin_dep_cutoff) + call nullify_small_elements(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10) + enddo + if (restore_symm) then + call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-12) + endif + SOFT_TOUCH mo_coef end diff --git a/src/scf_utils/EZFIO.cfg b/src/scf_utils/EZFIO.cfg index 4a56a35b..d69e278b 100644 --- a/src/scf_utils/EZFIO.cfg +++ b/src/scf_utils/EZFIO.cfg @@ -51,3 +51,8 @@ doc: If true, leave untouched all the orbitals defined as core and optimize all interface: ezfio,provider,ocaml default: False +[restore_symm] +type: logical +doc: If true, try to find symmetry in the MO coefficient matrices +interface: ezfio,provider,ocaml +default: True diff --git a/src/scf_utils/huckel.irp.f b/src/scf_utils/huckel.irp.f index 2d110e69..304d5144 100644 --- a/src/scf_utils/huckel.irp.f +++ b/src/scf_utils/huckel.irp.f @@ -25,6 +25,7 @@ subroutine huckel_guess TOUCH Fock_matrix_ao_alpha Fock_matrix_ao_beta mo_coef = eigenvectors_fock_matrix_mo + call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-12 ) SOFT_TOUCH mo_coef call save_mos deallocate(A) diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index dc44e262..e880fadf 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -165,6 +165,8 @@ END_DOC if(.not.frozen_orb_scf)then call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.) + call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-10) + call orthonormalize_mos call save_mos endif @@ -303,7 +305,7 @@ END_DOC Fock_matrix_AO_(i,j) = 0.d0 enddo do k=1,dim_DIIS - if (dabs(X_vector_DIIS(k)) < 1.d-10) cycle + if (dabs(X_vector_DIIS(k)) < 1.d-12) cycle do i=1,ao_num Fock_matrix_AO_(i,j) = Fock_matrix_AO_(i,j) + & X_vector_DIIS(k)*Fock_matrix_DIIS(i,j,dim_DIIS-k+1) diff --git a/src/tools/save_ortho_mos.irp.f b/src/tools/save_ortho_mos.irp.f index 6d998102..cb7a2d55 100644 --- a/src/tools/save_ortho_mos.irp.f +++ b/src/tools/save_ortho_mos.irp.f @@ -10,5 +10,7 @@ program save_ortho_mos ! Thanks to the Lowdin orthonormalization, the new MOs are the most similar to the guess MOs. END_DOC call orthonormalize_mos + mo_label = 'Orthonormalized' + SOFT_TOUCH mo_label call save_mos end diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index deb4f4ca..76e3a9f1 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1367,3 +1367,92 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) deallocate(A,eigenvalues) end +subroutine nullify_small_elements(m,n,A,LDA,thresh) + implicit none + integer, intent(in) :: m,n,LDA + double precision, intent(inout) :: A(LDA,n) + double precision, intent(in) :: thresh + integer :: i,j + double precision :: amax + + ! Find max value + amax = 0.d0 + do j=1,n + do i=1,m + amax = max(dabs(A(i,j)), amax) + enddo + enddo + amax = 1.d0/amax + + ! Remove tiny elements + do j=1,n + do i=1,m + if ( dabs(A(i,j) * amax) < thresh ) then + A(i,j) = 0.d0 + endif + enddo + enddo + +end + +subroutine restore_symmetry(m,n,A,LDA,thresh) + implicit none + integer, intent(in) :: m,n,LDA + double precision, intent(inout) :: A(LDA,n) + double precision, intent(in) :: thresh + integer :: i,j,k,l + logical, allocatable :: done(:,:) + double precision :: f, g, count, thresh2 + thresh2 = dsqrt(thresh) + call nullify_small_elements(m,n,A,LDA,thresh) + + allocate(done(m,n)) + + do j=1,n + do i=1,m + done(i,j) = A(i,j) == 0.d0 + enddo + enddo + + do j=1,n + do i=1,m + if ( done(i,j) ) cycle + done(i,j) = .True. + count = 1.d0 + f = 1.d0/A(i,j) + do l=1,n + do k=1,m + if ( done(k,l) ) cycle + g = f * A(k,l) + if ( dabs(dabs(g) - 1.d0) < thresh2 ) then + count = count + 1.d0 + if (g>0.d0) then + A(i,j) = A(i,j) + A(k,l) + else + A(i,j) = A(i,j) - A(k,l) + end if + endif + enddo + enddo + if (count > 1.d0) then + A(i,j) = A(i,j) / count + do l=1,n + do k=1,m + if ( done(k,l) ) cycle + g = f * A(k,l) + if ( dabs(dabs(g) - 1.d0) < thresh2 ) then + done(k,l) = .True. + if (g>0.d0) then + A(k,l) = A(i,j) + else + A(k,l) = -A(i,j) + end if + endif + enddo + enddo + endif + + enddo + enddo + +end