Merge pull request #233 from MohammadShakiba/devel

Bug fix for step 2 file saving and extension of states Slater determinants in step 3

Author: alexvakimov

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']:

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
 

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']