PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-23 12:55:37 +01:00
Code Releases Activity

more work on complex SCF

Browse Source
This commit is contained in:
Kevin Gasperich 2020-01-28 18:06:00 -06:00
parent a63ee551ef
commit 8bfcfe8f21
3 changed files with 462 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
src/bitmask/track_orb.irp.f
View File

@ -7,18 +7,32 @@ BEGIN_PROVIDER [ double precision, mo_coef_begin_iteration, (ao_num,mo_num) ]
END_DOC
END_PROVIDER
BEGIN_PROVIDER [ complex*16, mo_coef_begin_iteration_complex, (ao_num,mo_num) ]
implicit none
BEGIN_DOC
! Void provider to store the coefficients of the |MO| basis at the beginning of the SCF iteration
!
! Useful to track some orbitals
END_DOC
END_PROVIDER
subroutine initialize_mo_coef_begin_iteration
implicit none
BEGIN_DOC
!
! Initialize :c:data:`mo_coef_begin_iteration` to the current :c:data:`mo_coef`
END_DOC
mo_coef_begin_iteration = mo_coef
if (is_periodic) then
mo_coef_begin_iteration_complex = mo_coef_complex
else
mo_coef_begin_iteration = mo_coef
endif
end
subroutine reorder_core_orb
implicit none
BEGIN_DOC
! TODO: modify for complex
! routines that takes the current :c:data:`mo_coef` and reorder the core orbitals (see :c:data:`list_core` and :c:data:`n_core_orb`) according to the overlap with :c:data:`mo_coef_begin_iteration`
END_DOC
integer :: i,j,iorb

126
src/scf_utils/diis_complex.irp.f Normal file
View File

@ -0,0 +1,126 @@
BEGIN_PROVIDER [complex*16, FPS_SPF_Matrix_AO_complex, (AO_num, AO_num)]
implicit none
BEGIN_DOC
! Commutator FPS - SPF
END_DOC
complex*16, allocatable :: scratch(:,:)
allocate( &
scratch(AO_num, AO_num) &
)
! Compute FP
call zgemm('N','N',AO_num,AO_num,AO_num, &
(1.d0,0.d0), &
Fock_Matrix_AO_complex,Size(Fock_Matrix_AO_complex,1), &
SCF_Density_Matrix_AO_complex,Size(SCF_Density_Matrix_AO_complex,1), &
(0.d0,0.d0), &
scratch,Size(scratch,1))
! Compute FPS
call zgemm('N','N',AO_num,AO_num,AO_num, &
(1.d0,0.d0), &
scratch,Size(scratch,1), &
AO_Overlap_complex,Size(AO_Overlap_complex,1), &
(0.d0,0.d0), &
FPS_SPF_Matrix_AO_complex,Size(FPS_SPF_Matrix_AO_complex,1))
! Compute SP
call zgemm('N','N',AO_num,AO_num,AO_num, &
(1.d0,0.d0), &
AO_Overlap_complex,Size(AO_Overlap_complex,1), &
SCF_Density_Matrix_AO_complex,Size(SCF_Density_Matrix_AO_complex,1), &
(0.d0,0.d0), &
scratch,Size(scratch,1))
! Compute FPS - SPF
call zgemm('N','N',AO_num,AO_num,AO_num, &
(-1.d0,0.d0), &
scratch,Size(scratch,1), &
Fock_Matrix_AO_complex,Size(Fock_Matrix_AO_complex,1), &
(1.d0,0.d0), &
FPS_SPF_Matrix_AO_complex,Size(FPS_SPF_Matrix_AO_complex,1))
END_PROVIDER
BEGIN_PROVIDER [complex*16, FPS_SPF_Matrix_MO, (mo_num, mo_num)]
implicit none
begin_doc
! Commutator FPS - SPF in MO basis
end_doc
call ao_to_mo_complex(FPS_SPF_Matrix_AO_complex, size(FPS_SPF_Matrix_AO_complex,1), &
FPS_SPF_Matrix_MO_complex, size(FPS_SPF_Matrix_MO_complex,1))
END_PROVIDER
BEGIN_PROVIDER [ double precision, eigenvalues_Fock_matrix_AO_complex, (AO_num) ]
&BEGIN_PROVIDER [ complex*16, eigenvectors_Fock_matrix_AO_complex, (AO_num,AO_num) ]
!TODO: finish this provider; write provider for S_half_inv_complex
BEGIN_DOC
! Eigenvalues and eigenvectors of the Fock matrix over the AO basis
END_DOC
implicit none
double precision, allocatable :: scratch(:,:),work(:),Xt(:,:)
integer :: lwork,info
integer :: i,j
lwork = 3*AO_num - 1
allocate( &
scratch(AO_num,AO_num), &
work(lwork), &
Xt(AO_num,AO_num) &
)
! Calculate Xt
do i=1,AO_num
do j=1,AO_num
Xt(i,j) = S_half_inv(j,i)
enddo
enddo
! Calculate Fock matrix in orthogonal basis: F' = Xt.F.X
call dgemm('N','N',AO_num,AO_num,AO_num, &
1.d0, &
Fock_matrix_AO,size(Fock_matrix_AO,1), &
S_half_inv,size(S_half_inv,1), &
0.d0, &
eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1))
call dgemm('N','N',AO_num,AO_num,AO_num, &
1.d0, &
Xt,size(Xt,1), &
eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1), &
0.d0, &
scratch,size(scratch,1))
! Diagonalize F' to obtain eigenvectors in orthogonal basis C' and eigenvalues
call dsyev('V','U',AO_num, &
scratch,size(scratch,1), &
eigenvalues_Fock_matrix_AO, &
work,lwork,info)
if(info /= 0) then
print *, irp_here//' failed : ', info
stop 1
endif
! Back-transform eigenvectors: C =X.C'
call dgemm('N','N',AO_num,AO_num,AO_num, &
1.d0, &
S_half_inv,size(S_half_inv,1), &
scratch,size(scratch,1), &
0.d0, &
eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1))
END_PROVIDER

321
src/scf_utils/roothaan_hall_scf_complex.irp.f Normal file
View File

@ -0,0 +1,321 @@
subroutine Roothaan_Hall_SCF_complex
BEGIN_DOC
! Roothaan-Hall algorithm for SCF Hartree-Fock calculation
END_DOC
implicit none
double precision :: energy_SCF,energy_SCF_previous,Delta_energy_SCF
double precision :: max_error_DIIS,max_error_DIIS_alpha,max_error_DIIS_beta
complex*16, allocatable :: Fock_matrix_DIIS(:,:,:),error_matrix_DIIS(:,:,:)
integer :: iteration_SCF,dim_DIIS,index_dim_DIIS
integer :: i,j
logical, external :: qp_stop
complex*16, allocatable :: mo_coef_save(:,:)
PROVIDE ao_md5 mo_occ level_shift
allocate(mo_coef_save(ao_num,mo_num), &
Fock_matrix_DIIS (ao_num,ao_num,max_dim_DIIS), &
error_matrix_DIIS(ao_num,ao_num,max_dim_DIIS) &
)
call write_time(6)
print*,'Energy of the guess = ',SCF_energy
write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') &
'====','================','================','================','================'
write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') &
' N ', 'Energy ', 'Energy diff ', 'DIIS error ', 'Level shift '
write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') &
'====','================','================','================','================'
! Initialize energies and density matrices
energy_SCF_previous = SCF_energy
Delta_energy_SCF = 1.d0
iteration_SCF = 0
dim_DIIS = 0
max_error_DIIS = 1.d0
!
! Start of main SCF loop
!
PROVIDE FPS_SPF_matrix_AO_complex Fock_matrix_AO_complex
do while ( &
( (max_error_DIIS > threshold_DIIS_nonzero) .or. &
(dabs(Delta_energy_SCF) > thresh_SCF) &
) .and. (iteration_SCF < n_it_SCF_max) )
! Increment cycle number
iteration_SCF += 1
if(frozen_orb_scf)then
call initialize_mo_coef_begin_iteration
endif
! Current size of the DIIS space
dim_DIIS = min(dim_DIIS+1,max_dim_DIIS)
if (scf_algorithm == 'DIIS') then
! Store Fock and error matrices at each iteration
do j=1,ao_num
do i=1,ao_num
index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1
Fock_matrix_DIIS (i,j,index_dim_DIIS) = Fock_matrix_AO_complex(i,j)
error_matrix_DIIS(i,j,index_dim_DIIS) = FPS_SPF_matrix_AO_complex(i,j)
enddo
enddo
! Compute the extrapolated Fock matrix
call extrapolate_Fock_matrix_complex( &
error_matrix_DIIS,Fock_matrix_DIIS, &
Fock_matrix_AO_complex,size(Fock_matrix_AO_complex,1), &
iteration_SCF,dim_DIIS &
)
Fock_matrix_AO_alpha_complex = Fock_matrix_AO_complex*0.5d0
Fock_matrix_AO_beta_complex = Fock_matrix_AO_complex*0.5d0
TOUCH Fock_matrix_AO_alpha_complex Fock_matrix_AO_beta_complex
endif
mo_coef_complex = eigenvectors_fock_matrix_mo_complex
if(frozen_orb_scf)then
call reorder_core_orb
call initialize_mo_coef_begin_iteration
endif
TOUCH mo_coef_complex
! Calculate error vectors
max_error_DIIS = maxval(cdabs(FPS_SPF_Matrix_MO_complex))
! SCF energy
energy_SCF = SCF_energy
Delta_Energy_SCF = energy_SCF - energy_SCF_previous
if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0) ) then
Fock_matrix_AO_complex(1:ao_num,1:ao_num) = Fock_matrix_DIIS (1:ao_num,1:ao_num,index_dim_DIIS)
Fock_matrix_AO_alpha_complex = Fock_matrix_AO_complex*0.5d0
Fock_matrix_AO_beta_complex = Fock_matrix_AO_complex*0.5d0
TOUCH Fock_matrix_AO_alpha_complex Fock_matrix_AO_beta_complex
endif
double precision :: level_shift_save
level_shift_save = level_shift
mo_coef_save(1:ao_num,1:mo_num) = mo_coef_complex(1:ao_num,1:mo_num)
do while (Delta_energy_SCF > 0.d0)
mo_coef_complex(1:ao_num,1:mo_num) = mo_coef_save
if (level_shift <= .1d0) then
level_shift = 1.d0
else
level_shift = level_shift * 3.0d0
endif
TOUCH mo_coef_complex level_shift
mo_coef_complex(1:ao_num,1:mo_num) = eigenvectors_Fock_matrix_MO_complex(1:ao_num,1:mo_num)
if(frozen_orb_scf)then
call reorder_core_orb
call initialize_mo_coef_begin_iteration
endif
TOUCH mo_coef_complex
Delta_Energy_SCF = SCF_energy - energy_SCF_previous
energy_SCF = SCF_energy
if (level_shift-level_shift_save > 40.d0) then
level_shift = level_shift_save * 4.d0
SOFT_TOUCH level_shift
exit
endif
dim_DIIS=0
enddo
level_shift = level_shift * 0.5d0
SOFT_TOUCH level_shift
energy_SCF_previous = energy_SCF
! Print results at the end of each iteration
write(6,'(I4, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, I3)') &
iteration_SCF, energy_SCF, Delta_energy_SCF, max_error_DIIS, level_shift, dim_DIIS
if (Delta_energy_SCF < 0.d0) then
call save_mos
endif
if (qp_stop()) exit
enddo
if (iteration_SCF < n_it_SCF_max) then
mo_label = "Canonical"
endif
!
! End of Main SCF loop
!
write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') &
'====','================','================','================','================'
write(6,*)
if(.not.frozen_orb_scf)then
call mo_as_eigvectors_of_mo_matrix_complex(Fock_matrix_mo_complex,size(Fock_matrix_mo_complex,1),size(Fock_matrix_mo_complex,2),mo_label,1,.true.)
call save_mos
endif
call write_double(6, Energy_SCF, 'SCF energy')
call write_time(6)
end
subroutine extrapolate_Fock_matrix_complex( &
error_matrix_DIIS,Fock_matrix_DIIS, &
Fock_matrix_AO_,size_Fock_matrix_AO, &
iteration_SCF,dim_DIIS &
)
BEGIN_DOC
! Compute the extrapolated Fock matrix using the DIIS procedure
END_DOC
implicit none
complex*16,intent(in) :: Fock_matrix_DIIS(ao_num,ao_num,*),error_matrix_DIIS(ao_num,ao_num,*)
integer,intent(in) :: iteration_SCF, size_Fock_matrix_AO
complex*16,intent(inout):: Fock_matrix_AO_(size_Fock_matrix_AO,ao_num)
integer,intent(inout) :: dim_DIIS
double precision,allocatable :: B_matrix_DIIS(:,:),X_vector_DIIS(:)
double precision,allocatable :: C_vector_DIIS(:)
double precision :: accum_im, thr_im
complex*16,allocatable :: scratch(:,:)
integer :: i,j,k,i_DIIS,j_DIIS
thr_im = 1.0d-10
allocate( &
B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1), &
X_vector_DIIS(dim_DIIS+1), &
C_vector_DIIS(dim_DIIS+1), &
scratch(ao_num,ao_num) &
)
! Compute the matrices B and X
do j=1,dim_DIIS
do i=1,dim_DIIS
j_DIIS = mod(iteration_SCF-j,max_dim_DIIS)+1
i_DIIS = mod(iteration_SCF-i,max_dim_DIIS)+1
! Compute product of two errors vectors
call zgemm('N','N',ao_num,ao_num,ao_num, &
(1.d0,0.d0), &
error_matrix_DIIS(1,1,i_DIIS),size(error_matrix_DIIS,1), &
error_matrix_DIIS(1,1,j_DIIS),size(error_matrix_DIIS,1), &
(0.d0,0.d0), &
scratch,size(scratch,1))
! Compute Trace
B_matrix_DIIS(i,j) = 0.d0
accum_im = 0.d0
do k=1,ao_num
B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + dble(scratch(k,k))
accum_im = accum_im + dimag(scratch(k,k))
enddo
if (dabs(accum_im) .gt. thr_im) then
!stop 'problem with imaginary parts in DIIS B_matrix?'
print*, 'problem with imaginary parts in DIIS B_matrix?',accum_im
endif
enddo
enddo
! Pad B matrix and build the X matrix
do i=1,dim_DIIS
B_matrix_DIIS(i,dim_DIIS+1) = -1.d0
B_matrix_DIIS(dim_DIIS+1,i) = -1.d0
C_vector_DIIS(i) = 0.d0
enddo
B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1) = 0.d0
C_vector_DIIS(dim_DIIS+1) = -1.d0
! Solve the linear system C = B.X
integer :: info
integer,allocatable :: ipiv(:)
allocate( &
ipiv(dim_DIIS+1) &
)
double precision, allocatable :: AF(:,:)
allocate (AF(dim_DIIS+1,dim_DIIS+1))
double precision :: rcond, ferr, berr
integer :: iwork(dim_DIIS+1), lwork
call dsysvx('N','U',dim_DIIS+1,1, &
B_matrix_DIIS,size(B_matrix_DIIS,1), &
AF, size(AF,1), &
ipiv, &
C_vector_DIIS,size(C_vector_DIIS,1), &
X_vector_DIIS,size(X_vector_DIIS,1), &
rcond, &
ferr, &
berr, &
scratch,-1, &
iwork, &
info &
)
lwork = int(scratch(1,1))
deallocate(scratch)
allocate(scratch(lwork,1))
call dsysvx('N','U',dim_DIIS+1,1, &
B_matrix_DIIS,size(B_matrix_DIIS,1), &
AF, size(AF,1), &
ipiv, &
C_vector_DIIS,size(C_vector_DIIS,1), &
X_vector_DIIS,size(X_vector_DIIS,1), &
rcond, &
ferr, &
berr, &
scratch,size(scratch), &
iwork, &
info &
)
deallocate(scratch,ipiv)
if(info < 0) then
stop 'bug in DIIS'
endif
if (rcond > 1.d-12) then
! Compute extrapolated Fock matrix
!$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED) if (ao_num > 200)
do j=1,ao_num
do i=1,ao_num
Fock_matrix_AO_(i,j) = (0.d0,0.d0)
enddo
do k=1,dim_DIIS
do i=1,ao_num
Fock_matrix_AO_(i,j) = Fock_matrix_AO_(i,j) + &
X_vector_DIIS(k)*Fock_matrix_DIIS(i,j,dim_DIIS-k+1)
enddo
enddo
enddo
!$OMP END PARALLEL DO
else
dim_DIIS = 0
endif
end