From d0fe9aad4f97fb033602ab28fe4b312d6bb50b73 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Fri, 20 Mar 2020 12:22:10 -0500 Subject: [PATCH] scf kpts --- src/ao_one_e_ints/ao_overlap.irp.f | 108 ++++ src/bitmask/track_orb.irp.f | 104 +++- src/hartree_fock/hf_energy.irp.f | 16 +- src/hartree_fock/scf.irp.f | 3 +- src/mo_basis/utils_cplx.irp.f | 20 +- src/mo_one_e_ints/EZFIO.cfg | 12 + src/mo_one_e_ints/mo_overlap.irp.f | 54 ++ src/mo_one_e_ints/orthonormalize.irp.f | 12 +- src/scf_utils/diagonalize_fock_cplx.irp.f | 4 +- src/scf_utils/diis_cplx.irp.f | 152 +++++ src/scf_utils/fock_matrix.irp.f | 14 +- src/scf_utils/fock_matrix_cplx.irp.f | 80 +-- src/scf_utils/roothaan_hall_scf_cplx.irp.f | 334 +++++++++++ src/utils_complex/MolPyscfToQPkpts.py | 10 +- .../create_ezfio_complex_3idx.py | 553 ++++++++++++------ 15 files changed, 1193 insertions(+), 283 deletions(-) diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f index 2e1695a7..9877f882 100644 --- a/src/ao_one_e_ints/ao_overlap.irp.f +++ b/src/ao_one_e_ints/ao_overlap.irp.f @@ -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)) 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_DOC @@ -326,6 +338,66 @@ BEGIN_PROVIDER [ complex*16, S_half_inv_complex, (AO_num,AO_num) ] 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) ] implicit none @@ -395,3 +467,39 @@ BEGIN_PROVIDER [ complex*16, S_half_complex, (ao_num,ao_num) ] 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 + diff --git a/src/bitmask/track_orb.irp.f b/src/bitmask/track_orb.irp.f index 73bf78f3..9e96cca5 100644 --- a/src/bitmask/track_orb.irp.f +++ b/src/bitmask/track_orb.irp.f @@ -16,6 +16,15 @@ BEGIN_PROVIDER [ complex*16, mo_coef_begin_iteration_complex, (ao_num,mo_num) ] END_DOC 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 implicit none 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` END_DOC 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 mo_coef_begin_iteration = mo_coef endif @@ -42,37 +52,71 @@ subroutine reorder_core_orb integer :: i1,i2 if (is_complex) then complex*16, allocatable :: accu_c(:) - allocate(accu(mo_num),accu_c(mo_num),index_core_orb(n_core_orb),iorder(mo_num)) - do i = 1, n_core_orb - iorb = list_core(i) - do j = 1, mo_num - accu(j) = 0.d0 - accu_c(j) = (0.d0,0.d0) - iorder(j) = j - do k = 1, ao_num - do l = 1, ao_num - accu_c(j) += dconjg(mo_coef_begin_iteration_complex(k,iorb)) * & - mo_coef_complex(l,j) * ao_overlap_complex(k,l) - enddo - enddo - accu(j) = -cdabs(accu_c(j)) - enddo - call dsort(accu,iorder,mo_num) - index_core_orb(i) = iorder(1) - enddo + !allocate(accu(mo_num),accu_c(mo_num),index_core_orb(n_core_orb),iorder(mo_num)) + !do i = 1, n_core_orb + ! iorb = list_core(i) + ! do j = 1, mo_num + ! accu(j) = 0.d0 + ! accu_c(j) = (0.d0,0.d0) + ! iorder(j) = j + ! do k = 1, ao_num + ! do l = 1, ao_num + ! accu_c(j) += dconjg(mo_coef_begin_iteration_complex(k,iorb)) * & + ! mo_coef_complex(l,j) * ao_overlap_complex(k,l) + ! enddo + ! enddo + ! accu(j) = -cdabs(accu_c(j)) + ! enddo + ! call dsort(accu,iorder,mo_num) + ! index_core_orb(i) = iorder(1) + !enddo - complex*16 :: x_c - do j = 1, n_core_orb - i1 = list_core(j) - i2 = index_core_orb(j) - do i=1,ao_num - x_c = mo_coef_complex(i,i1) - mo_coef_complex(i,i1) = mo_coef_complex(i,i2) - mo_coef_complex(i,i2) = x_c - enddo - enddo - !call loc_cele_routine + !complex*16 :: x_c + !do j = 1, n_core_orb + ! i1 = list_core(j) + ! i2 = index_core_orb(j) + ! do i=1,ao_num + ! x_c = mo_coef_complex(i,i1) + ! mo_coef_complex(i,i1) = mo_coef_complex(i,i2) + ! mo_coef_complex(i,i2) = x_c + ! enddo + !enddo + !!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) else allocate(accu(mo_num),index_core_orb(n_core_orb),iorder(mo_num)) diff --git a/src/hartree_fock/hf_energy.irp.f b/src/hartree_fock/hf_energy.irp.f index 9a5e6d1d..5a68164f 100644 --- a/src/hartree_fock/hf_energy.irp.f +++ b/src/hartree_fock/hf_energy.irp.f @@ -18,7 +18,7 @@ END_PROVIDER BEGIN_DOC ! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components. END_DOC - integer :: i,j + integer :: i,j,k hf_energy = nuclear_repulsion hf_two_electron_energy = 0.d0 hf_one_electron_energy = 0.d0 @@ -26,12 +26,14 @@ END_PROVIDER complex*16 :: hf_1e_tmp, hf_2e_tmp hf_1e_tmp = (0.d0,0.d0) hf_2e_tmp = (0.d0,0.d0) - do j=1,ao_num - do i=1,ao_num - hf_2e_tmp += 0.5d0 * ( ao_two_e_integral_alpha_complex(i,j) * scf_density_matrix_ao_alpha_complex(j,i) & - +ao_two_e_integral_beta_complex(i,j) * scf_density_matrix_ao_beta_complex(j,i) ) - hf_1e_tmp += ao_one_e_integrals_complex(i,j) * (scf_density_matrix_ao_alpha_complex(j,i) & - + scf_density_matrix_ao_beta_complex (j,i) ) + do k=1,kpt_num + do j=1,ao_num_per_kpt + do i=1,ao_num_per_kpt + hf_2e_tmp += 0.5d0 * ( ao_two_e_integral_alpha_kpts(i,j,k) * scf_density_matrix_ao_alpha_kpts(j,i,k) & + +ao_two_e_integral_beta_kpts(i,j,k) * scf_density_matrix_ao_beta_kpts(j,i,k) ) + 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 if (dabs(dimag(hf_2e_tmp)).gt.1.d-10) then diff --git a/src/hartree_fock/scf.irp.f b/src/hartree_fock/scf.irp.f index 9d1671ec..8e438613 100644 --- a/src/hartree_fock/scf.irp.f +++ b/src/hartree_fock/scf.irp.f @@ -102,7 +102,8 @@ subroutine run mo_label = "Orthonormalized" if (is_complex) then - call roothaan_hall_scf_complex + !call roothaan_hall_scf_complex + call roothaan_hall_scf_kpts else call roothaan_hall_scf endif diff --git a/src/mo_basis/utils_cplx.irp.f b/src/mo_basis/utils_cplx.irp.f index 58230cdd..13327f57 100644 --- a/src/mo_basis/utils_cplx.irp.f +++ b/src/mo_basis/utils_cplx.irp.f @@ -52,11 +52,11 @@ subroutine mo_as_eigvectors_of_mo_matrix_complex(matrix,n,m,label,sign,output) enddo write (6,'(A)') '======== ================' write (6,'(A)') '' - write (6,'(A)') 'Fock Matrix' - write (6,'(A)') '-----------' - do i=1,n - write(*,'(200(E24.15))') A(i,:) - enddo + !write (6,'(A)') 'Fock Matrix' + !write (6,'(A)') '-----------' + !do i=1,n + ! write(*,'(200(E24.15))') A(i,:) + !enddo 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)) @@ -302,11 +302,11 @@ subroutine mo_as_eigvectors_of_mo_matrix_kpts(matrix,n,m,nk,label,sign,output) enddo write (6,'(A)') '======== ================' write (6,'(A)') '' - write (6,'(A)') 'Fock Matrix' - write (6,'(A)') '-----------' - do i=1,n - write(*,'(200(E24.15))') A(i,:) - enddo + !write (6,'(A)') 'Fock Matrix' + !write (6,'(A)') '-----------' + !do i=1,n + ! write(*,'(200(E24.15))') A(i,:) + !enddo endif call zgemm('N','N',ao_num_per_kpt,m,m,(1.d0,0.d0), & diff --git a/src/mo_one_e_ints/EZFIO.cfg b/src/mo_one_e_ints/EZFIO.cfg index bd60ca16..23b30aba 100644 --- a/src/mo_one_e_ints/EZFIO.cfg +++ b/src/mo_one_e_ints/EZFIO.cfg @@ -47,6 +47,18 @@ doc: Read/Write |MO| one-electron kinetic integrals from/to disk [ Write | Read interface: ezfio,provider,ocaml 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] type: double precision diff --git a/src/mo_one_e_ints/mo_overlap.irp.f b/src/mo_one_e_ints/mo_overlap.irp.f index 796c9fde..f004e1f4 100644 --- a/src/mo_one_e_ints/mo_overlap.irp.f +++ b/src/mo_one_e_ints/mo_overlap.irp.f @@ -74,3 +74,57 @@ BEGIN_PROVIDER [ complex*16, mo_overlap_complex,(mo_num,mo_num) ] 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 + diff --git a/src/mo_one_e_ints/orthonormalize.irp.f b/src/mo_one_e_ints/orthonormalize.irp.f index 11a09b4e..d7949ace 100644 --- a/src/mo_one_e_ints/orthonormalize.irp.f +++ b/src/mo_one_e_ints/orthonormalize.irp.f @@ -1,12 +1,14 @@ subroutine orthonormalize_mos implicit none - integer :: m,p,s + integer :: m,p,s,k if (is_complex) then - m = size(mo_coef_complex,1) - p = size(mo_overlap_complex,1) - call ortho_lowdin_complex(mo_overlap_complex,p,mo_num,mo_coef_complex,m,ao_num) + do k=1,kpt_num + m = size(mo_coef_kpts,1) + 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' - SOFT_TOUCH mo_coef_complex mo_label + SOFT_TOUCH mo_coef_kpts mo_label else m = size(mo_coef,1) p = size(mo_overlap,1) diff --git a/src/scf_utils/diagonalize_fock_cplx.irp.f b/src/scf_utils/diagonalize_fock_cplx.irp.f index 83d4b00f..82353ed0 100644 --- a/src/scf_utils/diagonalize_fock_cplx.irp.f +++ b/src/scf_utils/diagonalize_fock_cplx.irp.f @@ -72,8 +72,8 @@ BEGIN_PROVIDER [ complex*16, eigenvectors_Fock_matrix_mo_kpts, (ao_num_per_kpt,m allocate (diag(mo_num_per_kpt) ) do k=1,kpt_num - do j=1,mo_num - do i=1,mo_num + do j=1,mo_num_per_kpt + do i=1,mo_num_per_kpt !F(i,j) = fock_matrix_mo_complex(i,j) F(i,j) = fock_matrix_mo_kpts(i,j,k) enddo diff --git a/src/scf_utils/diis_cplx.irp.f b/src/scf_utils/diis_cplx.irp.f index 721a9751..4a0cdabf 100644 --- a/src/scf_utils/diis_cplx.irp.f +++ b/src/scf_utils/diis_cplx.irp.f @@ -140,3 +140,155 @@ END_PROVIDER deallocate(scratch) 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 diff --git a/src/scf_utils/fock_matrix.irp.f b/src/scf_utils/fock_matrix.irp.f index a77a78fc..efd64be1 100644 --- a/src/scf_utils/fock_matrix.irp.f +++ b/src/scf_utils/fock_matrix.irp.f @@ -157,15 +157,17 @@ BEGIN_PROVIDER [ double precision, SCF_energy ] END_DOC SCF_energy = nuclear_repulsion - integer :: i,j + integer :: i,j,k if (is_complex) then complex*16 :: scf_e_tmp scf_e_tmp = dcmplx(SCF_energy,0.d0) - do j=1,ao_num - do i=1,ao_num - scf_e_tmp += 0.5d0 * ( & - (ao_one_e_integrals_complex(i,j) + Fock_matrix_ao_alpha_complex(i,j) ) * SCF_density_matrix_ao_alpha_complex(j,i) +& - (ao_one_e_integrals_complex(i,j) + Fock_matrix_ao_beta_complex (i,j) ) * SCF_density_matrix_ao_beta_complex (j,i) ) + do k=1,kpt_num + do j=1,ao_num_per_kpt + do i=1,ao_num_per_kpt + scf_e_tmp += 0.5d0 * ( & + (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 !TODO: add check for imaginary part? (should be zero) diff --git a/src/scf_utils/fock_matrix_cplx.irp.f b/src/scf_utils/fock_matrix_cplx.irp.f index afa21072..cc0dc4af 100644 --- a/src/scf_utils/fock_matrix_cplx.irp.f +++ b/src/scf_utils/fock_matrix_cplx.irp.f @@ -593,14 +593,14 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) - kpt_i = (i-1)/kpt_num +1 - kpt_j = (j-1)/kpt_num +1 - kpt_k = (k-1)/kpt_num +1 - kpt_l = (l-1)/kpt_num +1 - idx_i = mod(i,kpt_num) - idx_j = mod(j,kpt_num) - idx_k = mod(k,kpt_num) - idx_l = mod(l,kpt_num) + kpt_i = (i-1)/ao_num_per_kpt +1 + kpt_j = (j-1)/ao_num_per_kpt +1 + kpt_k = (k-1)/ao_num_per_kpt +1 + kpt_l = (l-1)/ao_num_per_kpt +1 + idx_i = mod(i-1,ao_num_per_kpt)+1 + idx_j = mod(j-1,ao_num_per_kpt)+1 + idx_k = mod(k-1,ao_num_per_kpt)+1 + idx_l = mod(l-1,ao_num_per_kpt)+1 integral = i_sign(k2)*values(k1) !for klij and lkji, take complex conjugate !G_a(i,k) += D_{ab}(l,j)*() @@ -611,7 +611,7 @@ END_PROVIDER 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 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 endif 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_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 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 @@ -636,20 +636,20 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) - kpt_i = (i-1)/kpt_num +1 - kpt_j = (j-1)/kpt_num +1 - kpt_k = (k-1)/kpt_num +1 - kpt_l = (l-1)/kpt_num +1 - idx_i = mod(i,kpt_num) - idx_j = mod(j,kpt_num) - idx_k = mod(k,kpt_num) - idx_l = mod(l,kpt_num) + kpt_i = (i-1)/ao_num_per_kpt +1 + kpt_j = (j-1)/ao_num_per_kpt +1 + kpt_k = (k-1)/ao_num_per_kpt +1 + kpt_l = (l-1)/ao_num_per_kpt +1 + idx_i = mod(i-1,ao_num_per_kpt)+1 + idx_j = mod(j-1,ao_num_per_kpt)+1 + idx_k = mod(k-1,ao_num_per_kpt)+1 + idx_l = mod(l-1,ao_num_per_kpt)+1 integral = values(k1) 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 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 endif 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_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 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 @@ -714,14 +714,14 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) - kpt_i = (i-1)/kpt_num +1 - kpt_j = (j-1)/kpt_num +1 - kpt_k = (k-1)/kpt_num +1 - kpt_l = (l-1)/kpt_num +1 - idx_i = mod(i,kpt_num) - idx_j = mod(j,kpt_num) - idx_k = mod(k,kpt_num) - idx_l = mod(l,kpt_num) + kpt_i = (i-1)/ao_num_per_kpt +1 + kpt_j = (j-1)/ao_num_per_kpt +1 + kpt_k = (k-1)/ao_num_per_kpt +1 + kpt_l = (l-1)/ao_num_per_kpt +1 + idx_i = mod(i-1,ao_num_per_kpt)+1 + idx_j = mod(j-1,ao_num_per_kpt)+1 + idx_k = mod(k-1,ao_num_per_kpt)+1 + idx_l = mod(l-1,ao_num_per_kpt)+1 integral = i_sign(k2)*values(k1) ! for klij and lkji, take conjugate !G_a(i,k) += D_{ab}(l,j)*() @@ -732,7 +732,7 @@ END_PROVIDER 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 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 endif 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_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 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 @@ -757,20 +757,20 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) - kpt_i = (i-1)/kpt_num +1 - kpt_j = (j-1)/kpt_num +1 - kpt_k = (k-1)/kpt_num +1 - kpt_l = (l-1)/kpt_num +1 - idx_i = mod(i,kpt_num) - idx_j = mod(j,kpt_num) - idx_k = mod(k,kpt_num) - idx_l = mod(l,kpt_num) + kpt_i = (i-1)/ao_num_per_kpt +1 + kpt_j = (j-1)/ao_num_per_kpt +1 + kpt_k = (k-1)/ao_num_per_kpt +1 + kpt_l = (l-1)/ao_num_per_kpt +1 + idx_i = mod(i-1,ao_num_per_kpt)+1 + idx_j = mod(j-1,ao_num_per_kpt)+1 + idx_k = mod(k-1,ao_num_per_kpt)+1 + idx_l = mod(l-1,ao_num_per_kpt)+1 integral = values(k1) 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 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 endif 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_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 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 diff --git a/src/scf_utils/roothaan_hall_scf_cplx.irp.f b/src/scf_utils/roothaan_hall_scf_cplx.irp.f index 2a68a282..e074f884 100644 --- a/src/scf_utils/roothaan_hall_scf_cplx.irp.f +++ b/src/scf_utils/roothaan_hall_scf_cplx.irp.f @@ -319,3 +319,337 @@ END_DOC endif 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 diff --git a/src/utils_complex/MolPyscfToQPkpts.py b/src/utils_complex/MolPyscfToQPkpts.py index 55c41181..2f74c089 100644 --- a/src/utils_complex/MolPyscfToQPkpts.py +++ b/src/utils_complex/MolPyscfToQPkpts.py @@ -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_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_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) @@ -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_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_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'), (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_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_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_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'), (ovlp_mo, ne_mo, kin_mo)): print_kpts_unblocked_upper(ints,fname,thresh_mono) diff --git a/src/utils_complex/create_ezfio_complex_3idx.py b/src/utils_complex/create_ezfio_complex_3idx.py index 2947b565..ae6be312 100755 --- a/src/utils_complex/create_ezfio_complex_3idx.py +++ b/src/utils_complex/create_ezfio_complex_3idx.py @@ -4,198 +4,389 @@ import h5py import sys import numpy as np -filename = sys.argv[1] -qph5path = sys.argv[2] +fname = sys.argv[1] +qph5name = sys.argv[2] -ezfio.set_file(filename) #qph5=h5py.File(qph5path,'r') - -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 # -# # -########################################## +def convert_kpts(filename,qph5path): + ezfio.set_file(filename) + ezfio.set_nuclei_is_complex(True) - -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: + 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'] -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_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 # + # # + ########################################## - 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) + with h5py.File(qph5path,'r') as qph5: + mo_coef_kpts = qph5['mo_basis/mo_coef_kpts'][()].tolist() + mo_coef_cplx = qph5['mo_basis/mo_coef_complex'][()].tolist() + ezfio.set_mo_basis_mo_coef_kpts(mo_coef_kpts) + ezfio.set_mo_basis_mo_coef_complex(mo_coef_cplx) + #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 # + # # + ########################################## - 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') - + 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_kpts'][()].tolist() + 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) + 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') + 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_kpts'][()].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_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) + #ezfio.set_mo_one_e_ints_mo_integrals_n_e_complex(ne_mo_reim) + ezfio.set_mo_one_e_ints_mo_integrals_e_n_kpts(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 +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 #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() #ezfio.set_mo_two_e_ints_df_mo_integrals_real(dfmo_re) #ezfio.set_mo_two_e_ints_df_mo_integrals_imag(dfmo_im) + +convert_kpts(fname,qph5name)