PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-01-11 13:38:18 +01:00
Code Releases Activity

complex diis

Browse Source
This commit is contained in:
Kevin Gasperich 2020-01-29 15:41:23 -06:00
parent e64faf2845
commit b0d27f8503
1 changed files with 42 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
src/scf_utils/diis_complex.irp.f
View File

@ -47,7 +47,7 @@ BEGIN_PROVIDER [complex*16, FPS_SPF_Matrix_AO_complex, (AO_num, AO_num)]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [complex*16, FPS_SPF_Matrix_MO, (mo_num, mo_num)] BEGIN_PROVIDER [complex*16, FPS_SPF_Matrix_MO_complex, (mo_num, mo_num)]
implicit none implicit none
begin_doc begin_doc
! Commutator FPS - SPF in MO basis ! Commutator FPS - SPF in MO basis
@ -66,14 +66,13 @@ END_PROVIDER
implicit none implicit none
double precision, allocatable :: scratch(:,:),work(:),Xt(:,:) double precision, allocatable :: rwork(:)
integer :: lwork,info integer :: lwork,info,lrwork
complex*16, allocatable :: scratch(:,:),Xt(:,:),work(:)
integer :: i,j integer :: i,j
lwork = 3*AO_num - 1
allocate( & allocate( &
scratch(AO_num,AO_num), & scratch(AO_num,AO_num), &
work(lwork), &
Xt(AO_num,AO_num) & Xt(AO_num,AO_num) &
) )
@ -81,46 +80,63 @@ END_PROVIDER
do i=1,AO_num do i=1,AO_num
do j=1,AO_num do j=1,AO_num
Xt(i,j) = S_half_inv(j,i) Xt(i,j) = dconjg(S_half_inv_complex(j,i))
enddo enddo
enddo enddo
! Calculate Fock matrix in orthogonal basis: F' = Xt.F.X ! Calculate Fock matrix in orthogonal basis: F' = Xt.F.X
call dgemm('N','N',AO_num,AO_num,AO_num, & call zgemm('N','N',AO_num,AO_num,AO_num, &
1.d0, & (1.d0,0.d0), &
Fock_matrix_AO,size(Fock_matrix_AO,1), & Fock_matrix_AO_complex,size(Fock_matrix_AO_complex,1), &
S_half_inv,size(S_half_inv,1), & S_half_inv_complex,size(s_half_inv_complex,1), &
0.d0, & (0.d0,0.d0), &
eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1)) eigenvectors_Fock_matrix_AO_complex, &
size(eigenvectors_Fock_matrix_AO_complex,1))
call dgemm('N','N',AO_num,AO_num,AO_num, & call zgemm('N','N',AO_num,AO_num,AO_num, &
1.d0, & (1.d0,0.d0), &
Xt,size(Xt,1), & Xt,size(Xt,1), &
eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1), & eigenvectors_Fock_matrix_AO_complex, &
0.d0, & size(eigenvectors_Fock_matrix_AO_complex,1), &
(0.d0,0.d0), &
scratch,size(scratch,1)) scratch,size(scratch,1))
! Diagonalize F' to obtain eigenvectors in orthogonal basis C' and eigenvalues ! Diagonalize F' to obtain eigenvectors in orthogonal basis C' and eigenvalues
lrwork = 3*ao_num - 2
allocate(rwork(lrwork), work(1))
lwork = -1
call dsyev('V','U',AO_num, & call zheev('V','U',ao_num, &
scratch,size(scratch,1), & scratch,size(scratch,1), &
eigenvalues_Fock_matrix_AO, & eigenvalues_Fock_matrix_AO_complex, &
work,lwork,info) work,lwork,rwork,info)
lwork = int(work(1))
deallocate(work)
allocate(work(lwork))
call zheev('V','U',ao_num, &
scratch,size(scratch,1), &
eigenvalues_Fock_matrix_AO_complex, &
work,lwork,rwork,info)
if(info /= 0) then if(info /= 0) then
print *, irp_here//' failed : ', info print *, irp_here//' failed : ', info
stop 1 stop 1
endif endif
deallocate(work,rwork)
! Back-transform eigenvectors: C =X.C' ! Back-transform eigenvectors: C =X.C'
call dgemm('N','N',AO_num,AO_num,AO_num, & call zgemm('N','N',AO_num,AO_num,AO_num, &
1.d0, & (1.d0,0.d0), &
S_half_inv,size(S_half_inv,1), & S_half_inv_complex,size(S_half_inv_complex,1), &
scratch,size(scratch,1), & scratch,size(scratch,1), &
0.d0, & (0.d0,0.d0), &
eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1)) eigenvectors_Fock_matrix_AO_complex, &
size(eigenvectors_Fock_matrix_AO_complex,1))
deallocate(scratch)
END_PROVIDER END_PROVIDER