9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 19:43:32 +01:00
This commit is contained in:
Kevin Gasperich 2020-03-20 12:22:10 -05:00
parent a0eb1d34db
commit d0fe9aad4f
15 changed files with 1193 additions and 283 deletions

View File

@ -209,6 +209,18 @@ BEGIN_PROVIDER [ complex*16, S_inv_complex,(ao_num,ao_num) ]
size(ao_overlap_complex,1),ao_num,ao_num,S_inv_complex,size(S_inv_complex,1)) size(ao_overlap_complex,1),ao_num,ao_num,S_inv_complex,size(S_inv_complex,1))
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ complex*16, S_inv_kpts,(ao_num_per_kpt,ao_num_per_kpt,kpt_num) ]
implicit none
BEGIN_DOC
! Inverse of the overlap matrix
END_DOC
integer :: k
do k=1,kpt_num
call get_pseudo_inverse_complex(ao_overlap_kpts(1,1,k), &
size(ao_overlap_kpts,1),ao_num_per_kpt,ao_num_per_kpt,S_inv_kpts(1,1,k),size(S_inv_kpts,1))
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ] BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ]
BEGIN_DOC BEGIN_DOC
@ -326,6 +338,66 @@ BEGIN_PROVIDER [ complex*16, S_half_inv_complex, (AO_num,AO_num) ]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ complex*16, S_half_inv_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_num) ]
BEGIN_DOC
! :math:`X = S^{-1/2}` obtained by SVD
END_DOC
implicit none
integer :: num_linear_dependencies
integer :: LDA, LDC
double precision, allocatable :: D(:)
complex*16, allocatable :: U(:,:),Vt(:,:)
integer :: info, i, j, k,kk
double precision, parameter :: threshold_overlap_AO_eigenvalues = 1.d-6
LDA = size(ao_overlap_kpts,1)
LDC = size(s_half_inv_kpts,1)
allocate( &
U(LDC,ao_num_per_kpt), &
Vt(LDA,ao_num_per_kpt), &
D(ao_num_per_kpt))
do kk=1,kpt_num
call svd_complex( &
ao_overlap_kpts(1,1,kk),LDA, &
U,LDC, &
D, &
Vt,LDA, &
ao_num_per_kpt,ao_num_per_kpt)
num_linear_dependencies = 0
do i=1,ao_num_per_kpt
print*,D(i)
if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then
D(i) = 0.d0
num_linear_dependencies += 1
else
ASSERT (D(i) > 0.d0)
D(i) = 1.d0/sqrt(D(i))
endif
do j=1,ao_num_per_kpt
S_half_inv_kpts(j,i,kk) = 0.d0
enddo
enddo
write(*,*) 'linear dependencies, k: ',num_linear_dependencies,', ',kk
do k=1,ao_num_per_kpt
if(D(k) /= 0.d0) then
do j=1,ao_num_per_kpt
do i=1,ao_num_per_kpt
S_half_inv_kpts(i,j,kk) = S_half_inv_kpts(i,j,kk) + U(i,k)*D(k)*Vt(k,j)
enddo
enddo
endif
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, S_half, (ao_num,ao_num) ] BEGIN_PROVIDER [ double precision, S_half, (ao_num,ao_num) ]
implicit none implicit none
@ -395,3 +467,39 @@ BEGIN_PROVIDER [ complex*16, S_half_complex, (ao_num,ao_num) ]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ complex*16, S_half_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_num) ]
implicit none
BEGIN_DOC
! :math:`S^{1/2}`
END_DOC
integer :: i,j,k,kk
complex*16, allocatable :: U(:,:)
complex*16, allocatable :: Vt(:,:)
double precision, allocatable :: D(:)
allocate(U(ao_num_per_kpt,ao_num_per_kpt),Vt(ao_num_per_kpt,ao_num_per_kpt),D(ao_num_per_kpt))
do kk=1,kpt_num
call svd_complex(ao_overlap_kpts(1,1,k),size(ao_overlap_kpts,1),U,size(U,1),D,Vt,size(Vt,1),ao_num_per_kpt,ao_num_per_kpt)
do i=1,ao_num_per_kpt
D(i) = dsqrt(D(i))
do j=1,ao_num_per_kpt
S_half_kpts(j,i,kk) = (0.d0,0.d0)
enddo
enddo
do k=1,ao_num_per_kpt
do j=1,ao_num_per_kpt
do i=1,ao_num_per_kpt
S_half_kpts(i,j,kk) = S_half_kpts(i,j,kk) + U(i,k)*D(k)*Vt(k,j)
enddo
enddo
enddo
enddo
deallocate(U,Vt,D)
END_PROVIDER

View File

@ -16,6 +16,15 @@ BEGIN_PROVIDER [ complex*16, mo_coef_begin_iteration_complex, (ao_num,mo_num) ]
END_DOC END_DOC
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ complex*16, mo_coef_begin_iteration_kpts, (ao_num_per_kpt,mo_num_per_kpt,kpt_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 subroutine initialize_mo_coef_begin_iteration
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -23,7 +32,8 @@ subroutine initialize_mo_coef_begin_iteration
! Initialize :c:data:`mo_coef_begin_iteration` to the current :c:data:`mo_coef` ! Initialize :c:data:`mo_coef_begin_iteration` to the current :c:data:`mo_coef`
END_DOC END_DOC
if (is_complex) then if (is_complex) then
mo_coef_begin_iteration_complex = mo_coef_complex !mo_coef_begin_iteration_complex = mo_coef_complex
mo_coef_begin_iteration_kpts = mo_coef_kpts
else else
mo_coef_begin_iteration = mo_coef mo_coef_begin_iteration = mo_coef
endif endif
@ -42,37 +52,71 @@ subroutine reorder_core_orb
integer :: i1,i2 integer :: i1,i2
if (is_complex) then if (is_complex) then
complex*16, allocatable :: accu_c(:) complex*16, allocatable :: accu_c(:)
allocate(accu(mo_num),accu_c(mo_num),index_core_orb(n_core_orb),iorder(mo_num)) !allocate(accu(mo_num),accu_c(mo_num),index_core_orb(n_core_orb),iorder(mo_num))
do i = 1, n_core_orb !do i = 1, n_core_orb
iorb = list_core(i) ! iorb = list_core(i)
do j = 1, mo_num ! do j = 1, mo_num
accu(j) = 0.d0 ! accu(j) = 0.d0
accu_c(j) = (0.d0,0.d0) ! accu_c(j) = (0.d0,0.d0)
iorder(j) = j ! iorder(j) = j
do k = 1, ao_num ! do k = 1, ao_num
do l = 1, ao_num ! do l = 1, ao_num
accu_c(j) += dconjg(mo_coef_begin_iteration_complex(k,iorb)) * & ! accu_c(j) += dconjg(mo_coef_begin_iteration_complex(k,iorb)) * &
mo_coef_complex(l,j) * ao_overlap_complex(k,l) ! mo_coef_complex(l,j) * ao_overlap_complex(k,l)
enddo ! enddo
enddo ! enddo
accu(j) = -cdabs(accu_c(j)) ! accu(j) = -cdabs(accu_c(j))
enddo ! enddo
call dsort(accu,iorder,mo_num) ! call dsort(accu,iorder,mo_num)
index_core_orb(i) = iorder(1) ! index_core_orb(i) = iorder(1)
enddo !enddo
complex*16 :: x_c !complex*16 :: x_c
do j = 1, n_core_orb !do j = 1, n_core_orb
i1 = list_core(j) ! i1 = list_core(j)
i2 = index_core_orb(j) ! i2 = index_core_orb(j)
do i=1,ao_num ! do i=1,ao_num
x_c = mo_coef_complex(i,i1) ! x_c = mo_coef_complex(i,i1)
mo_coef_complex(i,i1) = mo_coef_complex(i,i2) ! mo_coef_complex(i,i1) = mo_coef_complex(i,i2)
mo_coef_complex(i,i2) = x_c ! mo_coef_complex(i,i2) = x_c
enddo ! enddo
enddo !enddo
!call loc_cele_routine !!call loc_cele_routine
!deallocate(accu,accu_c,index_core_orb, iorder)
allocate(accu(mo_num_per_kpt),accu_c(mo_num_per_kpt),index_core_orb(n_core_orb),iorder(mo_num_per_kpt))
integer :: kk
do kk=1,kpt_num
do i = 1, n_core_orb_kpts(kk)
iorb = list_core_kpts(i,kk)
do j = 1, mo_num_per_kpt
accu(j) = 0.d0
accu_c(j) = (0.d0,0.d0)
iorder(j) = j
do k = 1, ao_num_per_kpt
do l = 1, ao_num_per_kpt
accu_c(j) += dconjg(mo_coef_begin_iteration_kpts(k,iorb,kk)) * &
mo_coef_kpts(l,j,kk) * ao_overlap_kpts(k,l,kk)
enddo
enddo
accu(j) = -cdabs(accu_c(j))
enddo
call dsort(accu,iorder,mo_num_per_kpt)
index_core_orb(i) = iorder(1)
enddo
complex*16 :: x_c
do j = 1, n_core_orb
i1 = list_core_kpts(j,kk)
i2 = index_core_orb(j)
do i=1,ao_num_per_kpt
x_c = mo_coef_kpts(i,i1,kk)
mo_coef_kpts(i,i1,kk) = mo_coef_kpts(i,i2,kk)
mo_coef_kpts(i,i2,kk) = x_c
enddo
enddo
!call loc_cele_routine
enddo
deallocate(accu,accu_c,index_core_orb, iorder) deallocate(accu,accu_c,index_core_orb, iorder)
else else
allocate(accu(mo_num),index_core_orb(n_core_orb),iorder(mo_num)) allocate(accu(mo_num),index_core_orb(n_core_orb),iorder(mo_num))

View File

@ -18,7 +18,7 @@ END_PROVIDER
BEGIN_DOC BEGIN_DOC
! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components. ! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components.
END_DOC END_DOC
integer :: i,j integer :: i,j,k
hf_energy = nuclear_repulsion hf_energy = nuclear_repulsion
hf_two_electron_energy = 0.d0 hf_two_electron_energy = 0.d0
hf_one_electron_energy = 0.d0 hf_one_electron_energy = 0.d0
@ -26,12 +26,14 @@ END_PROVIDER
complex*16 :: hf_1e_tmp, hf_2e_tmp complex*16 :: hf_1e_tmp, hf_2e_tmp
hf_1e_tmp = (0.d0,0.d0) hf_1e_tmp = (0.d0,0.d0)
hf_2e_tmp = (0.d0,0.d0) hf_2e_tmp = (0.d0,0.d0)
do j=1,ao_num do k=1,kpt_num
do i=1,ao_num do j=1,ao_num_per_kpt
hf_2e_tmp += 0.5d0 * ( ao_two_e_integral_alpha_complex(i,j) * scf_density_matrix_ao_alpha_complex(j,i) & do i=1,ao_num_per_kpt
+ao_two_e_integral_beta_complex(i,j) * scf_density_matrix_ao_beta_complex(j,i) ) hf_2e_tmp += 0.5d0 * ( ao_two_e_integral_alpha_kpts(i,j,k) * scf_density_matrix_ao_alpha_kpts(j,i,k) &
hf_1e_tmp += ao_one_e_integrals_complex(i,j) * (scf_density_matrix_ao_alpha_complex(j,i) & +ao_two_e_integral_beta_kpts(i,j,k) * scf_density_matrix_ao_beta_kpts(j,i,k) )
+ scf_density_matrix_ao_beta_complex (j,i) ) hf_1e_tmp += ao_one_e_integrals_kpts(i,j,k) * (scf_density_matrix_ao_alpha_kpts(j,i,k) &
+ scf_density_matrix_ao_beta_kpts (j,i,k) )
enddo
enddo enddo
enddo enddo
if (dabs(dimag(hf_2e_tmp)).gt.1.d-10) then if (dabs(dimag(hf_2e_tmp)).gt.1.d-10) then

View File

@ -102,7 +102,8 @@ subroutine run
mo_label = "Orthonormalized" mo_label = "Orthonormalized"
if (is_complex) then if (is_complex) then
call roothaan_hall_scf_complex !call roothaan_hall_scf_complex
call roothaan_hall_scf_kpts
else else
call roothaan_hall_scf call roothaan_hall_scf
endif endif

View File

@ -52,11 +52,11 @@ subroutine mo_as_eigvectors_of_mo_matrix_complex(matrix,n,m,label,sign,output)
enddo enddo
write (6,'(A)') '======== ================' write (6,'(A)') '======== ================'
write (6,'(A)') '' write (6,'(A)') ''
write (6,'(A)') 'Fock Matrix' !write (6,'(A)') 'Fock Matrix'
write (6,'(A)') '-----------' !write (6,'(A)') '-----------'
do i=1,n !do i=1,n
write(*,'(200(E24.15))') A(i,:) ! write(*,'(200(E24.15))') A(i,:)
enddo !enddo
endif endif
call zgemm('N','N',ao_num,m,m,(1.d0,0.d0),mo_coef_new,size(mo_coef_new,1),R,size(R,1),(0.d0,0.d0),mo_coef_complex,size(mo_coef_complex,1)) call zgemm('N','N',ao_num,m,m,(1.d0,0.d0),mo_coef_new,size(mo_coef_new,1),R,size(R,1),(0.d0,0.d0),mo_coef_complex,size(mo_coef_complex,1))
@ -302,11 +302,11 @@ subroutine mo_as_eigvectors_of_mo_matrix_kpts(matrix,n,m,nk,label,sign,output)
enddo enddo
write (6,'(A)') '======== ================' write (6,'(A)') '======== ================'
write (6,'(A)') '' write (6,'(A)') ''
write (6,'(A)') 'Fock Matrix' !write (6,'(A)') 'Fock Matrix'
write (6,'(A)') '-----------' !write (6,'(A)') '-----------'
do i=1,n !do i=1,n
write(*,'(200(E24.15))') A(i,:) ! write(*,'(200(E24.15))') A(i,:)
enddo !enddo
endif endif
call zgemm('N','N',ao_num_per_kpt,m,m,(1.d0,0.d0), & call zgemm('N','N',ao_num_per_kpt,m,m,(1.d0,0.d0), &

View File

@ -47,6 +47,18 @@ doc: Read/Write |MO| one-electron kinetic integrals from/to disk [ Write | Read
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: None default: None
[mo_integrals_overlap_kpts]
type: double precision
doc: Complex overlap integrals in |MO| basis set
size: (2,mo_basis.mo_num_per_kpt,mo_basis.mo_num_per_kpt,nuclei.kpt_num)
interface: ezfio
[io_mo_integrals_overlap]
type: Disk_access
doc: Read/Write |MO| one-electron overlap integrals from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[mo_integrals_pseudo] [mo_integrals_pseudo]
type: double precision type: double precision

View File

@ -74,3 +74,57 @@ BEGIN_PROVIDER [ complex*16, mo_overlap_complex,(mo_num,mo_num) ]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ complex*16, mo_overlap_kpts,(mo_num_per_kpt,mo_num_per_kpt,kpt_num) ]
implicit none
BEGIN_DOC
! Provider to check that the MOs are indeed orthonormal.
END_DOC
integer :: i,j,n,l,k
integer :: lmax
print *, 'Providing MO overlap integrals'
if (read_mo_integrals_overlap) then
call ezfio_get_mo_one_e_ints_mo_integrals_overlap_kpts(mo_overlap_kpts)
print *, 'MO overlap integrals read from disk'
else
print *, 'Providing MO overlap integrals from AO overlap integrals'
! call ao_to_mo_kpts( &
! ao_kinetic_integrals_kpts, &
! size(ao_kinetic_integrals_kpts,1), &
! mo_kinetic_integrals_kpts, &
! size(mo_kinetic_integrals_kpts,1) &
! )
!endif
lmax = (ao_num_per_kpt/4) * 4
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) &
!$OMP PRIVATE(i,j,n,l,k) &
!$OMP SHARED(mo_overlap_kpts,mo_coef_kpts,ao_overlap_kpts, &
!$OMP mo_num_per_kpt,ao_num_per_kpt,lmax,kpt_num)
do k=1,kpt_num
do j=1,mo_num_per_kpt
do i= 1,mo_num_per_kpt
mo_overlap_kpts(i,j,k) = (0.d0,0.d0)
do n = 1, lmax,4
do l = 1, ao_num_per_kpt
mo_overlap_kpts(i,j,k) = mo_overlap_kpts(i,j,k) + dconjg(mo_coef_kpts(l,i,k)) * &
( mo_coef_kpts(n ,j,k) * ao_overlap_kpts(l,n ,k) &
+ mo_coef_kpts(n+1,j,k) * ao_overlap_kpts(l,n+1,k) &
+ mo_coef_kpts(n+2,j,k) * ao_overlap_kpts(l,n+2,k) &
+ mo_coef_kpts(n+3,j,k) * ao_overlap_kpts(l,n+3,k) )
enddo
enddo
do n = lmax+1, ao_num_per_kpt
do l = 1, ao_num_per_kpt
mo_overlap_kpts(i,j,k) = mo_overlap_kpts(i,j,k) + mo_coef_kpts(n,j,k) * &
dconjg(mo_coef_kpts(l,i,k)) * ao_overlap_kpts(l,n,k)
enddo
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
endif
END_PROVIDER

View File

@ -1,12 +1,14 @@
subroutine orthonormalize_mos subroutine orthonormalize_mos
implicit none implicit none
integer :: m,p,s integer :: m,p,s,k
if (is_complex) then if (is_complex) then
m = size(mo_coef_complex,1) do k=1,kpt_num
p = size(mo_overlap_complex,1) m = size(mo_coef_kpts,1)
call ortho_lowdin_complex(mo_overlap_complex,p,mo_num,mo_coef_complex,m,ao_num) p = size(mo_overlap_kpts,1)
call ortho_lowdin_complex(mo_overlap_kpts(1,1,k),p,mo_num_per_kpt,mo_coef_kpts(1,1,k),m,ao_num_per_kpt)
enddo
mo_label = 'Orthonormalized' mo_label = 'Orthonormalized'
SOFT_TOUCH mo_coef_complex mo_label SOFT_TOUCH mo_coef_kpts mo_label
else else
m = size(mo_coef,1) m = size(mo_coef,1)
p = size(mo_overlap,1) p = size(mo_overlap,1)

View File

@ -72,8 +72,8 @@ BEGIN_PROVIDER [ complex*16, eigenvectors_Fock_matrix_mo_kpts, (ao_num_per_kpt,m
allocate (diag(mo_num_per_kpt) ) allocate (diag(mo_num_per_kpt) )
do k=1,kpt_num do k=1,kpt_num
do j=1,mo_num do j=1,mo_num_per_kpt
do i=1,mo_num do i=1,mo_num_per_kpt
!F(i,j) = fock_matrix_mo_complex(i,j) !F(i,j) = fock_matrix_mo_complex(i,j)
F(i,j) = fock_matrix_mo_kpts(i,j,k) F(i,j) = fock_matrix_mo_kpts(i,j,k)
enddo enddo

View File

@ -140,3 +140,155 @@ END_PROVIDER
deallocate(scratch) deallocate(scratch)
END_PROVIDER END_PROVIDER
!============================================!
! !
! kpts !
! !
!============================================!
BEGIN_PROVIDER [complex*16, FPS_SPF_Matrix_AO_kpts, (AO_num_per_kpt, AO_num_per_kpt,kpt_num)]
implicit none
BEGIN_DOC
! Commutator FPS - SPF
END_DOC
complex*16, allocatable :: scratch(:,:)
integer :: k
allocate( &
scratch(ao_num_per_kpt, ao_num_per_kpt) &
)
do k=1,kpt_num
! Compute FP
call zgemm('N','N',AO_num_per_kpt,AO_num_per_kpt,AO_num_per_kpt, &
(1.d0,0.d0), &
Fock_Matrix_AO_kpts(1,1,k),Size(Fock_Matrix_AO_kpts,1), &
SCF_Density_Matrix_AO_kpts(1,1,k),Size(SCF_Density_Matrix_AO_kpts,1), &
(0.d0,0.d0), &
scratch,Size(scratch,1))
! Compute FPS
call zgemm('N','N',AO_num_per_kpt,AO_num_per_kpt,AO_num_per_kpt, &
(1.d0,0.d0), &
scratch,Size(scratch,1), &
AO_Overlap_kpts(1,1,k),Size(AO_Overlap_kpts,1), &
(0.d0,0.d0), &
FPS_SPF_Matrix_AO_kpts(1,1,k),Size(FPS_SPF_Matrix_AO_kpts,1))
! Compute SP
call zgemm('N','N',AO_num_per_kpt,AO_num_per_kpt,AO_num_per_kpt, &
(1.d0,0.d0), &
AO_Overlap_kpts(1,1,k),Size(AO_Overlap_kpts,1), &
SCF_Density_Matrix_AO_kpts(1,1,k),Size(SCF_Density_Matrix_AO_kpts,1), &
(0.d0,0.d0), &
scratch,Size(scratch,1))
! Compute FPS - SPF
call zgemm('N','N',AO_num_per_kpt,AO_num_per_kpt,AO_num_per_kpt, &
(-1.d0,0.d0), &
scratch,Size(scratch,1), &
Fock_Matrix_AO_kpts(1,1,k),Size(Fock_Matrix_AO_kpts,1), &
(1.d0,0.d0), &
FPS_SPF_Matrix_AO_kpts(1,1,k),Size(FPS_SPF_Matrix_AO_kpts,1))
enddo
END_PROVIDER
BEGIN_PROVIDER [complex*16, FPS_SPF_Matrix_MO_kpts, (mo_num_per_kpt, mo_num_per_kpt,kpt_num)]
implicit none
begin_doc
! Commutator FPS - SPF in MO basis
end_doc
call ao_to_mo_kpts(FPS_SPF_Matrix_AO_kpts, size(FPS_SPF_Matrix_AO_kpts,1), &
FPS_SPF_Matrix_MO_kpts, size(FPS_SPF_Matrix_MO_kpts,1))
END_PROVIDER
BEGIN_PROVIDER [ double precision, eigenvalues_fock_matrix_ao_kpts, (ao_num_per_kpt,kpt_num) ]
&BEGIN_PROVIDER [ complex*16, eigenvectors_fock_matrix_ao_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_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 :: rwork(:)
integer :: lwork,info,lrwork
complex*16, allocatable :: scratch(:,:),Xt(:,:),work(:)
integer :: i,j,k
allocate( &
scratch(ao_num_per_kpt,ao_num_per_kpt), &
Xt(ao_num_per_kpt,ao_num_per_kpt) &
)
do k=1,kpt_num
! Calculate Xt
do i=1,ao_num_per_kpt
do j=1,ao_num_per_kpt
! Xt(i,j) = dconjg(s_half_inv_complex(j,i,k))
Xt(i,j) = dconjg(S_half_inv_kpts(j,i,k))
enddo
enddo
! Calculate Fock matrix in orthogonal basis: F' = Xt.F.X
call zgemm('N','N',ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt, &
(1.d0,0.d0), &
fock_matrix_ao_kpts(1,1,k),size(fock_matrix_ao_kpts,1), &
s_half_inv_kpts(1,1,k),size(s_half_inv_kpts,1), &
(0.d0,0.d0), &
eigenvectors_fock_matrix_ao_kpts(1,1,k), &
size(eigenvectors_fock_matrix_ao_kpts,1))
call zgemm('N','N',ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt, &
(1.d0,0.d0), &
Xt,size(Xt,1), &
eigenvectors_fock_matrix_ao_kpts(1,1,k), &
size(eigenvectors_fock_matrix_ao_kpts,1), &
(0.d0,0.d0), &
scratch,size(scratch,1))
! Diagonalize F' to obtain eigenvectors in orthogonal basis C' and eigenvalues
lrwork = 3*ao_num_per_kpt - 2
allocate(rwork(lrwork), work(1))
lwork = -1
call zheev('V','U',ao_num_per_kpt, &
scratch,size(scratch,1), &
eigenvalues_fock_matrix_ao_kpts(1,k), &
work,lwork,rwork,info)
lwork = int(work(1))
deallocate(work)
allocate(work(lwork))
call zheev('V','U',ao_num_per_kpt, &
scratch,size(scratch,1), &
eigenvalues_fock_matrix_ao_kpts(1,k), &
work,lwork,rwork,info)
if(info /= 0) then
print *, irp_here//' failed : ', info
stop 1
endif
deallocate(work,rwork)
! Back-transform eigenvectors: C =X.C'
call zgemm('N','N',ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt, &
(1.d0,0.d0), &
s_half_inv_kpts(1,1,k),size(s_half_inv_kpts,1), &
scratch,size(scratch,1), &
(0.d0,0.d0), &
eigenvectors_fock_matrix_ao_kpts(1,1,k), &
size(eigenvectors_fock_matrix_ao_kpts,1))
enddo
deallocate(scratch)
END_PROVIDER

View File

@ -157,15 +157,17 @@ BEGIN_PROVIDER [ double precision, SCF_energy ]
END_DOC END_DOC
SCF_energy = nuclear_repulsion SCF_energy = nuclear_repulsion
integer :: i,j integer :: i,j,k
if (is_complex) then if (is_complex) then
complex*16 :: scf_e_tmp complex*16 :: scf_e_tmp
scf_e_tmp = dcmplx(SCF_energy,0.d0) scf_e_tmp = dcmplx(SCF_energy,0.d0)
do j=1,ao_num do k=1,kpt_num
do i=1,ao_num do j=1,ao_num_per_kpt
scf_e_tmp += 0.5d0 * ( & do i=1,ao_num_per_kpt
(ao_one_e_integrals_complex(i,j) + Fock_matrix_ao_alpha_complex(i,j) ) * SCF_density_matrix_ao_alpha_complex(j,i) +& scf_e_tmp += 0.5d0 * ( &
(ao_one_e_integrals_complex(i,j) + Fock_matrix_ao_beta_complex (i,j) ) * SCF_density_matrix_ao_beta_complex (j,i) ) (ao_one_e_integrals_kpts(i,j,k) + Fock_matrix_ao_alpha_kpts(i,j,k) ) * SCF_density_matrix_ao_alpha_kpts(j,i,k) +&
(ao_one_e_integrals_kpts(i,j,k) + Fock_matrix_ao_beta_kpts (i,j,k) ) * SCF_density_matrix_ao_beta_kpts (j,i,k) )
enddo
enddo enddo
enddo enddo
!TODO: add check for imaginary part? (should be zero) !TODO: add check for imaginary part? (should be zero)

View File

@ -593,14 +593,14 @@ END_PROVIDER
j = jj(k2) j = jj(k2)
k = kk(k2) k = kk(k2)
l = ll(k2) l = ll(k2)
kpt_i = (i-1)/kpt_num +1 kpt_i = (i-1)/ao_num_per_kpt +1
kpt_j = (j-1)/kpt_num +1 kpt_j = (j-1)/ao_num_per_kpt +1
kpt_k = (k-1)/kpt_num +1 kpt_k = (k-1)/ao_num_per_kpt +1
kpt_l = (l-1)/kpt_num +1 kpt_l = (l-1)/ao_num_per_kpt +1
idx_i = mod(i,kpt_num) idx_i = mod(i-1,ao_num_per_kpt)+1
idx_j = mod(j,kpt_num) idx_j = mod(j-1,ao_num_per_kpt)+1
idx_k = mod(k,kpt_num) idx_k = mod(k-1,ao_num_per_kpt)+1
idx_l = mod(l,kpt_num) idx_l = mod(l-1,ao_num_per_kpt)+1
integral = i_sign(k2)*values(k1) !for klij and lkji, take complex conjugate integral = i_sign(k2)*values(k1) !for klij and lkji, take complex conjugate
!G_a(i,k) += D_{ab}(l,j)*(<ij|kl>) !G_a(i,k) += D_{ab}(l,j)*(<ij|kl>)
@ -611,7 +611,7 @@ END_PROVIDER
if (kpt_l.eq.kpt_j) then if (kpt_l.eq.kpt_j) then
c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral
if(kpt_i.ne.kpt_k) then if(kpt_i.ne.kpt_k) then
print*,'problem in ',irp_here print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l
stop 1 stop 1
endif endif
ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i) += c0 ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i) += c0
@ -620,7 +620,7 @@ END_PROVIDER
if (kpt_l.eq.kpt_i) then if (kpt_l.eq.kpt_i) then
if(kpt_j.ne.kpt_k) then if(kpt_j.ne.kpt_k) then
print*,'problem in ',irp_here print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l
stop 1 stop 1
endif endif
ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral
@ -636,20 +636,20 @@ END_PROVIDER
j = jj(k2) j = jj(k2)
k = kk(k2) k = kk(k2)
l = ll(k2) l = ll(k2)
kpt_i = (i-1)/kpt_num +1 kpt_i = (i-1)/ao_num_per_kpt +1
kpt_j = (j-1)/kpt_num +1 kpt_j = (j-1)/ao_num_per_kpt +1
kpt_k = (k-1)/kpt_num +1 kpt_k = (k-1)/ao_num_per_kpt +1
kpt_l = (l-1)/kpt_num +1 kpt_l = (l-1)/ao_num_per_kpt +1
idx_i = mod(i,kpt_num) idx_i = mod(i-1,ao_num_per_kpt)+1
idx_j = mod(j,kpt_num) idx_j = mod(j-1,ao_num_per_kpt)+1
idx_k = mod(k,kpt_num) idx_k = mod(k-1,ao_num_per_kpt)+1
idx_l = mod(l,kpt_num) idx_l = mod(l-1,ao_num_per_kpt)+1
integral = values(k1) integral = values(k1)
if (kpt_l.eq.kpt_j) then if (kpt_l.eq.kpt_j) then
c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral
if(kpt_i.ne.kpt_k) then if(kpt_i.ne.kpt_k) then
print*,'problem in ',irp_here print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l
stop 1 stop 1
endif endif
ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i) += c0 ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i) += c0
@ -658,7 +658,7 @@ END_PROVIDER
if (kpt_l.eq.kpt_i) then if (kpt_l.eq.kpt_i) then
if(kpt_j.ne.kpt_k) then if(kpt_j.ne.kpt_k) then
print*,'problem in ',irp_here print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l
stop 1 stop 1
endif endif
ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral
@ -714,14 +714,14 @@ END_PROVIDER
j = jj(k2) j = jj(k2)
k = kk(k2) k = kk(k2)
l = ll(k2) l = ll(k2)
kpt_i = (i-1)/kpt_num +1 kpt_i = (i-1)/ao_num_per_kpt +1
kpt_j = (j-1)/kpt_num +1 kpt_j = (j-1)/ao_num_per_kpt +1
kpt_k = (k-1)/kpt_num +1 kpt_k = (k-1)/ao_num_per_kpt +1
kpt_l = (l-1)/kpt_num +1 kpt_l = (l-1)/ao_num_per_kpt +1
idx_i = mod(i,kpt_num) idx_i = mod(i-1,ao_num_per_kpt)+1
idx_j = mod(j,kpt_num) idx_j = mod(j-1,ao_num_per_kpt)+1
idx_k = mod(k,kpt_num) idx_k = mod(k-1,ao_num_per_kpt)+1
idx_l = mod(l,kpt_num) idx_l = mod(l-1,ao_num_per_kpt)+1
integral = i_sign(k2)*values(k1) ! for klij and lkji, take conjugate integral = i_sign(k2)*values(k1) ! for klij and lkji, take conjugate
!G_a(i,k) += D_{ab}(l,j)*(<ij|kl>) !G_a(i,k) += D_{ab}(l,j)*(<ij|kl>)
@ -732,7 +732,7 @@ END_PROVIDER
if (kpt_l.eq.kpt_j) then if (kpt_l.eq.kpt_j) then
c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral
if(kpt_i.ne.kpt_k) then if(kpt_i.ne.kpt_k) then
print*,'problem in ',irp_here print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l
stop 1 stop 1
endif endif
ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i) += c0 ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i) += c0
@ -741,7 +741,7 @@ END_PROVIDER
if (kpt_l.eq.kpt_i) then if (kpt_l.eq.kpt_i) then
if(kpt_j.ne.kpt_k) then if(kpt_j.ne.kpt_k) then
print*,'problem in ',irp_here print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l
stop 1 stop 1
endif endif
ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral
@ -757,20 +757,20 @@ END_PROVIDER
j = jj(k2) j = jj(k2)
k = kk(k2) k = kk(k2)
l = ll(k2) l = ll(k2)
kpt_i = (i-1)/kpt_num +1 kpt_i = (i-1)/ao_num_per_kpt +1
kpt_j = (j-1)/kpt_num +1 kpt_j = (j-1)/ao_num_per_kpt +1
kpt_k = (k-1)/kpt_num +1 kpt_k = (k-1)/ao_num_per_kpt +1
kpt_l = (l-1)/kpt_num +1 kpt_l = (l-1)/ao_num_per_kpt +1
idx_i = mod(i,kpt_num) idx_i = mod(i-1,ao_num_per_kpt)+1
idx_j = mod(j,kpt_num) idx_j = mod(j-1,ao_num_per_kpt)+1
idx_k = mod(k,kpt_num) idx_k = mod(k-1,ao_num_per_kpt)+1
idx_l = mod(l,kpt_num) idx_l = mod(l-1,ao_num_per_kpt)+1
integral = values(k1) integral = values(k1)
if (kpt_l.eq.kpt_j) then if (kpt_l.eq.kpt_j) then
c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral
if(kpt_i.ne.kpt_k) then if(kpt_i.ne.kpt_k) then
print*,'problem in ',irp_here print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l
stop 1 stop 1
endif endif
ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i) += c0 ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i) += c0
@ -779,7 +779,7 @@ END_PROVIDER
if (kpt_l.eq.kpt_i) then if (kpt_l.eq.kpt_i) then
if(kpt_j.ne.kpt_k) then if(kpt_j.ne.kpt_k) then
print*,'problem in ',irp_here print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l
stop 1 stop 1
endif endif
ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral

View File

@ -319,3 +319,337 @@ END_DOC
endif endif
end end
!============================================!
! !
! kpts !
! !
!============================================!
subroutine Roothaan_Hall_SCF_kpts
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,k,kk
logical, external :: qp_stop
complex*16, allocatable :: mo_coef_save(:,:,:)
PROVIDE ao_md5 mo_occ level_shift
allocate(mo_coef_save(ao_num_per_kpt,mo_num_per_kpt,kpt_num), &
Fock_matrix_DIIS (ao_num_per_kpt,ao_num_per_kpt,max_dim_DIIS,kpt_num), &
error_matrix_DIIS(ao_num_per_kpt,ao_num_per_kpt,max_dim_DIIS,kpt_num) &
)
!todo: add kpt_num dim to diis mats? (3 or 4)
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
PROVIDE fps_spf_matrix_ao_kpts fock_matrix_ao_kpts
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
do kk=1,kpt_num
! Store Fock and error matrices at each iteration
do j=1,ao_num_per_kpt
do i=1,ao_num_per_kpt
index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1
Fock_matrix_DIIS (i,j,index_dim_DIIS,kk) = fock_matrix_ao_kpts(i,j,kk)
error_matrix_DIIS(i,j,index_dim_DIIS,kk) = fps_spf_matrix_ao_kpts(i,j,kk)
enddo
enddo
! Compute the extrapolated Fock matrix
call extrapolate_fock_matrix_kpts( &
error_matrix_DIIS(1,1,1,kk),Fock_matrix_DIIS(1,1,1,kk), &
Fock_matrix_AO_kpts(1,1,kk),size(Fock_matrix_AO_kpts,1), &
iteration_SCF,dim_DIIS &
)
enddo
Fock_matrix_AO_alpha_kpts = Fock_matrix_AO_kpts*0.5d0
Fock_matrix_AO_beta_kpts = Fock_matrix_AO_kpts*0.5d0
TOUCH Fock_matrix_AO_alpha_kpts Fock_matrix_AO_beta_kpts
endif
mo_coef_kpts = eigenvectors_fock_matrix_mo_kpts
if(frozen_orb_scf)then
call reorder_core_orb
call initialize_mo_coef_begin_iteration
endif
TOUCH mo_coef_kpts
! Calculate error vectors
max_error_DIIS = maxval(cdabs(FPS_SPF_Matrix_MO_kpts))
! SCF energy
! call print_debug_scf_complex
energy_SCF = scf_energy
Delta_Energy_SCF = energy_SCF - energy_SCF_previous
if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0) ) then
do kk=1,kpt_num
Fock_matrix_AO_kpts(1:ao_num_per_kpt,1:ao_num_per_kpt,kk) = &
Fock_matrix_DIIS (1:ao_num_per_kpt,1:ao_num_per_kpt,index_dim_DIIS,kk)
enddo
Fock_matrix_AO_alpha_kpts = Fock_matrix_AO_kpts*0.5d0
Fock_matrix_AO_beta_kpts = Fock_matrix_AO_kpts*0.5d0
TOUCH Fock_matrix_AO_alpha_kpts Fock_matrix_AO_beta_kpts
endif
double precision :: level_shift_save
level_shift_save = level_shift
mo_coef_save(1:ao_num_per_kpt,1:mo_num_per_kpt,1:kpt_num) = mo_coef_kpts(1:ao_num_per_kpt,1:mo_num_per_kpt,1:kpt_num)
do while (Delta_energy_SCF > 0.d0)
mo_coef_kpts(1:ao_num_per_kpt,1:mo_num_per_kpt,1:kpt_num) = mo_coef_save
if (level_shift <= .1d0) then
level_shift = 1.d0
else
level_shift = level_shift * 3.0d0
endif
TOUCH mo_coef_kpts level_shift
mo_coef_kpts(1:ao_num_per_kpt,1:mo_num_per_kpt,1:kpt_num) = &
eigenvectors_fock_matrix_mo_kpts(1:ao_num_per_kpt,1:mo_num_per_kpt,1:kpt_num)
if(frozen_orb_scf)then
call reorder_core_orb
call initialize_mo_coef_begin_iteration
endif
TOUCH mo_coef_kpts
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_kpts(Fock_matrix_mo_kpts,size(Fock_matrix_mo_kpts,1),size(Fock_matrix_mo_kpts,2),size(Fock_matrix_mo_kpts,3),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_kpts( &
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_per_kpt,ao_num_per_kpt,*),error_matrix_DIIS(ao_num_per_kpt,ao_num_per_kpt,*)
integer,intent(in) :: iteration_SCF, size_Fock_matrix_AO
complex*16,intent(inout):: Fock_matrix_AO_(size_Fock_matrix_AO,ao_num_per_kpt)
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_per_kpt,ao_num_per_kpt,ao_num_per_kpt, &
(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_per_kpt
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
deallocate(scratch)
! 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(:,:),scratch_d1(:)
allocate (AF(dim_DIIS+1,dim_DIIS+1),scratch_d1(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_d1,-1, &
iwork, &
info &
)
lwork = int(scratch_d1(1))
deallocate(scratch_d1)
allocate(scratch_d1(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_d1,size(scratch_d1), &
iwork, &
info &
)
deallocate(scratch_d1,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_per_kpt > 200)
do j=1,ao_num_per_kpt
do i=1,ao_num_per_kpt
Fock_matrix_AO_(i,j) = (0.d0,0.d0)
enddo
do k=1,dim_DIIS
do i=1,ao_num_per_kpt
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

View File

@ -690,7 +690,7 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
qph5.create_dataset('mo_basis/mo_coef_kpts_real',data=mo_k.real) qph5.create_dataset('mo_basis/mo_coef_kpts_real',data=mo_k.real)
qph5.create_dataset('mo_basis/mo_coef_kpts_imag',data=mo_k.imag) qph5.create_dataset('mo_basis/mo_coef_kpts_imag',data=mo_k.imag)
qph5.create_dataset('mo_basis/mo_coef_complex',data=mo_coef_blocked_f.view(dtype=np.float64).reshape((Nk*nmo,Nk*nao,2))) qph5.create_dataset('mo_basis/mo_coef_complex',data=mo_coef_blocked_f.view(dtype=np.float64).reshape((Nk*nmo,Nk*nao,2)))
qph5.create_dataset('mo_basis/mo_coef_complex_kpts',data=mo_coef_f.view(dtype=np.float64).reshape((Nk,nmo,nao,2))) qph5.create_dataset('mo_basis/mo_coef_kpts',data=mo_coef_f.view(dtype=np.float64).reshape((Nk,nmo,nao,2)))
print_kpts_unblocked(mo_k,'C.qp',mo_coef_threshold) print_kpts_unblocked(mo_k,'C.qp',mo_coef_threshold)
@ -729,6 +729,9 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
qph5.create_dataset('ao_one_e_ints/ao_integrals_kinetic',data=kin_ao_blocked_f.view(dtype=np.float64).reshape((Nk*nao,Nk*nao,2))) qph5.create_dataset('ao_one_e_ints/ao_integrals_kinetic',data=kin_ao_blocked_f.view(dtype=np.float64).reshape((Nk*nao,Nk*nao,2)))
qph5.create_dataset('ao_one_e_ints/ao_integrals_overlap',data=ovlp_ao_blocked_f.view(dtype=np.float64).reshape((Nk*nao,Nk*nao,2))) qph5.create_dataset('ao_one_e_ints/ao_integrals_overlap',data=ovlp_ao_blocked_f.view(dtype=np.float64).reshape((Nk*nao,Nk*nao,2)))
qph5.create_dataset('ao_one_e_ints/ao_integrals_n_e', data=ne_ao_blocked_f.view(dtype=np.float64).reshape((Nk*nao,Nk*nao,2))) qph5.create_dataset('ao_one_e_ints/ao_integrals_n_e', data=ne_ao_blocked_f.view(dtype=np.float64).reshape((Nk*nao,Nk*nao,2)))
qph5.create_dataset('ao_one_e_ints/ao_integrals_kinetic_kpts',data=kin_ao_f.view(dtype=np.float64).reshape((Nk,nao,nao,2)))
qph5.create_dataset('ao_one_e_ints/ao_integrals_overlap_kpts',data=ovlp_ao_f.view(dtype=np.float64).reshape((Nk,nao,nao,2)))
qph5.create_dataset('ao_one_e_ints/ao_integrals_n_e_kpts', data=ne_ao_f.view(dtype=np.float64).reshape((Nk,nao,nao,2)))
for fname,ints in zip(('S.qp','V.qp','T.qp'), for fname,ints in zip(('S.qp','V.qp','T.qp'),
(ovlp_ao, ne_ao, kin_ao)): (ovlp_ao, ne_ao, kin_ao)):
@ -757,9 +760,14 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
qph5.create_dataset('mo_one_e_ints/mo_integrals_overlap_imag',data=ovlp_mo_blocked.imag) qph5.create_dataset('mo_one_e_ints/mo_integrals_overlap_imag',data=ovlp_mo_blocked.imag)
qph5.create_dataset('mo_one_e_ints/mo_integrals_n_e_real', data=ne_mo_blocked.real) qph5.create_dataset('mo_one_e_ints/mo_integrals_n_e_real', data=ne_mo_blocked.real)
qph5.create_dataset('mo_one_e_ints/mo_integrals_n_e_imag', data=ne_mo_blocked.imag) qph5.create_dataset('mo_one_e_ints/mo_integrals_n_e_imag', data=ne_mo_blocked.imag)
qph5.create_dataset('mo_one_e_ints/mo_integrals_kinetic',data=kin_mo_blocked_f.view(dtype=np.float64).reshape((Nk*nmo,Nk*nmo,2))) qph5.create_dataset('mo_one_e_ints/mo_integrals_kinetic',data=kin_mo_blocked_f.view(dtype=np.float64).reshape((Nk*nmo,Nk*nmo,2)))
qph5.create_dataset('mo_one_e_ints/mo_integrals_overlap',data=ovlp_mo_blocked_f.view(dtype=np.float64).reshape((Nk*nmo,Nk*nmo,2))) qph5.create_dataset('mo_one_e_ints/mo_integrals_overlap',data=ovlp_mo_blocked_f.view(dtype=np.float64).reshape((Nk*nmo,Nk*nmo,2)))
qph5.create_dataset('mo_one_e_ints/mo_integrals_n_e', data=ne_mo_blocked_f.view(dtype=np.float64).reshape((Nk*nmo,Nk*nmo,2))) qph5.create_dataset('mo_one_e_ints/mo_integrals_n_e', data=ne_mo_blocked_f.view(dtype=np.float64).reshape((Nk*nmo,Nk*nmo,2)))
qph5.create_dataset('mo_one_e_ints/mo_integrals_kinetic_kpts',data=kin_mo_f.view(dtype=np.float64).reshape((Nk,nmo,nmo,2)))
qph5.create_dataset('mo_one_e_ints/mo_integrals_overlap_kpts',data=ovlp_mo_f.view(dtype=np.float64).reshape((Nk,nmo,nmo,2)))
qph5.create_dataset('mo_one_e_ints/mo_integrals_n_e_kpts', data=ne_mo_f.view(dtype=np.float64).reshape((Nk,nmo,nmo,2)))
for fname,ints in zip(('S.mo.qp','V.mo.qp','T.mo.qp'), for fname,ints in zip(('S.mo.qp','V.mo.qp','T.mo.qp'),
(ovlp_mo, ne_mo, kin_mo)): (ovlp_mo, ne_mo, kin_mo)):
print_kpts_unblocked_upper(ints,fname,thresh_mono) print_kpts_unblocked_upper(ints,fname,thresh_mono)

View File

@ -4,198 +4,389 @@ import h5py
import sys import sys
import numpy as np import numpy as np
filename = sys.argv[1] fname = sys.argv[1]
qph5path = sys.argv[2] qph5name = sys.argv[2]
ezfio.set_file(filename)
#qph5=h5py.File(qph5path,'r') #qph5=h5py.File(qph5path,'r')
def convert_kpts(filename,qph5path):
ezfio.set_nuclei_is_complex(True) ezfio.set_file(filename)
ezfio.set_nuclei_is_complex(True)
with h5py.File(qph5path,'r') as qph5:
kpt_num = qph5['nuclei'].attrs['kpt_num']
nucl_num = qph5['nuclei'].attrs['nucl_num']
ao_num = qph5['ao_basis'].attrs['ao_num']
mo_num = qph5['mo_basis'].attrs['mo_num']
elec_alpha_num = qph5['electrons'].attrs['elec_alpha_num']
elec_beta_num = qph5['electrons'].attrs['elec_beta_num']
ezfio.set_nuclei_kpt_num(kpt_num)
kpt_pair_num = (kpt_num*kpt_num + kpt_num)//2
ezfio.set_nuclei_kpt_pair_num(kpt_pair_num)
# don't multiply nuclei by kpt_num
# work in k-space, not in equivalent supercell
nucl_num_per_kpt = nucl_num
ezfio.set_nuclei_nucl_num(nucl_num_per_kpt)
# these are totals (kpt_num * num_per_kpt)
# need to change if we want to truncate orbital space within pyscf
ezfio.set_ao_basis_ao_num(ao_num)
ezfio.set_mo_basis_mo_num(mo_num)
ezfio.electrons_elec_alpha_num = elec_alpha_num
ezfio.electrons_elec_beta_num = elec_beta_num
##ao_num = mo_num
##Important !
#import math
#nelec_per_kpt = num_elec // n_kpts
#nelec_alpha_per_kpt = int(math.ceil(nelec_per_kpt / 2.))
#nelec_beta_per_kpt = int(math.floor(nelec_per_kpt / 2.))
#
#ezfio.electrons_elec_alpha_num = int(nelec_alpha_per_kpt * n_kpts)
#ezfio.electrons_elec_beta_num = int(nelec_beta_per_kpt * n_kpts)
#ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.))
#ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.))
#ezfio.set_utils_num_kpts(n_kpts)
#ezfio.set_integrals_bielec_df_num(n_aux)
#(old)Important
#ezfio.set_nuclei_nucl_num(nucl_num)
#ezfio.set_nuclei_nucl_charge([0.]*nucl_num)
#ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num )
#ezfio.set_nuclei_nucl_label( ['He'] * nucl_num )
with h5py.File(qph5path,'r') as qph5:
nucl_charge=qph5['nuclei/nucl_charge'][()].tolist()
nucl_coord=qph5['nuclei/nucl_coord'][()].T.tolist()
nucl_label=qph5['nuclei/nucl_label'][()].tolist()
nuclear_repulsion = qph5['nuclei'].attrs['nuclear_repulsion']
ezfio.set_nuclei_nucl_charge(nucl_charge)
ezfio.set_nuclei_nucl_coord(nucl_coord)
ezfio.set_nuclei_nucl_label(nucl_label)
ezfio.set_nuclei_io_nuclear_repulsion('Read')
ezfio.set_nuclei_nuclear_repulsion(nuclear_repulsion)
##########################################
# #
# Basis #
# #
##########################################
with h5py.File(qph5path,'r') as qph5:
ezfio.set_ao_basis_ao_basis(qph5['ao_basis'].attrs['ao_basis'])
ezfio.set_ao_basis_ao_nucl(qph5['ao_basis/ao_nucl'][()].tolist())
#Just need one (can clean this up later)
ao_prim_num_max = 5
d = [ [0] *ao_prim_num_max]*ao_num
ezfio.set_ao_basis_ao_prim_num([ao_prim_num_max]*ao_num)
ezfio.set_ao_basis_ao_power(d)
ezfio.set_ao_basis_ao_coef(d)
ezfio.set_ao_basis_ao_expo(d)
##########################################
# #
# MO Coef #
# #
##########################################
with h5py.File(qph5path,'r') as qph5:
with h5py.File(qph5path,'r') as qph5: kpt_num = qph5['nuclei'].attrs['kpt_num']
mo_coef_reim = qph5['mo_basis/mo_coef_complex'][()].tolist() nucl_num = qph5['nuclei'].attrs['nucl_num']
ezfio.set_mo_basis_mo_coef_complex(mo_coef_reim) ao_num = qph5['ao_basis'].attrs['ao_num']
#maybe fix qp so we don't need this? mo_num = qph5['mo_basis'].attrs['mo_num']
#ezfio.set_mo_basis_mo_coef([[i for i in range(mo_num)] * ao_num]) elec_alpha_num = qph5['electrons'].attrs['elec_alpha_num']
elec_beta_num = qph5['electrons'].attrs['elec_beta_num']
##########################################
# #
# Integrals Mono #
# #
##########################################
with h5py.File(qph5path,'r') as qph5: ezfio.set_nuclei_kpt_num(kpt_num)
if 'ao_one_e_ints' in qph5.keys(): kpt_pair_num = (kpt_num*kpt_num + kpt_num)//2
kin_ao_reim=qph5['ao_one_e_ints/ao_integrals_kinetic'][()].tolist() ezfio.set_nuclei_kpt_pair_num(kpt_pair_num)
ovlp_ao_reim=qph5['ao_one_e_ints/ao_integrals_overlap'][()].tolist()
ne_ao_reim=qph5['ao_one_e_ints/ao_integrals_n_e'][()].tolist() # don't multiply nuclei by kpt_num
# work in k-space, not in equivalent supercell
ezfio.set_ao_one_e_ints_ao_integrals_kinetic_complex(kin_ao_reim) nucl_num_per_kpt = nucl_num
ezfio.set_ao_one_e_ints_ao_integrals_overlap_complex(ovlp_ao_reim) ezfio.set_nuclei_nucl_num(nucl_num_per_kpt)
ezfio.set_ao_one_e_ints_ao_integrals_n_e_complex(ne_ao_reim)
# these are totals (kpt_num * num_per_kpt)
# need to change if we want to truncate orbital space within pyscf
ezfio.set_ao_basis_ao_num(ao_num)
ezfio.set_mo_basis_mo_num(mo_num)
ezfio.electrons_elec_alpha_num = elec_alpha_num
ezfio.electrons_elec_beta_num = elec_beta_num
##ao_num = mo_num
##Important !
#import math
#nelec_per_kpt = num_elec // n_kpts
#nelec_alpha_per_kpt = int(math.ceil(nelec_per_kpt / 2.))
#nelec_beta_per_kpt = int(math.floor(nelec_per_kpt / 2.))
#
#ezfio.electrons_elec_alpha_num = int(nelec_alpha_per_kpt * n_kpts)
#ezfio.electrons_elec_beta_num = int(nelec_beta_per_kpt * n_kpts)
#ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.))
#ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.))
#ezfio.set_utils_num_kpts(n_kpts)
#ezfio.set_integrals_bielec_df_num(n_aux)
#(old)Important
#ezfio.set_nuclei_nucl_num(nucl_num)
#ezfio.set_nuclei_nucl_charge([0.]*nucl_num)
#ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num )
#ezfio.set_nuclei_nucl_label( ['He'] * nucl_num )
with h5py.File(qph5path,'r') as qph5:
nucl_charge=qph5['nuclei/nucl_charge'][()].tolist()
nucl_coord=qph5['nuclei/nucl_coord'][()].T.tolist()
nucl_label=qph5['nuclei/nucl_label'][()].tolist()
nuclear_repulsion = qph5['nuclei'].attrs['nuclear_repulsion']
ezfio.set_nuclei_nucl_charge(nucl_charge)
ezfio.set_nuclei_nucl_coord(nucl_coord)
ezfio.set_nuclei_nucl_label(nucl_label)
ezfio.set_nuclei_io_nuclear_repulsion('Read')
ezfio.set_nuclei_nuclear_repulsion(nuclear_repulsion)
##########################################
# #
# Basis #
# #
##########################################
with h5py.File(qph5path,'r') as qph5:
ezfio.set_ao_basis_ao_basis(qph5['ao_basis'].attrs['ao_basis'])
ezfio.set_ao_basis_ao_nucl(qph5['ao_basis/ao_nucl'][()].tolist())
#Just need one (can clean this up later)
ao_prim_num_max = 5
d = [ [0] *ao_prim_num_max]*ao_num
ezfio.set_ao_basis_ao_prim_num([ao_prim_num_max]*ao_num)
ezfio.set_ao_basis_ao_power(d)
ezfio.set_ao_basis_ao_coef(d)
ezfio.set_ao_basis_ao_expo(d)
##########################################
# #
# MO Coef #
# #
##########################################
ezfio.set_ao_one_e_ints_io_ao_integrals_kinetic('Read')
ezfio.set_ao_one_e_ints_io_ao_integrals_overlap('Read')
ezfio.set_ao_one_e_ints_io_ao_integrals_n_e('Read')
with h5py.File(qph5path,'r') as qph5: with h5py.File(qph5path,'r') as qph5:
if 'mo_one_e_ints' in qph5.keys(): mo_coef_kpts = qph5['mo_basis/mo_coef_kpts'][()].tolist()
kin_mo_reim=qph5['mo_one_e_ints/mo_integrals_kinetic'][()].tolist() mo_coef_cplx = qph5['mo_basis/mo_coef_complex'][()].tolist()
#ovlp_mo_reim=qph5['mo_one_e_ints/mo_integrals_overlap'][()].tolist() ezfio.set_mo_basis_mo_coef_kpts(mo_coef_kpts)
ne_mo_reim=qph5['mo_one_e_ints/mo_integrals_n_e'][()].tolist() ezfio.set_mo_basis_mo_coef_complex(mo_coef_cplx)
#maybe fix qp so we don't need this?
ezfio.set_mo_one_e_ints_mo_integrals_kinetic_complex(kin_mo_reim) #ezfio.set_mo_basis_mo_coef([[i for i in range(mo_num)] * ao_num])
#ezfio.set_mo_one_e_ints_mo_integrals_overlap_complex(ovlp_mo_reim)
#ezfio.set_mo_one_e_ints_mo_integrals_n_e_complex(ne_mo_reim)
ezfio.set_mo_one_e_ints_mo_integrals_e_n_complex(ne_mo_reim) ##########################################
# #
# Integrals Mono #
# #
##########################################
ezfio.set_mo_one_e_ints_io_mo_integrals_kinetic('Read') with h5py.File(qph5path,'r') as qph5:
#ezfio.set_mo_one_e_ints_io_mo_integrals_overlap('Read') if 'ao_one_e_ints' in qph5.keys():
#ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read') kin_ao_reim=qph5['ao_one_e_ints/ao_integrals_kinetic_kpts'][()].tolist()
ezfio.set_mo_one_e_ints_io_mo_integrals_e_n('Read') ovlp_ao_reim=qph5['ao_one_e_ints/ao_integrals_overlap_kpts'][()].tolist()
ne_ao_reim=qph5['ao_one_e_ints/ao_integrals_n_e_kpts'][()].tolist()
##########################################
# # ezfio.set_ao_one_e_ints_ao_integrals_kinetic_kpts(kin_ao_reim)
# k-points # ezfio.set_ao_one_e_ints_ao_integrals_overlap_kpts(ovlp_ao_reim)
# # ezfio.set_ao_one_e_ints_ao_integrals_n_e_kpts(ne_ao_reim)
##########################################
ezfio.set_ao_one_e_ints_io_ao_integrals_kinetic('Read')
with h5py.File(qph5path,'r') as qph5: ezfio.set_ao_one_e_ints_io_ao_integrals_overlap('Read')
kconserv = qph5['nuclei/kconserv'][()].tolist() ezfio.set_ao_one_e_ints_io_ao_integrals_n_e('Read')
ezfio.set_nuclei_kconserv(kconserv)
ezfio.set_nuclei_io_kconserv('Read') with h5py.File(qph5path,'r') as qph5:
if 'mo_one_e_ints' in qph5.keys():
########################################## kin_mo_reim=qph5['mo_one_e_ints/mo_integrals_kinetic_kpts'][()].tolist()
# # ovlp_mo_reim=qph5['mo_one_e_ints/mo_integrals_overlap'][()].tolist()
# Integrals Bi # ne_mo_reim=qph5['mo_one_e_ints/mo_integrals_n_e_kpts'][()].tolist()
# #
########################################## ezfio.set_mo_one_e_ints_mo_integrals_kinetic_kpts(kin_mo_reim)
ezfio.set_mo_one_e_ints_mo_integrals_overlap_kpts(ovlp_mo_reim)
# should this be in ao_basis? ao_two_e_ints? #ezfio.set_mo_one_e_ints_mo_integrals_n_e_complex(ne_mo_reim)
with h5py.File(qph5path,'r') as qph5: ezfio.set_mo_one_e_ints_mo_integrals_e_n_kpts(ne_mo_reim)
if 'ao_two_e_ints' in qph5.keys():
df_num = qph5['ao_two_e_ints'].attrs['df_num'] ezfio.set_mo_one_e_ints_io_mo_integrals_kinetic('Read')
ezfio.set_ao_two_e_ints_df_num(df_num) ezfio.set_mo_one_e_ints_io_mo_integrals_overlap('Read')
if 'df_ao_integrals' in qph5['ao_two_e_ints'].keys(): #ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read')
# dfao_re0=qph5['ao_two_e_ints/df_ao_integrals_real'][()].transpose((3,2,1,0)) ezfio.set_mo_one_e_ints_io_mo_integrals_e_n('Read')
# dfao_im0=qph5['ao_two_e_ints/df_ao_integrals_imag'][()].transpose((3,2,1,0))
# dfao_cmplx0 = np.stack((dfao_re0,dfao_im0),axis=-1).tolist() ##########################################
# ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_cmplx0) # #
dfao_reim=qph5['ao_two_e_ints/df_ao_integrals'][()].tolist() # k-points #
ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_reim) # #
ezfio.set_ao_two_e_ints_io_df_ao_integrals('Read') ##########################################
if 'mo_two_e_ints' in qph5.keys(): with h5py.File(qph5path,'r') as qph5:
df_num = qph5['ao_two_e_ints'].attrs['df_num'] kconserv = qph5['nuclei/kconserv'][()].tolist()
ezfio.set_ao_two_e_ints_df_num(df_num)
# dfmo_re0=qph5['mo_two_e_ints/df_mo_integrals_real'][()].transpose((3,2,1,0)) ezfio.set_nuclei_kconserv(kconserv)
# dfmo_im0=qph5['mo_two_e_ints/df_mo_integrals_imag'][()].transpose((3,2,1,0)) ezfio.set_nuclei_io_kconserv('Read')
# dfmo_cmplx0 = np.stack((dfmo_re0,dfmo_im0),axis=-1).tolist()
# ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_cmplx0) ##########################################
dfmo_reim=qph5['mo_two_e_ints/df_mo_integrals'][()].tolist() # #
ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_reim) # Integrals Bi #
ezfio.set_mo_two_e_ints_io_df_mo_integrals('Read') # #
##########################################
# should this be in ao_basis? ao_two_e_ints?
with h5py.File(qph5path,'r') as qph5:
if 'ao_two_e_ints' in qph5.keys():
df_num = qph5['ao_two_e_ints'].attrs['df_num']
ezfio.set_ao_two_e_ints_df_num(df_num)
if 'df_ao_integrals' in qph5['ao_two_e_ints'].keys():
# dfao_re0=qph5['ao_two_e_ints/df_ao_integrals_real'][()].transpose((3,2,1,0))
# dfao_im0=qph5['ao_two_e_ints/df_ao_integrals_imag'][()].transpose((3,2,1,0))
# dfao_cmplx0 = np.stack((dfao_re0,dfao_im0),axis=-1).tolist()
# ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_cmplx0)
dfao_reim=qph5['ao_two_e_ints/df_ao_integrals'][()].tolist()
ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_reim)
ezfio.set_ao_two_e_ints_io_df_ao_integrals('Read')
if 'mo_two_e_ints' in qph5.keys():
df_num = qph5['ao_two_e_ints'].attrs['df_num']
ezfio.set_ao_two_e_ints_df_num(df_num)
# dfmo_re0=qph5['mo_two_e_ints/df_mo_integrals_real'][()].transpose((3,2,1,0))
# dfmo_im0=qph5['mo_two_e_ints/df_mo_integrals_imag'][()].transpose((3,2,1,0))
# dfmo_cmplx0 = np.stack((dfmo_re0,dfmo_im0),axis=-1).tolist()
# ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_cmplx0)
dfmo_reim=qph5['mo_two_e_ints/df_mo_integrals'][()].tolist()
ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_reim)
ezfio.set_mo_two_e_ints_io_df_mo_integrals('Read')
return
def convert_cplx(filename,qph5path):
ezfio.set_file(filename)
ezfio.set_nuclei_is_complex(True)
with h5py.File(qph5path,'r') as qph5:
kpt_num = qph5['nuclei'].attrs['kpt_num']
nucl_num = qph5['nuclei'].attrs['nucl_num']
ao_num = qph5['ao_basis'].attrs['ao_num']
mo_num = qph5['mo_basis'].attrs['mo_num']
elec_alpha_num = qph5['electrons'].attrs['elec_alpha_num']
elec_beta_num = qph5['electrons'].attrs['elec_beta_num']
ezfio.set_nuclei_kpt_num(kpt_num)
kpt_pair_num = (kpt_num*kpt_num + kpt_num)//2
ezfio.set_nuclei_kpt_pair_num(kpt_pair_num)
# don't multiply nuclei by kpt_num
# work in k-space, not in equivalent supercell
nucl_num_per_kpt = nucl_num
ezfio.set_nuclei_nucl_num(nucl_num_per_kpt)
# these are totals (kpt_num * num_per_kpt)
# need to change if we want to truncate orbital space within pyscf
ezfio.set_ao_basis_ao_num(ao_num)
ezfio.set_mo_basis_mo_num(mo_num)
ezfio.electrons_elec_alpha_num = elec_alpha_num
ezfio.electrons_elec_beta_num = elec_beta_num
##ao_num = mo_num
##Important !
#import math
#nelec_per_kpt = num_elec // n_kpts
#nelec_alpha_per_kpt = int(math.ceil(nelec_per_kpt / 2.))
#nelec_beta_per_kpt = int(math.floor(nelec_per_kpt / 2.))
#
#ezfio.electrons_elec_alpha_num = int(nelec_alpha_per_kpt * n_kpts)
#ezfio.electrons_elec_beta_num = int(nelec_beta_per_kpt * n_kpts)
#ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.))
#ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.))
#ezfio.set_utils_num_kpts(n_kpts)
#ezfio.set_integrals_bielec_df_num(n_aux)
#(old)Important
#ezfio.set_nuclei_nucl_num(nucl_num)
#ezfio.set_nuclei_nucl_charge([0.]*nucl_num)
#ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num )
#ezfio.set_nuclei_nucl_label( ['He'] * nucl_num )
with h5py.File(qph5path,'r') as qph5:
nucl_charge=qph5['nuclei/nucl_charge'][()].tolist()
nucl_coord=qph5['nuclei/nucl_coord'][()].T.tolist()
nucl_label=qph5['nuclei/nucl_label'][()].tolist()
nuclear_repulsion = qph5['nuclei'].attrs['nuclear_repulsion']
ezfio.set_nuclei_nucl_charge(nucl_charge)
ezfio.set_nuclei_nucl_coord(nucl_coord)
ezfio.set_nuclei_nucl_label(nucl_label)
ezfio.set_nuclei_io_nuclear_repulsion('Read')
ezfio.set_nuclei_nuclear_repulsion(nuclear_repulsion)
##########################################
# #
# Basis #
# #
##########################################
with h5py.File(qph5path,'r') as qph5:
ezfio.set_ao_basis_ao_basis(qph5['ao_basis'].attrs['ao_basis'])
ezfio.set_ao_basis_ao_nucl(qph5['ao_basis/ao_nucl'][()].tolist())
#Just need one (can clean this up later)
ao_prim_num_max = 5
d = [ [0] *ao_prim_num_max]*ao_num
ezfio.set_ao_basis_ao_prim_num([ao_prim_num_max]*ao_num)
ezfio.set_ao_basis_ao_power(d)
ezfio.set_ao_basis_ao_coef(d)
ezfio.set_ao_basis_ao_expo(d)
##########################################
# #
# MO Coef #
# #
##########################################
with h5py.File(qph5path,'r') as qph5:
mo_coef_reim = qph5['mo_basis/mo_coef_complex'][()].tolist()
ezfio.set_mo_basis_mo_coef_complex(mo_coef_reim)
#maybe fix qp so we don't need this?
#ezfio.set_mo_basis_mo_coef([[i for i in range(mo_num)] * ao_num])
##########################################
# #
# Integrals Mono #
# #
##########################################
with h5py.File(qph5path,'r') as qph5:
if 'ao_one_e_ints' in qph5.keys():
kin_ao_reim=qph5['ao_one_e_ints/ao_integrals_kinetic'][()].tolist()
ovlp_ao_reim=qph5['ao_one_e_ints/ao_integrals_overlap'][()].tolist()
ne_ao_reim=qph5['ao_one_e_ints/ao_integrals_n_e'][()].tolist()
ezfio.set_ao_one_e_ints_ao_integrals_kinetic_complex(kin_ao_reim)
ezfio.set_ao_one_e_ints_ao_integrals_overlap_complex(ovlp_ao_reim)
ezfio.set_ao_one_e_ints_ao_integrals_n_e_complex(ne_ao_reim)
ezfio.set_ao_one_e_ints_io_ao_integrals_kinetic('Read')
ezfio.set_ao_one_e_ints_io_ao_integrals_overlap('Read')
ezfio.set_ao_one_e_ints_io_ao_integrals_n_e('Read')
with h5py.File(qph5path,'r') as qph5:
if 'mo_one_e_ints' in qph5.keys():
kin_mo_reim=qph5['mo_one_e_ints/mo_integrals_kinetic'][()].tolist()
#ovlp_mo_reim=qph5['mo_one_e_ints/mo_integrals_overlap'][()].tolist()
ne_mo_reim=qph5['mo_one_e_ints/mo_integrals_n_e'][()].tolist()
ezfio.set_mo_one_e_ints_mo_integrals_kinetic_complex(kin_mo_reim)
#ezfio.set_mo_one_e_ints_mo_integrals_overlap_complex(ovlp_mo_reim)
#ezfio.set_mo_one_e_ints_mo_integrals_n_e_complex(ne_mo_reim)
ezfio.set_mo_one_e_ints_mo_integrals_e_n_complex(ne_mo_reim)
ezfio.set_mo_one_e_ints_io_mo_integrals_kinetic('Read')
#ezfio.set_mo_one_e_ints_io_mo_integrals_overlap('Read')
#ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read')
ezfio.set_mo_one_e_ints_io_mo_integrals_e_n('Read')
##########################################
# #
# k-points #
# #
##########################################
with h5py.File(qph5path,'r') as qph5:
kconserv = qph5['nuclei/kconserv'][()].tolist()
ezfio.set_nuclei_kconserv(kconserv)
ezfio.set_nuclei_io_kconserv('Read')
##########################################
# #
# Integrals Bi #
# #
##########################################
# should this be in ao_basis? ao_two_e_ints?
with h5py.File(qph5path,'r') as qph5:
if 'ao_two_e_ints' in qph5.keys():
df_num = qph5['ao_two_e_ints'].attrs['df_num']
ezfio.set_ao_two_e_ints_df_num(df_num)
if 'df_ao_integrals' in qph5['ao_two_e_ints'].keys():
# dfao_re0=qph5['ao_two_e_ints/df_ao_integrals_real'][()].transpose((3,2,1,0))
# dfao_im0=qph5['ao_two_e_ints/df_ao_integrals_imag'][()].transpose((3,2,1,0))
# dfao_cmplx0 = np.stack((dfao_re0,dfao_im0),axis=-1).tolist()
# ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_cmplx0)
dfao_reim=qph5['ao_two_e_ints/df_ao_integrals'][()].tolist()
ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_reim)
ezfio.set_ao_two_e_ints_io_df_ao_integrals('Read')
if 'mo_two_e_ints' in qph5.keys():
df_num = qph5['ao_two_e_ints'].attrs['df_num']
ezfio.set_ao_two_e_ints_df_num(df_num)
# dfmo_re0=qph5['mo_two_e_ints/df_mo_integrals_real'][()].transpose((3,2,1,0))
# dfmo_im0=qph5['mo_two_e_ints/df_mo_integrals_imag'][()].transpose((3,2,1,0))
# dfmo_cmplx0 = np.stack((dfmo_re0,dfmo_im0),axis=-1).tolist()
# ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_cmplx0)
dfmo_reim=qph5['mo_two_e_ints/df_mo_integrals'][()].tolist()
ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_reim)
ezfio.set_mo_two_e_ints_io_df_mo_integrals('Read')
return
#TODO: add check and only do this if ints exist #TODO: add check and only do this if ints exist
#dfmo_re=qph5['mo_two_e_ints/df_mo_integrals_real'][()].transpose((3,2,1,0)).tolist() #dfmo_re=qph5['mo_two_e_ints/df_mo_integrals_real'][()].transpose((3,2,1,0)).tolist()
#dfmo_im=qph5['mo_two_e_ints/df_mo_integrals_imag'][()].transpose((3,2,1,0)).tolist() #dfmo_im=qph5['mo_two_e_ints/df_mo_integrals_imag'][()].transpose((3,2,1,0)).tolist()
#ezfio.set_mo_two_e_ints_df_mo_integrals_real(dfmo_re) #ezfio.set_mo_two_e_ints_df_mo_integrals_real(dfmo_re)
#ezfio.set_mo_two_e_ints_df_mo_integrals_imag(dfmo_im) #ezfio.set_mo_two_e_ints_df_mo_integrals_imag(dfmo_im)
convert_kpts(fname,qph5name)