LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-30 15:15:38 +01:00
Code Releases Activity

Symmetry

Browse Source
This commit is contained in:
Anthony Scemama 2020-11-11 15:51:19 +01:00
parent 1b6abd2601
commit 7c100a1a2d
8 changed files with 195 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/ao_one_e_ints/ao_ortho_canonical.irp.f
View File

@ -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)]

2
src/cisd/30.cisd.bats
View File

@ -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

3
src/hartree_fock/scf.irp.f
View File

@ -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

1
src/mo_basis/mos.irp.f
View File

@ -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

1
src/mo_one_e_ints/mo_one_e_ints.irp.f
View File

@ -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

4
src/mo_one_e_ints/orthonormalize.irp.f
View File

@ -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

2
src/scf_utils/huckel.irp.f
View File

@ -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

191
src/utils/linear_algebra.irp.f
View File

@ -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