LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-08-28 23:32:04 +02:00
Code Releases Activity

optim in AdAt

Browse Source
This commit is contained in:
Abdallah Ammar 2024-08-21 13:51:12 +02:00
parent a1786b2ade
commit 992e3dff4b
4 changed files with 32 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
src/RPA/sort_ppRPA.f90
View File

@ -39,9 +39,10 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
double precision,intent(out) :: X2(nVV,nOO)
double precision,intent(out) :: Y2(nOO,nOO)
! Memory allocation
allocate(M(nOO+nVV,nOO+nVV),Z1(nOO+nVV,nVV),Z2(nOO+nVV,nOO),order1(nVV),order2(nOO))
allocate(M(nOO+nVV,nOO+nVV),Z1(nOO+nVV,nVV),Z2(nOO+nVV,nOO),order1(nVV),order2(nOO))
! Initializatiom
@ -86,7 +87,7 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
end if
end do
end do
if(minval(Om1) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!')
if(maxval(Om2) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!')
@ -111,7 +112,8 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
call quick_sort(Om2,order2,nOO)
call set_order(Z2,order2,nOO+nVV,nOO)
end if
end if
! Orthogonalize eigenvectors
@ -202,7 +204,6 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
if(nVV > 0) call dgemm ('N', 'N', nOO+nVV, nVV, nOO+nVV, 1d0, M, nOO+nVV, Z1, nOO+nVV, 0d0, tmp1, nOO+nVV)
if(nVV > 0) call dgemm ('T', 'N', nVV , nVV, nOO+nVV, 1d0, Z1, nOO+nVV, tmp1, nOO+nVV, 0d0, S1, nVV)
if(nOO > 0) call dgemm ('N', 'N', nOO+nVV, nOO, nOO+nVV, 1d0, M, nOO+nVV, -1d0*Z2, nOO+nVV, 0d0, tmp2, nOO+nVV)
if(nOO > 0) call dgemm ('T', 'N', nOO , nOO, nOO+nVV, 1d0, Z2, nOO+nVV, tmp2, nOO+nVV, 0d0, S2, nOO)

2
src/utils/orthogonalization_matrix.f90
View File

@ -54,7 +54,7 @@ subroutine orthogonalization_matrix(nBas,S,X)
end do
call ADAt(nBas,Uvec,Uval,X)
call ADAt(nBas, Uvec(1,1), Uval(1), X(1,1))
elseif(ortho_type == 2) then

30
src/utils/utils.f90
View File

@ -375,15 +375,29 @@ subroutine ADAt(N,A,D,B)
double precision,intent(out) :: B(N,N)
B = 0d0
double precision, allocatable :: tmp(:,:)
do i=1,N
do j=1,N
do k=1,N
B(i,k) = B(i,k) + A(i,j)*D(j)*A(k,j)
end do
end do
end do
allocate(tmp(N,N))
!$OMP PARALLEL DEFAULT(NONE) PRIVATE(i, j) SHARED(N, A, D, tmp)
!$OMP DO
do i = 1, N
do j = 1, N
tmp(i,j) = D(i) * A(j,i)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm("N", "N", N, N, N, 1.d0, A(1,1), N, tmp(1,1), N, 0.d0, B(1,1), N)
deallocate(tmp)
! B = 0d0
! do i=1,N
! do j=1,N
! do k=1,N
! B(i,k) = B(i,k) + A(i,j)*D(j)*A(k,j)
! end do
! end do
! end do
end subroutine
!------------------------------------------------------------------------

4
src/utils/wrap_lapack.f90
View File

@ -31,6 +31,8 @@ subroutine diagonalize_general_matrix(N,A,WR,VR)
call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info)
deallocate(work, WI, VL)
if(info /= 0) then
print*,'Problem in diagonalize_general_matrix (dgeev)!!'
end if
@ -67,6 +69,8 @@ subroutine diagonalize_matrix(N,A,e)
allocate(work(lwork))
call dsyev('V','U',N,A,N,e,work,lwork,info)
deallocate(work)
if(info /= 0) then
print*,'Problem in diagonalize_matrix (dsyev)!!'