Merge pull request #122 from kgasperich/dev-real-kpts

bin/qp_convert_h5_to_ezfio

@@ -76,6 +76,8 @@ def convert_mol(filename,qph5path):
     
     ezfio.set_nuclei_nucl_charge(nucl_charge)
     ezfio.set_nuclei_nucl_coord(nucl_coord)
+    if isinstance(nucl_label[0],bytes):
+        nucl_label = list(map(lambda x:x.decode(),nucl_label))
     ezfio.set_nuclei_nucl_label(nucl_label)
     
     ezfio.set_nuclei_io_nuclear_repulsion('Read')
@@ -135,6 +137,387 @@ def convert_mol(filename,qph5path):
     
     return
 
+def convert_kpts(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.set_ao_basis_ao_num_per_kpt(ao_num//kpt_num)
+    ezfio.set_mo_basis_mo_num_per_kpt(mo_num//kpt_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)
+    if isinstance(nucl_label[0],bytes):
+        nucl_label = list(map(lambda x:x.decode(),nucl_label))
+    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_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              #
+    #                                        #
+    ##########################################
+        
+    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_kpts'][()].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_n_e_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_n_e('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)
+    if isinstance(nucl_label[0],bytes):
+        nucl_label = list(map(lambda x:x.decode(),nucl_label))
+    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_n_e_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_n_e('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
+    
+
 if __name__ == '__main__':
     ARGUMENTS = docopt(__doc__)
 
@@ -145,15 +528,21 @@ if __name__ == '__main__':
     else:
         EZFIO_FILE = "{0}.ezfio".format(FILE)
 
+    with h5py.File(FILE,'r') as qph5:
+        do_kpts = ('kconserv' in qph5['nuclei'].keys())
+    if (do_kpts):
+        print("converting HDF5 to EZFIO for periodic system")
+        convert_kpts(EZFIO_FILE,FILE)
+    else:
+        print("converting HDF5 to EZFIO for molecular system")
+        convert_mol(EZFIO_FILE,FILE)
 
-    convert_mol(EZFIO_FILE,FILE)
+#    sys.stdout.flush()
-
+#    if os.system("qp_run save_ortho_mos "+EZFIO_FILE) != 0:
-    sys.stdout.flush()
+#        print("""Warning: You need to run
-    if os.system("qp_run save_ortho_mos "+EZFIO_FILE) != 0:
+#
-        print("""Warning: You need to run
+#         qp run save_ortho_mos
-
+#
-         qp run save_ortho_mos
+#to be sure your MOs will be orthogonal, which is not the case when
-
+#the MOs are read from output files (not enough precision in output).""")
-to be sure your MOs will be orthogonal, which is not the case when
-the MOs are read from output files (not enough precision in output).""")
 

src/davidson/diagonalize_ci.irp.f

@@ -313,7 +313,7 @@ END_PROVIDER
          H_prime(j,j) = H_prime(j,j) + alpha*(s_z2_sz - expected_s2)
        enddo
        call lapack_diag_complex(eigenvalues,eigenvectors,H_prime,size(H_prime,1),N_det)
-       ci_electronic_energy_complex(:) = (0.d0,0.d0)
+       ci_electronic_energy_complex(:) = 0.d0
        i_state = 0
        allocate (s2_eigvalues(N_det))
        allocate(index_good_state_array(N_det),good_state_array(N_det))

src/hartree_fock/print_e_scf.irp.f

@@ -7,7 +7,9 @@ subroutine run
   use bitmasks
   implicit none
 
-  call print_debug_scf_complex
+  !if (is_complex) then
+  !  call print_debug_scf_complex
+  !endif
 
   print*,'hf 1e,2e,total energy'
   print*,hf_one_electron_energy

src/mo_guess/pot_mo_ortho_lowdin_ints.irp.f

@@ -28,7 +28,7 @@ BEGIN_PROVIDER [complex*16, ao_ortho_lowdin_n_e_ints_cplx, (mo_num,mo_num)]
  integer :: i1,j1,i,j
  complex*16 :: c_i1,c_j1
 
- ao_ortho_lowdin_nucl_elec_integrals = (0.d0,0.d0)
+ ao_ortho_lowdin_n_e_ints_cplx = (0.d0,0.d0)
  !$OMP PARALLEL DO DEFAULT(none) &
  !$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) &
  !$OMP SHARED(mo_num,ao_num,ao_ortho_lowdin_coef_complex, &