From 8794296f37a4704d99c34fe0d7f8ecadb3520aee Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Mon, 17 Feb 2020 16:16:46 -0600 Subject: [PATCH] updated converters and fixed ao df ints --- src/ao_one_e_ints/ao_overlap.irp.f | 2 +- src/ao_two_e_ints/df_ao_ints.irp.f | 19 ++-- src/hartree_fock/hf_energy.irp.f | 32 +++---- src/hartree_fock/print_e_scf.irp.f | 19 ++++ src/hartree_fock/scf.irp.f | 4 + src/scf_utils/print_debug_scf_complex.irp.f | 24 +++-- .../scf_density_matrix_ao_complex.irp.f | 18 ++-- .../Gen_Ezfio_from_integral_complex_3idx.sh | 3 + src/utils_complex/MolPyscfToQPkpts.py | 17 +++- src/utils_complex/import_ao_2e_complex.irp.f | 89 +++++++++++++++++++ src/utils_complex/import_kconserv.irp.f | 2 +- 11 files changed, 187 insertions(+), 42 deletions(-) create mode 100644 src/hartree_fock/print_e_scf.irp.f create mode 100644 src/utils_complex/import_ao_2e_complex.irp.f diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f index 5afadbbe..ad9fcff5 100644 --- a/src/ao_one_e_ints/ao_overlap.irp.f +++ b/src/ao_one_e_ints/ao_overlap.irp.f @@ -260,7 +260,7 @@ BEGIN_PROVIDER [ complex*16, S_half_inv_complex, (AO_num,AO_num) ] integer :: info, i, j, k double precision, parameter :: threshold_overlap_AO_eigenvalues = 1.d-6 - LDA = size(AO_overlap,1) + LDA = size(AO_overlap_complex,1) LDC = size(S_half_inv_complex,1) allocate( & diff --git a/src/ao_two_e_ints/df_ao_ints.irp.f b/src/ao_two_e_ints/df_ao_ints.irp.f index 6cd508ca..2200fbb5 100644 --- a/src/ao_two_e_ints/df_ao_ints.irp.f +++ b/src/ao_two_e_ints/df_ao_ints.irp.f @@ -116,8 +116,18 @@ subroutine ao_map_fill_from_df do kl=1, kpt_num do kj=1, kl call idx2_tri_int(kj,kl,kjkl2) - print*,'kj,kl,kjkl2',kj,kl,kjkl2 - ints_jl = df_ao_integrals_complex(:,:,:,kjkl2) + !print*,'kj,kl,kjkl2',kj,kl,kjkl2 + if (kj < kl) then + do i_ao=1,ao_num_per_kpt + do j_ao=1,ao_num_per_kpt + do i_df=1,df_num + ints_jl(i_ao,j_ao,i_df) = dconjg(df_ao_integrals_complex(j_ao,i_ao,i_df,kjkl2)) + enddo + enddo + enddo + else + ints_jl = df_ao_integrals_complex(:,:,:,kjkl2) + endif !$OMP PARALLEL PRIVATE(i,k,j,l,ki,kk,ii,ik,ij,il,kikk2,jl2,ik2, & !$OMP ints_ik, ints_ikjl, i_ao, j_ao, i_df, & @@ -143,15 +153,12 @@ subroutine ao_map_fill_from_df ki=kconserv(kl,kk,kj) if ((kl == kj) .and. (ki > kk)) cycle call idx2_tri_int(ki,kk,kikk2) - print*,'ki,kk,kikk2',ki,kk,kikk2 if (kikk2 > kjkl2) cycle - if (ki >= kk) then - !if (ki < kk) then !this didn't fix the problem + if (ki < kk) then do i_ao=1,ao_num_per_kpt do j_ao=1,ao_num_per_kpt do i_df=1,df_num ints_ik(i_ao,j_ao,i_df) = dconjg(df_ao_integrals_complex(j_ao,i_ao,i_df,kikk2)) - !ints_ik(j_ao,i_ao,i_df) = dconjg(df_ao_integrals_complex(j_ao,i_ao,i_df,kikk2)) enddo enddo enddo diff --git a/src/hartree_fock/hf_energy.irp.f b/src/hartree_fock/hf_energy.irp.f index db723600..9a5e6d1d 100644 --- a/src/hartree_fock/hf_energy.irp.f +++ b/src/hartree_fock/hf_energy.irp.f @@ -11,50 +11,50 @@ BEGIN_PROVIDER [double precision, extra_e_contrib_density] END_PROVIDER - BEGIN_PROVIDER [ double precision, HF_energy] -&BEGIN_PROVIDER [ double precision, HF_two_electron_energy] -&BEGIN_PROVIDER [ double precision, HF_one_electron_energy] + BEGIN_PROVIDER [ double precision, hf_energy] +&BEGIN_PROVIDER [ double precision, hf_two_electron_energy] +&BEGIN_PROVIDER [ double precision, hf_one_electron_energy] implicit none BEGIN_DOC ! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components. END_DOC integer :: i,j - HF_energy = nuclear_repulsion - HF_two_electron_energy = 0.d0 - HF_one_electron_energy = 0.d0 + hf_energy = nuclear_repulsion + hf_two_electron_energy = 0.d0 + hf_one_electron_energy = 0.d0 if (is_complex) then 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) ) + 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) ) enddo enddo if (dabs(dimag(hf_2e_tmp)).gt.1.d-10) then print*,'HF_2e energy should be real:',irp_here stop -1 else - HF_two_electron_energy = dble(hf_2e_tmp) + hf_two_electron_energy = dble(hf_2e_tmp) endif if (dabs(dimag(hf_1e_tmp)).gt.1.d-10) then print*,'HF_1e energy should be real:',irp_here stop -1 else - HF_one_electron_energy = dble(hf_1e_tmp) + hf_one_electron_energy = dble(hf_1e_tmp) endif else do j=1,ao_num do i=1,ao_num - HF_two_electron_energy += 0.5d0 * ( ao_two_e_integral_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) & - +ao_two_e_integral_beta(i,j) * SCF_density_matrix_ao_beta(i,j) ) - HF_one_electron_energy += ao_one_e_integrals(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) ) + hf_two_electron_energy += 0.5d0 * ( ao_two_e_integral_alpha(i,j) * scf_density_matrix_ao_alpha(i,j) & + +ao_two_e_integral_beta(i,j) * scf_density_matrix_ao_beta(i,j) ) + hf_one_electron_energy += ao_one_e_integrals(i,j) * (scf_density_matrix_ao_alpha(i,j) + scf_density_matrix_ao_beta (i,j) ) enddo enddo endif - HF_energy += HF_two_electron_energy + HF_one_electron_energy + hf_energy += hf_two_electron_energy + hf_one_electron_energy END_PROVIDER diff --git a/src/hartree_fock/print_e_scf.irp.f b/src/hartree_fock/print_e_scf.irp.f new file mode 100644 index 00000000..65e97a56 --- /dev/null +++ b/src/hartree_fock/print_e_scf.irp.f @@ -0,0 +1,19 @@ +program print_e_scf + call run +end + +subroutine run + + use bitmasks + implicit none + + call print_debug_scf_complex + + print*,'hf 1e,2e,total energy' + print*,hf_one_electron_energy + print*,hf_two_electron_energy + print*,hf_energy + +end + + diff --git a/src/hartree_fock/scf.irp.f b/src/hartree_fock/scf.irp.f index 8dddda92..276b4e65 100644 --- a/src/hartree_fock/scf.irp.f +++ b/src/hartree_fock/scf.irp.f @@ -98,6 +98,10 @@ subroutine run call roothaan_hall_scf endif call ezfio_set_hartree_fock_energy(SCF_energy) + print*,'hf 1e,2e,total energy' + print*,hf_one_electron_energy + print*,hf_two_electron_energy + print*,hf_energy end diff --git a/src/scf_utils/print_debug_scf_complex.irp.f b/src/scf_utils/print_debug_scf_complex.irp.f index 91311c58..65a047c3 100644 --- a/src/scf_utils/print_debug_scf_complex.irp.f +++ b/src/scf_utils/print_debug_scf_complex.irp.f @@ -15,26 +15,36 @@ subroutine print_debug_scf_complex do i=1,ao_num write(*,'(200(E24.15))') scf_density_matrix_ao_alpha_complex(i,:) enddo - write(*,'(A)') 'scf_density_matrix_ao_beta_complex' + write(*,'(A)') 'ao_one_e_integrals_complex' write(*,'(A)') '---------------' do i=1,ao_num - write(*,'(200(E24.15))') scf_density_matrix_ao_beta_complex(i,:) + write(*,'(200(E24.15))') ao_one_e_integrals_complex(i,:) enddo write(*,'(A)') 'ao_two_e_integral_alpha_complex' write(*,'(A)') '---------------' do i=1,ao_num write(*,'(200(E24.15))') ao_two_e_integral_alpha_complex(i,:) enddo - write(*,'(A)') 'ao_two_e_integral_beta_complex' - write(*,'(A)') '---------------' - do i=1,ao_num - write(*,'(200(E24.15))') ao_two_e_integral_beta_complex(i,:) - enddo write(*,'(A)') 'fock_matrix_ao_alpha_complex' write(*,'(A)') '---------------' do i=1,ao_num write(*,'(200(E24.15))') fock_matrix_ao_alpha_complex(i,:) enddo + write(*,'(A)') 'ao_overlap_complex' + write(*,'(A)') '---------------' + do i=1,ao_num + write(*,'(200(E24.15))') ao_overlap_complex(i,:) + enddo + write(*,'(A)') 'scf_density_matrix_ao_beta_complex' + write(*,'(A)') '---------------' + do i=1,ao_num + write(*,'(200(E24.15))') scf_density_matrix_ao_beta_complex(i,:) + enddo + write(*,'(A)') 'ao_two_e_integral_beta_complex' + write(*,'(A)') '---------------' + do i=1,ao_num + write(*,'(200(E24.15))') ao_two_e_integral_beta_complex(i,:) + enddo write(*,'(A)') 'fock_matrix_ao_beta_complex' write(*,'(A)') '---------------' do i=1,ao_num diff --git a/src/scf_utils/scf_density_matrix_ao_complex.irp.f b/src/scf_utils/scf_density_matrix_ao_complex.irp.f index 2bf7b77e..6e22e209 100644 --- a/src/scf_utils/scf_density_matrix_ao_complex.irp.f +++ b/src/scf_utils/scf_density_matrix_ao_complex.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [complex*16, SCF_density_matrix_ao_alpha_complex, (ao_num,ao_num) ] +BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_alpha_complex, (ao_num,ao_num) ] implicit none BEGIN_DOC ! $C.C^t$ over $\alpha$ MOs @@ -7,11 +7,11 @@ BEGIN_PROVIDER [complex*16, SCF_density_matrix_ao_alpha_complex, (ao_num,ao_num) call zgemm('N','C',ao_num,ao_num,elec_alpha_num,(1.d0,0.d0), & mo_coef_complex, size(mo_coef_complex,1), & mo_coef_complex, size(mo_coef_complex,1), (0.d0,0.d0), & - SCF_density_matrix_ao_alpha_complex, size(SCF_density_matrix_ao_alpha_complex,1)) + scf_density_matrix_ao_alpha_complex, size(scf_density_matrix_ao_alpha_complex,1)) END_PROVIDER -BEGIN_PROVIDER [ complex*16, SCF_density_matrix_ao_beta_complex, (ao_num,ao_num) ] +BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_beta_complex, (ao_num,ao_num) ] implicit none BEGIN_DOC ! $C.C^t$ over $\beta$ MOs @@ -20,21 +20,21 @@ BEGIN_PROVIDER [ complex*16, SCF_density_matrix_ao_beta_complex, (ao_num,ao_num call zgemm('N','C',ao_num,ao_num,elec_beta_num,(1.d0,0.d0), & mo_coef_complex, size(mo_coef_complex,1), & mo_coef_complex, size(mo_coef_complex,1), (0.d0,0.d0), & - SCF_density_matrix_ao_beta_complex, size(SCF_density_matrix_ao_beta_complex,1)) + scf_density_matrix_ao_beta_complex, size(scf_density_matrix_ao_beta_complex,1)) END_PROVIDER -BEGIN_PROVIDER [ complex*16, SCF_density_matrix_ao_complex, (ao_num,ao_num) ] +BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_complex, (ao_num,ao_num) ] implicit none BEGIN_DOC ! Sum of $\alpha$ and $\beta$ density matrices END_DOC - ASSERT (size(SCF_density_matrix_ao_complex,1) == size(SCF_density_matrix_ao_alpha_complex,1)) + ASSERT (size(scf_density_matrix_ao_complex,1) == size(scf_density_matrix_ao_alpha_complex,1)) if (elec_alpha_num== elec_beta_num) then - SCF_density_matrix_ao_complex = SCF_density_matrix_ao_alpha_complex + SCF_density_matrix_ao_alpha_complex + scf_density_matrix_ao_complex = scf_density_matrix_ao_alpha_complex + scf_density_matrix_ao_alpha_complex else - ASSERT (size(SCF_density_matrix_ao_complex,1) == size(SCF_density_matrix_ao_beta_complex ,1)) - SCF_density_matrix_ao_complex = SCF_density_matrix_ao_alpha_complex + SCF_density_matrix_ao_beta_complex + ASSERT (size(scf_density_matrix_ao_complex,1) == size(scf_density_matrix_ao_beta_complex ,1)) + scf_density_matrix_ao_complex = scf_density_matrix_ao_alpha_complex + scf_density_matrix_ao_beta_complex endif END_PROVIDER diff --git a/src/utils_complex/Gen_Ezfio_from_integral_complex_3idx.sh b/src/utils_complex/Gen_Ezfio_from_integral_complex_3idx.sh index 31d273c1..9bce8816 100755 --- a/src/utils_complex/Gen_Ezfio_from_integral_complex_3idx.sh +++ b/src/utils_complex/Gen_Ezfio_from_integral_complex_3idx.sh @@ -17,6 +17,9 @@ echo 'Create EZFIO' qp_edit -c $ezfio &> /dev/null #cp $ezfio/{ao,mo}_basis/ao_md5 +qp_run import_kconserv $ezfio +#qp_run import_ao_2e_complex $ezfio +#qp_run dump_ao_2e_from_df $ezfio #Read the integral #echo 'Read Integral' diff --git a/src/utils_complex/MolPyscfToQPkpts.py b/src/utils_complex/MolPyscfToQPkpts.py index d909c99e..e0284096 100644 --- a/src/utils_complex/MolPyscfToQPkpts.py +++ b/src/utils_complex/MolPyscfToQPkpts.py @@ -508,6 +508,7 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, natom = cell.natm nelec = cell.nelectron + neleca,nelecb = cell.nelec atom_xyz = mf.cell.atom_coords() if not(mf.cell.unit.startswith(('B','b','au','AU'))): atom_xyz /= nist.BOHR # always convert to au @@ -537,8 +538,8 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, print("n df fitting functions", naux) #in old version: param << nelec*Nk, nmo*Nk, natom*Nk - qph5['electrons'].attrs['elec_alpha_num']=nelec*Nk - qph5['electrons'].attrs['elec_beta_num']=nelec*Nk + qph5['electrons'].attrs['elec_alpha_num']=neleca*Nk + qph5['electrons'].attrs['elec_beta_num']=nelecb*Nk qph5['mo_basis'].attrs['mo_num']=Nk*nmo qph5['ao_basis'].attrs['ao_num']=Nk*nao qph5['nuclei'].attrs['nucl_num']=Nk*natom @@ -572,6 +573,18 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, qph5.create_dataset('mo_basis/mo_coef_imag',data=mo_coef_blocked.imag) 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) + + with open('C.qp','w') as outfile: + c_kpts = np.reshape(mo_k,(Nk,nao,nmo)) + + for ik in range(Nk): + shift1=ik*nao+1 + shift2=ik*nmo+1 + for i in range(nao): + for j in range(nmo): + cij = c_kpts[ik,i,j] + if abs(cij) > mo_coef_threshold: + outfile.write('%s %s %s %s\n' % (i+shift1, j+shift2, cij.real, cij.imag)) # ___ # | ._ _|_ _ _ ._ _. | _ |\/| _ ._ _ diff --git a/src/utils_complex/import_ao_2e_complex.irp.f b/src/utils_complex/import_ao_2e_complex.irp.f new file mode 100644 index 00000000..f22e6c50 --- /dev/null +++ b/src/utils_complex/import_ao_2e_complex.irp.f @@ -0,0 +1,89 @@ +program import_ao_2e_complex + call run +end + +subroutine run + use map_module + implicit none + + integer :: iunit + integer :: getunitandopen + + integer ::i,j,k,l + double precision :: integral + complex*16, allocatable :: C(:,:) + double precision :: tmp_re, tmp_im + + integer :: n_integrals_1, n_integrals_2 + integer(key_kind), allocatable :: buffer_i_1(:), buffer_i_2(:) + real(integral_kind), allocatable :: buffer_values_1(:), buffer_values_2(:) + logical :: use_map1 + integer(key_kind) :: idx_tmp + double precision :: sign + + +! call ezfio_set_ao_basis_ao_num(ao_num) + + allocate(buffer_i_1(ao_num**3), buffer_values_1(ao_num**3)) + allocate(buffer_i_2(ao_num**3), buffer_values_2(ao_num**3)) + iunit = getunitandopen('W.qp','r') + n_integrals_1=0 + n_integrals_2=0 + buffer_values_1 = 0.d0 + buffer_values_2 = 0.d0 + do + read (iunit,*,end=13) i,j,k,l, tmp_re, tmp_im + call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign) +! print'(4(I4),(L3),(I6),(F7.1))',i,j,k,l,use_map1,idx_tmp,sign + if (use_map1) then + n_integrals_1 += 1 + buffer_i_1(n_integrals_1)=idx_tmp + buffer_values_1(n_integrals_1)=tmp_re +! print'(A,4(I4),(I6),(E15.7))','map1',i,j,k,l,idx_tmp,tmp_re + if (sign.ne.0.d0) then + n_integrals_1 += 1 + buffer_i_1(n_integrals_1)=idx_tmp+1 + buffer_values_1(n_integrals_1)=tmp_im*sign +! print'(A,4(I4),(I6),(E15.7))','map1',i,j,k,l,idx_tmp+1,tmp_im*sign + endif + if (n_integrals_1 >= size(buffer_i_1)-1) then + call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + n_integrals_1 = 0 + endif + else + n_integrals_2 += 1 + buffer_i_2(n_integrals_2)=idx_tmp + buffer_values_2(n_integrals_2)=tmp_re +! print'(A,4(I4),(I6),(E15.7))','map2',i,j,k,l,idx_tmp,tmp_re + if (sign.ne.0.d0) then + n_integrals_2 += 1 + buffer_i_2(n_integrals_2)=idx_tmp+1 + buffer_values_2(n_integrals_2)=tmp_im*sign +! print'(A,4(I4),(I6),(E15.7))','map2',i,j,k,l,idx_tmp+1,tmp_im*sign + endif + if (n_integrals_2 >= size(buffer_i_2)-1) then + call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + n_integrals_2 = 0 + endif + endif + enddo + 13 continue + close(iunit) + + if (n_integrals_1 > 0) then + call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + endif + if (n_integrals_2 > 0) then + call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + endif + + call map_sort(ao_integrals_map) + call map_unique(ao_integrals_map) + call map_sort(ao_integrals_map_2) + call map_unique(ao_integrals_map_2) + + call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_1',ao_integrals_map) + call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_2',ao_integrals_map_2) + call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + +end diff --git a/src/utils_complex/import_kconserv.irp.f b/src/utils_complex/import_kconserv.irp.f index 0dff268b..d15b1eed 100644 --- a/src/utils_complex/import_kconserv.irp.f +++ b/src/utils_complex/import_kconserv.irp.f @@ -22,7 +22,7 @@ subroutine run allocate(A(kpt_num,kpt_num,kpt_num)) A = 0 - iunit = getunitandopen('kconserv','r') + iunit = getunitandopen('K.qp','r') do read (iunit,*,end=10) i,j,k,l A(i,j,k) = l