diff --git a/src/libra_py/workflows/nbra/step2.py b/src/libra_py/workflows/nbra/step2.py index e42ff4b54..426a18c74 100755 --- a/src/libra_py/workflows/nbra/step2.py +++ b/src/libra_py/workflows/nbra/step2.py @@ -419,7 +419,7 @@ def run_cp2k_libint_step2(params): # Now try to get parameters from the input critical_params = [ "cp2k_ot_input_template", "cp2k_diag_input_template", "trajectory_xyz_filename" ] - default_params = {"res_dir": os.getcwd() + "/res", "all_logfiles": os.getcwd() + "/all_logfiles", "all_pdosfiles": os.getcwd() + "/all_pdosfiles", "all_images": os.getcwd() + "/all_images", "image_format": 'bmp', "istep": 0, "fstep": 2, "lowest_orbital": 1, "highest_orbital": 2, "is_spherical": True, "isXTB": False, "isUKS": False, "remove_molden": True, "nprocs": 2, "cp2k_exe": "cp2k.psmp", "mpi_executable": "mpirun", "cube_visualization": False, "vmd_input_template": "vmd.tcl", "states_to_plot": [1], "plot_phase_corrected": True, "vmd_exe": "vmd", "tachyon_exe": "tachyon_LINIXAMD64", "x_pixels": 1024, "y_pixels": 1024, "remove_cube": True, 'together_mode': False, 'restart_file': True} + default_params = {"res_dir": os.getcwd() + "/res", "all_logfiles": os.getcwd() + "/all_logfiles", "all_pdosfiles": os.getcwd() + "/all_pdosfiles", "all_images": os.getcwd() + "/all_images", "image_format": 'bmp', "istep": 0, "fstep": 2, "lowest_orbital": 1, "highest_orbital": 2, "is_spherical": True, "isXTB": False, "isUKS": False, "remove_molden": True, "nprocs": 2, "cp2k_exe": "cp2k.psmp", "mpi_executable": "mpirun", "cube_visualization": False, "vmd_input_template": "vmd.tcl", "states_to_plot": [1], "plot_phase_corrected": True, "vmd_exe": "vmd", "tachyon_exe": "tachyon_LINIXAMD64", "x_pixels": 1024, "y_pixels": 1024, "remove_cube": True, 'together_mode': False, 'restart_file': True, 'use_emultipole': False, 'emultipole_components': ['S', 'x', 'y', 'z']} comn.check_input(params, default_params, critical_params) @@ -616,6 +616,7 @@ def run_cp2k_libint_step2(params): zero_mat,zero_mat.T,\ np.diag(beta_energies_1)[lowest_orbital-1:highest_orbital,lowest_orbital-1:highest_orbital]) E_step_sparse = scipy.sparse.csc_matrix(E_step) + scipy.sparse.save_npz(params['res_dir']+F'/E_ks_{step}.npz', E_step_sparse) else: if params['use_emultipole']: diff --git a/src/libra_py/workflows/nbra/step2_many_body.py b/src/libra_py/workflows/nbra/step2_many_body.py index 19082faf0..1451d0053 100755 --- a/src/libra_py/workflows/nbra/step2_many_body.py +++ b/src/libra_py/workflows/nbra/step2_many_body.py @@ -392,7 +392,7 @@ def compute_cube_ks_overlaps( cubefiles_prev, params): -def reindex_cp2k_sd_states( ks_orbital_homo_index, ks_orbital_indicies, sd_basis_states, sd_format=2, ks_beta_homo_index=0): +def reindex_cp2k_sd_states( ks_orbital_homo_index, ks_orbital_indicies, sd_basis_states, sd_format=2, ks_beta_homo_index=0, active_space_num_occ_orbitals=0): """ ks_orbital_homo_index: Index of the homo ks orbital, from 1 ks_orbital_indicies: Range of the considered ks orbtials. Ex) [8,9,10,11], where 9 is homo orbtial index (from 1) @@ -503,6 +503,32 @@ def reindex_cp2k_sd_states( ks_orbital_homo_index, ks_orbital_indicies, sd_basis else: sd_excitation.append(sd_state) sd_basis.append( sd_excitation ) + + # Here we extend the Slater determinants based on the number + # of occupied orbitals specified by the user using + # `active_space_num_occ_orbitals` keyword + # We first compute the difference between the current number of occupied orbitals + # and the number of orbitals specified by the user + active_space_diff = int(active_space_num_occ_orbitals - len(sd_basis[0])/2) + if active_space_diff>0: + sd_basis_ = [] + # Then, for each SD we shift the orbitals so that all occupied + # orbitals specified by the user are involved in them + for i in range(len(sd_basis)): + sd_i = [] + for j in range(1, active_space_diff): + sd_i.append(j) + sd_i.append(-j) + for j in range(len(sd_basis[i])): + # For alpha spin channel + if sd_basis[i][j]>0: + sd_i.append(sd_basis[i][j]+active_space_diff) + # For beta spin channel + else: + sd_i.append(sd_basis[i][j]-active_space_diff) + sd_basis_.append(sd_i) + sd_basis = sd_basis_ + return sd_basis diff --git a/src/libra_py/workflows/nbra/step3.py b/src/libra_py/workflows/nbra/step3.py index 75d85e553..17b7ab375 100755 --- a/src/libra_py/workflows/nbra/step3.py +++ b/src/libra_py/workflows/nbra/step3.py @@ -2258,7 +2258,7 @@ def run_step3_sd_nacs_libint(params): 'apply_orthonormalization': False, 'do_state_reordering': 0, 'state_reordering_alpha': 0, 'is_many_body': False, 'num_occ_states': 1, 'num_occ_alpha': 1, 'num_unocc_alpha': 1, - 'num_occ_beta': 1, 'num_unocc_beta': 1, + 'num_occ_beta': 1, 'num_unocc_beta': 1, 'active_space_num_occ_orbitals': 0, 'num_unocc_states': 1, 'verbosity': 0, 'isUKS': 0, 'es_software': 'cp2k', 'use_multiprocessing': False, 'logfile_directory': os.getcwd()+'/all_logfiles', 'nac_algo': 0 } @@ -2440,10 +2440,12 @@ def run_step3_sd_nacs_libint(params): if params['isUKS']: sd_states_reindexed = step2_many_body.reindex_cp2k_sd_states( ks_homo_index_alpha, ks_orbital_indicies, sd_unique_basis, sd_format=1, - ks_beta_homo_index=ks_homo_index_beta) + ks_beta_homo_index=ks_homo_index_beta, active_space_num_occ_orbitals=params['active_space_num_occ_orbitals']) elif not params['isUKS']: sd_states_reindexed = step2_many_body.reindex_cp2k_sd_states( ks_homo_index, ks_orbital_indicies, - sd_unique_basis, sd_format=2 ) + sd_unique_basis, sd_format=2 , active_space_num_occ_orbitals=params['active_space_num_occ_orbitals']) + + # Some printings print('sd_unique_basis is:', sd_unique_basis) print('sd_states_reindexed is:', sd_states_reindexed) @@ -2452,6 +2454,21 @@ def run_step3_sd_nacs_libint(params): params['isnap'] = start_time params['fsnap'] = finish_time + if ks_homo_index - min(ks_active_space) < params['active_space_num_occ_orbitals']: + ks_active_space_ = [] + for i in range(ks_homo_index-params['active_space_num_occ_orbitals'], ks_homo_index): + ks_active_space_.append(i) + ks_active_space_.append(i+int(data_dim/2)) + for i in ks_active_space: + if i not in ks_active_space_: + ks_active_space_.append(i) + ks_active_space_ = np.sort(ks_active_space_).tolist() + ks_active_space = ks_active_space_ + params['active_space'] = ks_active_space_ + + + print('Flag ks_active_space:', ks_active_space) + # The flag for phase-correction algorithm apply_phase_correction = params['apply_phase_correction']