From 7c100a1a2d71114687c544ccd52b8279ab303823 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 11 Nov 2020 15:51:19 +0100 Subject: [PATCH] Symmetry --- src/ao_one_e_ints/ao_ortho_canonical.irp.f | 1 + src/cisd/30.cisd.bats | 2 +- src/hartree_fock/scf.irp.f | 3 +- src/mo_basis/mos.irp.f | 1 + src/mo_one_e_ints/mo_one_e_ints.irp.f | 1 + src/mo_one_e_ints/orthonormalize.irp.f | 4 +- src/scf_utils/huckel.irp.f | 2 +- src/utils/linear_algebra.irp.f | 191 ++++++++++++++++++++- 8 files changed, 195 insertions(+), 10 deletions(-) diff --git a/src/ao_one_e_ints/ao_ortho_canonical.irp.f b/src/ao_one_e_ints/ao_ortho_canonical.irp.f index 45275a06..668b920d 100644 --- a/src/ao_one_e_ints/ao_ortho_canonical.irp.f +++ b/src/ao_one_e_ints/ao_ortho_canonical.irp.f @@ -137,6 +137,7 @@ END_PROVIDER deallocate(S) endif + END_PROVIDER BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)] diff --git a/src/cisd/30.cisd.bats b/src/cisd/30.cisd.bats index 44365174..ba39d446 100644 --- a/src/cisd/30.cisd.bats +++ b/src/cisd/30.cisd.bats @@ -110,7 +110,7 @@ function run() { [[ -n $TRAVIS ]] && skip qp set_file cu_nh3_4_2plus.ezfio qp set_mo_class --core="[1-24]" --act="[25-45]" --del="[46-87]" - run -1862.98611018932 -1862.68751252590 + run -1862.98676507415 -1862.68809564033 } @test "ClF" { # 30.3225s diff --git a/src/hartree_fock/scf.irp.f b/src/hartree_fock/scf.irp.f index 0368c347..2ed02bee 100644 --- a/src/hartree_fock/scf.irp.f +++ b/src/hartree_fock/scf.irp.f @@ -50,12 +50,13 @@ subroutine create_guess mo_label = 'Guess' if (mo_guess_type == "HCore") then mo_coef = ao_ortho_lowdin_coef + call restore_symmetry(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10) TOUCH mo_coef call mo_as_eigvectors_of_mo_matrix(mo_one_e_integrals, & 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) + call restore_symmetry(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10) SOFT_TOUCH mo_coef else if (mo_guess_type == "Huckel") then call huckel_guess diff --git a/src/mo_basis/mos.irp.f b/src/mo_basis/mos.irp.f index 73d33901..57ebb533 100644 --- a/src/mo_basis/mos.irp.f +++ b/src/mo_basis/mos.irp.f @@ -277,6 +277,7 @@ subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo) T, ao_num, & 0.d0, A_mo, size(A_mo,1)) + call restore_symmetry(mo_num,mo_num,A_mo,size(A_mo,1),1.d-12) deallocate(T) end diff --git a/src/mo_one_e_ints/mo_one_e_ints.irp.f b/src/mo_one_e_ints/mo_one_e_ints.irp.f index a6a614ab..926ac7bd 100644 --- a/src/mo_one_e_ints/mo_one_e_ints.irp.f +++ b/src/mo_one_e_ints/mo_one_e_ints.irp.f @@ -18,5 +18,6 @@ BEGIN_PROVIDER [ double precision, mo_one_e_integrals,(mo_num,mo_num)] call ezfio_set_mo_one_e_ints_mo_one_e_integrals(mo_one_e_integrals) print *, 'MO one-e integrals written to disk' ENDIF + call nullify_small_elements(mo_num,mo_num,mo_one_e_integrals,size(mo_one_e_integrals,1),1.d-10) END_PROVIDER diff --git a/src/mo_one_e_ints/orthonormalize.irp.f b/src/mo_one_e_ints/orthonormalize.irp.f index ca382535..94ece41f 100644 --- a/src/mo_one_e_ints/orthonormalize.irp.f +++ b/src/mo_one_e_ints/orthonormalize.irp.f @@ -7,9 +7,7 @@ subroutine orthonormalize_mos 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-10) - endif + call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-10) SOFT_TOUCH mo_coef end diff --git a/src/scf_utils/huckel.irp.f b/src/scf_utils/huckel.irp.f index a6087da3..676cc345 100644 --- a/src/scf_utils/huckel.irp.f +++ b/src/scf_utils/huckel.irp.f @@ -25,7 +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 ) + call restore_symmetry(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10) call orthonormalize_mos SOFT_TOUCH mo_coef call save_mos diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 76e3a9f1..26d0b78c 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -5,7 +5,7 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) ! ! LDx : leftmost dimension of x ! - ! Dimsneion of A is m x n + ! Dimension of A is m x n ! END_DOC @@ -24,22 +24,190 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) ! Find optimal size for temp arrays allocate(work(1)) lwork = -1 - call dgesvd('S','S', m, n, A_tmp, LDA, & + call dgesvd('A','A', m, n, A_tmp, LDA, & D, U, LDU, Vt, LDVt, work, lwork, info) ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648 lwork = max(int(work(1)), 5*MIN(M,N)) deallocate(work) allocate(work(lwork)) - call dgesvd('S','S', m, n, A_tmp, LDA, & + call dgesvd('A','A', m, n, A_tmp, LDA, & D, U, LDU, Vt, LDVt, work, lwork, info) - deallocate(work,A_tmp) + deallocate(A_tmp,work) if (info /= 0) then print *, info, ': SVD failed' stop endif + do j=1,m + do i=1,m + if (dabs(U(i,j)) < 1.d-14) U(i,j) = 0.d0 + enddo + enddo + + do j=1,n + do i=1,n + if (dabs(Vt(i,j)) < 1.d-14) Vt(i,j) = 0.d0 + enddo + enddo + +end + +subroutine svd_symm(A,LDA,U,LDU,D,Vt,LDVt,m,n) + implicit none + BEGIN_DOC + ! Compute A = U.D.Vt + ! + ! LDx : leftmost dimension of x + ! + ! Dimension of A is m x n + ! + END_DOC + + integer, intent(in) :: LDA, LDU, LDVt, m, n + double precision, intent(in) :: A(LDA,n) + double precision, intent(out) :: U(LDU,min(m,n)) + double precision,intent(out) :: Vt(LDVt,n) + double precision,intent(out) :: D(min(m,n)) + double precision,allocatable :: work(:) + integer :: info, lwork, i, j, k + + double precision,allocatable :: A_tmp(:,:) + allocate (A_tmp(LDA,n)) + A_tmp(:,:) = A(:,:) + + ! Find optimal size for temp arrays + allocate(work(1)) + lwork = -1 + call dgesvd('A','A', m, n, A_tmp, LDA, & + D, U, LDU, Vt, LDVt, work, lwork, info) + ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648 + lwork = max(int(work(1)), 5*MIN(M,N)) + deallocate(work) + + allocate(work(lwork)) + call dgesvd('A','A', m, n, A_tmp, LDA, & + D, U, LDU, Vt, LDVt, work, lwork, info) + deallocate(A_tmp,work) + + if (info /= 0) then + print *, info, ': SVD failed' + stop + endif + + ! Iterative refinement + ! -------------------- + ! https://doi.org/10.1016/j.cam.2019.112512 + + integer :: iter + double precision,allocatable :: R(:,:), S(:,:), T(:,:) + double precision,allocatable :: sigma(:), F(:,:), G(:,:) + double precision :: alpha, beta, x, thresh + + allocate (R(m,m), S(n,n), T(m,n), sigma(m), F(m,m), G(n,n), A_tmp(m,n)) + sigma = 0.d0 + R = 0.d0 + S = 0.d0 + T = 0.d0 + F = 0.d0 + G = 0.d0 + + thresh = 1.d-8 + call restore_symmetry(m,m,U,size(U,1),thresh) + call restore_symmetry(n,n,Vt,size(Vt,1),thresh) + + do iter=1,4 + do k=1,n + A_tmp(1:m,k) = D(k) * U(1:m,k) + enddo + + call dgemm('N','N',m,n,n,1.d0,A_tmp,size(A_tmp,1), & + Vt,size(Vt,1),0.d0,R,size(R,1)) + + print *, maxval(dabs(R(1:m,1:n) - A(1:m,1:n))) + + + call dgemm('T','N',m,m,m,-1.d0,U,size(U,1), & + U,size(U,1),0.d0,R,size(R,1)) + do i=1,m + R(i,i) = R(i,i) + 1.d0 + enddo + + call dgemm('N','T',n,n,n,-1.d0,Vt,size(Vt,1), & + Vt,size(Vt,1),0.d0,S,size(S,1)) + do i=1,n + S(i,i) = S(i,i) + 1.d0 + enddo + + + call dgemm('T','N',m,n,m,1.d0,U,size(U,1), & + A,size(A,1),0.d0,A_tmp,size(A_tmp,1)) + + call dgemm('N','T',m,n,n,1.d0,A_tmp,size(A_tmp,1), & + Vt,size(Vt,1),0.d0,T,size(T,1)) + + + do i=1,n + sigma(i) = T(i,i)/(1.d0 - (R(i,i)+S(i,i))*0.5d0) + F(i,i) = 0.5d0*R(i,i) + G(i,i) = 0.5d0*S(i,i) + enddo + + do j=1,n + do i=1,n + if (i == j) cycle + alpha = T(i,j) + sigma(j) * R(i,j) + beta = T(i,j) + sigma(j) * S(i,j) + x = 1.d0 / (sigma(j)*sigma(j) - sigma(i)*sigma(i)) + F(i,j) = (alpha * sigma(j) + beta * sigma(i)) * x + G(i,j) = (alpha * sigma(i) + beta * sigma(j)) * x + enddo + enddo + + do i=1,n + x = 1.d0/sigma(i) + do j=n+1,m + F(i,j) = -T(j,i) * x + enddo + enddo + + do i=n+1,m + do j=1,n + F(i,j) = R(i,j) - F(j,i) + enddo + enddo + + do i=n+1,m + do j=n+1,m + F(i,j) = R(i,j)*0.5d0 + enddo + enddo + + D(1:min(n,m)) = sigma(1:min(n,m)) + call dgemm('N','N',m,m,m,1.d0,U,size(U,1),F,size(F,1), & + 0.d0, A_tmp, size(A_tmp,1)) + do j=1,m + do i=1,m + U(i,j) = U(i,j) + A_tmp(i,j) + enddo + enddo + + call dgemm('T','N',n,n,n,1.d0,G,size(G,1),Vt,size(Vt,1), & + 0.d0, A_tmp, size(A_tmp,1)) + do j=1,n + do i=1,n + Vt(i,j) = Vt(i,j) + A_tmp(i,j) + enddo + enddo + + thresh = 0.01d0 * thresh + call restore_symmetry(m,m,U,size(U,1),thresh) + call restore_symmetry(n,n,Vt,size(Vt,1),thresh) + + enddo + + deallocate(A_tmp,R,S,F,G,sigma) end subroutine eigSVD(A,LDA,U,LDU,D,Vt,LDVt,m,n) @@ -1397,6 +1565,11 @@ end subroutine restore_symmetry(m,n,A,LDA,thresh) implicit none + BEGIN_DOC +! Tries to find the matrix elements that are the same, and sets them +! to the average value. +! If restore_symm is False, only nullify small elements + END_DOC integer, intent(in) :: m,n,LDA double precision, intent(inout) :: A(LDA,n) double precision, intent(in) :: thresh @@ -1406,6 +1579,16 @@ subroutine restore_symmetry(m,n,A,LDA,thresh) thresh2 = dsqrt(thresh) call nullify_small_elements(m,n,A,LDA,thresh) + if (.not.restore_symm) then + return + endif + + ! TODO: Costs O(n^4), but can be improved to (2 n^2 * log(n)): + ! - copy all values in a 1D array + ! - sort 1D array + ! - average nearby elements + ! - for all elements, find matching value in the sorted 1D array + allocate(done(m,n)) do j=1,n