PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-05-06 23:24:58 +02:00
Code Releases Activity

added complex aotomo routines

Browse Source
This commit is contained in:
Loris Burth 2025-03-26 19:13:31 +01:00
parent a7d94e8193
commit 6d878d7a65
5 changed files with 141 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
src/AOtoMO/complex_AOtoMO.f90 Normal file
View File

@ -0,0 +1,24 @@
subroutine complex_AOtoMO(nBas, nOrb, C, M_AOs, M_MOs)
! Perform AO to MO transformation of a matrix M_AOs for given coefficients c
! M_MOs = C.T M_AOs C
implicit none
integer, intent(in) :: nBas, nOrb
complex*16, intent(in) :: C(nBas,nOrb)
double precision, intent(in) :: M_AOs(nBas,nBas)
complex*16, intent(out) :: M_MOs(nOrb,nOrb)
complex*16, allocatable :: AC(:,:)
complex*16, allocatable :: complex_C(:,:)
allocate(AC(nBas,nOrb))
AC = matmul(M_AOs, C)
M_MOs = matmul(transpose(C), AC)
deallocate(AC)
end subroutine

58
src/AOtoMO/complex_AOtoMO_ERI_RHF.f90 Normal file
View File

@ -0,0 +1,58 @@
subroutine complex_AOtoMO_ERI_RHF(nBas,nOrb,c,ERI_AO,ERI_MO)
! AO to MO transformation of two-electron integrals via the semi-direct O(N^5) algorithm
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
complex*16,intent(in) :: c(nBas,nOrb)
! Local variables
complex*16,allocatable :: a1(:,:,:,:)
complex*16,allocatable :: a2(:,:,:,:)
complex*16,allocatable :: complex_ERI_AO(:,:,:,:)
! Output variables
complex*16,intent(out) :: ERI_MO(nOrb,nOrb,nOrb,nOrb)
! Memory allocation
allocate(a2(nBas,nBas,nBas,nOrb))
allocate(a1(nBas,nBas,nOrb,nOrb))
allocate(complex_ERI_AO(nBas,nBas,nBas,nBas))
complex_ERI_AO = (1d0,0d0)*ERI_AO
! Four-index transform via semi-direct O(N^5) algorithm
call zgemm( 'T', 'N', nBas*nBas*nBas, nOrb, nBas, 1.d0 &
, complex_ERI_AO(1,1,1,1), nBas, c(1,1), nBas&
, 0.d0, a2(1,1,1,1), nBas*nBas*nBas)
deallocate(complex_ERI_AO)
call zgemm( 'T', 'N', nBas*nBas*nOrb, nOrb, nBas, 1.d0 &
, a2(1,1,1,1), nBas, c(1,1), nBas &
, 0.d0, a1(1,1,1,1), nBas*nBas*nOrb)
deallocate(a2)
allocate(a2(nBas,nOrb,nOrb,nOrb))
call zgemm( 'T', 'N', nBas*nOrb*nOrb, nOrb, nBas, 1.d0 &
, a1(1,1,1,1), nBas, c(1,1), nBas &
, 0.d0, a2(1,1,1,1), nBas*nOrb*nOrb)
deallocate(a1)
call zgemm( 'T', 'N', nOrb*nOrb*nOrb, nOrb, nBas, 1.d0 &
, a2(1,1,1,1), nBas, c(1,1), nBas &
, 0.d0, ERI_MO(1,1,1,1), nOrb*nOrb*nOrb)
deallocate(a2)
end subroutine

9
src/HF/cRHF.f90
View File

@ -39,6 +39,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
double precision :: Conv
double precision :: rcond
double precision,external :: trace_matrix
complex*16,external :: complex_trace_matrix
complex*16,allocatable :: J(:,:)
complex*16,allocatable :: K(:,:)
@ -132,18 +133,18 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
! CAP energy
EW = cmplx(trace_matrix(nBas,real(matmul(P,CAP))),trace_matrix(nBas,aimag(matmul(P,CAP))),kind=8)
EW = complex_trace_matrix(nBas,matmul(P,(0d0,1d0)*CAP))
write(*,*) "EW", EW
! Hartree energy
EJ = 0.5d0*cmplx(trace_matrix(nBas,real(matmul(P,J))),trace_matrix(nBas,aimag(matmul(P,J))),kind=8)
! Exchange energy
EK = 0.25d0*cmplx(trace_matrix(nBas,real(matmul(P,K))),trace_matrix(nBas,aimag(matmul(P,K))),kind=8)
! Total energy
ERHF = ET + EV + EW + EJ + EK

45
src/QuAcK/RQuAcK.f90
View File

@ -113,9 +113,12 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
double precision :: ERHF
complex*16 :: complex_ERHF
double precision,allocatable :: CAP_MO(:,:)
complex*16,allocatable :: complex_CAP_MO(:,:)
double precision,allocatable :: dipole_int_MO(:,:,:)
complex*16,allocatable :: complex_dipole_int_MO(:,:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
double precision,allocatable :: ERI_MO(:,:,:,:)
complex*16,allocatable :: complex_ERI_MO(:,:,:,:)
integer :: ixyz
integer :: nS
@ -133,16 +136,20 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
allocate(complex_eHF(nOrb))
allocate(complex_cHF(nBas,nOrb))
allocate(complex_FHF(nBas,nBas))
allocate(complex_dipole_int_MO(nOrb,nOrb,ncart))
allocate(complex_ERI_MO(nOrb,nOrb,nOrb,nOrb))
if (doCAP) allocate(complex_CAP_MO(nOrb,nOrb))
else
allocate(PHF(nBas,nBas))
allocate(eHF(nOrb))
allocate(cHF(nBas,nOrb))
allocate(FHF(nBas,nBas))
allocate(dipole_int_MO(nOrb,nOrb,ncart))
allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb))
if (doCAP) allocate(CAP_MO(nOrb,nOrb))
end if
allocate(ERI_AO(nBas,nBas,nBas,nBas))
allocate(dipole_int_MO(nOrb,nOrb,ncart))
allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb))
if (doCAP) allocate(CAP_MO(nOrb,nOrb))
call wall_time(start_int)
call read_2e_integrals(working_dir,nBas,ERI_AO)
call wall_time(end_int)
@ -205,16 +212,20 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
if (docRHF) then
! ! Transform from to complex MOs
!
! ! Read and transform dipole-related integrals
! do ixyz=1,ncart
! call AOtoMO(nBas,nOrb,complex_cHF,dipole_int_AO(1,1,ixyz),dipole_int_MO(1,1,ixyz))
! end do
! ! 4-index transform
! call AOtoMO_ERI_RHF(nBas,nOrb,complex_cHF,ERI_AO,ERI_MO)
! ! Transform CAP integrals
! if (doCAP) call AOtoMO(nBas,nOrb,complex_cHF,CAP_AO,CAP_MO)
! Transform from to complex MOs
! Read and transform dipole-related integrals
do ixyz=1,ncart
call complex_AOtoMO(nBas,nOrb,complex_cHF,dipole_int_AO(1,1,ixyz),complex_dipole_int_MO(1,1,ixyz))
end do
! 4-index transform
call complex_AOtoMO_ERI_RHF(nBas,nOrb,complex_cHF,ERI_AO,complex_ERI_MO)
! Transform CAP integrals
if (doCAP) then
call complex_AOtoMO(nBas,nOrb,complex_cHF,CAP_AO,complex_CAP_MO)
complex_CAP_MO = (0d0,1d0)*complex_CAP_MO
call complex_matout(nBas,nOrb,complex_CAP_MO)
end if
else
! Transform to real MOs
@ -423,10 +434,12 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3
if (allocated(cHF)) deallocate(cHF)
if (allocated(PHF)) deallocate(PHF)
if (allocated(FHF)) deallocate(FHF)
deallocate(dipole_int_MO)
deallocate(ERI_MO)
deallocate(ERI_AO)
if (allocated(dipole_int_MO)) deallocate(dipole_int_MO)
if (allocated(ERI_MO)) deallocate(ERI_MO)
if (allocated(complex_ERI_MO)) deallocate(complex_ERI_MO)
if (allocated(ERI_AO)) deallocate(ERI_AO)
if (allocated(CAP_MO)) deallocate(CAP_MO)
if (allocated(complex_CAP_MO)) deallocate(complex_CAP_MO)
if (allocated(complex_eHF)) deallocate(complex_eHF)
if (allocated(complex_cHF)) deallocate(complex_cHF)
if (allocated(complex_PHF)) deallocate(complex_PHF)

25
src/utils/utils.f90
View File

@ -214,6 +214,31 @@ function trace_matrix(n,A) result(Tr)
end function
function complex_trace_matrix(n,A) result(Tr)
! Calculate the trace of the square matrix A
implicit none
! Input variables
integer,intent(in) :: n
complex*16,intent(in) :: A(n,n)
! Local variables
integer :: i
! Output variables
complex*16 :: Tr
Tr = (0d0,0d0)
do i=1,n
Tr = Tr + A(i,i)
end do
end function
!------------------------------------------------------------------------
subroutine compute_error(nData,Mean,Var,Error)