Included the mapping3 stuff in the cp2k and mopac workflows

Author: alexvakimov

src/libra_py/packages/mopac/methods.py

@@ -34,9 +34,11 @@
 from libra_py import units
 from libra_py import scan
 from libra_py import regexlib as rgl
+from libra_py import data_conv
 
 import libra_py.packages.cp2k.methods as CP2K_methods
 import libra_py.workflows.nbra.mapping2 as mapping2
+import libra_py.workflows.nbra.mapping3 as mapping3
 import libra_py.workflows.nbra.step3 as step3
 
 
@@ -516,6 +518,7 @@ def mopac_compute_adi(q, params, full_id):
     mopac_jobid = params[itraj]["mopac_jobid"] = F"timestep_{timestep}_traj_{itraj}"
     mopac_input_prefix = params[itraj]["mopac_input_prefix"]
     mopac_output_prefix = params[itraj]["mopac_output_prefix"]
+    active_space = params[itraj]["active_space"]
     dt = params[itraj]["dt"]
     do_Lowdin = params[itraj]["do_Lowdin"]
 
@@ -594,13 +597,26 @@ def mopac_compute_adi(q, params, full_id):
     ci_ovlp_curr = None
     if do_Lowdin:
         ident_curr = U_curr.T() * S_curr.real() * U_curr
-        ovlp_sd_curr = mapping2.ovlp_mat_arb(configs_curr, configs_curr, ident_curr, False).real()
+        # The original (older) approach
+        #ovlp_sd_curr = mapping2.ovlp_mat_arb(configs_curr, configs_curr, ident_curr, False).real()
+
+        # The new way:
+        # I don't like this way - we need to fix it later
+        ident_curr = data_conv.MATRIX2nparray(ident_curr).real  # MATRIX -> real np array
+        ovlp_sd_curr = mapping3.ovlp_mat_arb(configs_curr, configs_curr, ident_curr, active_space)  # do the calculations
+        ovlp_sd_curr = data_conv.nparray2MATRIX( ovlp_sd_curr )  # real np array -> MATRIX
+
         ovlp_sd_curr.scale(-1, 0, sqt2)
         ovlp_sd_curr.scale(0, -1, sqt2)
         ovlp_sd_curr.scale(0, 0, 0.5)
 
         ident_prev = U_prev.T() * S_prev.real() * U_prev
-        ovlp_sd_prev = mapping2.ovlp_mat_arb(configs_prev, configs_prev, ident_prev, False).real()
+        #ovlp_sd_prev = mapping2.ovlp_mat_arb(configs_prev, configs_prev, ident_prev, False).real()
+
+        ident_prev = data_conv.MATRIX2nparray(ident_prev).real  # MATRIX -> real np array
+        ovlp_sd_prev = mapping3.ovlp_mat_arb(configs_prev, configs_prev, ident_prev, active_space)  # do the calculations
+        ovlp_sd_prev = data_conv.nparray2MATRIX( ovlp_sd_prev )  # complex np array -> CMATRIX
+
         ovlp_sd_prev.scale(-1, 0, sqt2)
         ovlp_sd_prev.scale(0, -1, sqt2)
         ovlp_sd_prev.scale(0, 0, 0.5)
@@ -620,14 +636,25 @@ def mopac_compute_adi(q, params, full_id):
         ci_ovlp_curr = U_curr.T() * ovlp_ci_curr.real() * U_curr
     else:
         S_curr = MO_curr.T() * MO_curr
-        ovlp_sd_curr = mapping2.ovlp_mat_arb(configs_curr, configs_curr, S_curr, False).real()
+        #ovlp_sd_curr = mapping2.ovlp_mat_arb(configs_curr, configs_curr, S_curr, False).real()
+
+        S_curr = data_conv.MATRIX2nparray(S_curr).real  # MATRIX -> real np array
+        ovlp_sd_curr = mapping3.ovlp_mat_arb(configs_curr, configs_curr, S_curr, active_space)  # do the calculations
+        ovlp_sd_curr = data_conv.nparray2MATRIX( ovlp_sd_curr )  # complex np array -> MATRIX
+
         ovlp_sd_curr.scale(-1, 0, sqt2)
         ovlp_sd_curr.scale(0, -1, sqt2)
         ovlp_sd_curr.scale(0, 0, 0.5)
         ci_ovlp_curr = CI_curr.T() * ovlp_sd_curr * CI_curr
 
     # Time-overlap in the SD basis
-    time_ovlp_sd = mapping2.ovlp_mat_arb(configs_prev, configs_curr, mo_st, False).real()
+    #time_ovlp_sd = mapping2.ovlp_mat_arb(configs_prev, configs_curr, mo_st, False).real()
+
+    mo_st = data_conv.MATRIX2nparray(mo_st).real  # MATRIX -> real np array
+    time_ovlp_sd = mapping3.ovlp_mat_arb(configs_prev, configs_curr, mo_st, active_space)  # do the calculations
+    time_ovlp_sd = data_conv.nparray2MATRIX( time_ovlp_sd )  # complex np array -> MATRIX
+
+
     # Scaling to account for SAC prefactors:
     time_ovlp_sd.scale(-1, 0, sqt2)
     time_ovlp_sd.scale(0, -1, sqt2)

src/libra_py/workflows/nbra/mapping3.py

@@ -19,9 +19,8 @@
 
 """
 
-import os
-import sys
-import math
+import os, sys, math
+import numpy as np
 
 # Fisrt, we add the location of the library to test to the PYTHON path
 if sys.platform == "cygwin":
@@ -54,7 +53,7 @@ def ovlp_arb(_SD1, _SD2, S, active_space=None, verbose=False):
 
     """
 
-    nbasis = S.num_of_rows
+    nbasis = S.shape[0] #num_of_rows
 
     SD1, SD2 = [], []
 
@@ -87,7 +86,7 @@ def ovlp_arb(_SD1, _SD2, S, active_space=None, verbose=False):
             if (I * J) > 0:
                 i = abs(I) - 1
                 j = abs(J) - 1
-                s[indx_I, indx_J] = S.get(i, j).real
+                s[indx_I, indx_J] = S[i, j] # .real #S.get(i, j).real
             else:
                 s[indx_I, indx_J] = 0.0
 
@@ -100,7 +99,7 @@ def ovlp_arb(_SD1, _SD2, S, active_space=None, verbose=False):
 
 
 
-def ovlp_mat_arb(SD1, SD2, S, active_space):
+def ovlp_mat_arb(SD1, SD2, S, active_space=None):
     """Compute a matrix of overlaps in the SD basis
 
     Args:
@@ -126,12 +125,13 @@ def ovlp_mat_arb(SD1, SD2, S, active_space):
     """
 
     N, M = len(SD1), len(SD2)
-    res = CMATRIX(N, M)
+    res = np.zeros( (N, M), dtype=np.float64 ) #CMATRIX(N, M)
 
     for n in range(0, N):
         for m in range(0, M):
-            val = ovlp_arb(SD1[n], SD2[m], S, active_space, 0) 
-            res.set(n, m, val)
+            res[n, m] = ovlp_arb(SD1[n], SD2[m], S, active_space, 0) 
+            #res.set(n, m, val)
+            #res[n, m] = val
 
     return res
 

src/libra_py/workflows/nbra/step3.py

@@ -43,7 +43,7 @@
 elif sys.platform == "linux" or sys.platform == "linux2":
     from liblibra_core import *
 
-from . import mapping, step2_many_body, step3_many_body
+from . import mapping, mapping3, step2_many_body, step3_many_body
 import util.libutil as comn
 import libra_py.packages.cp2k.methods as CP2K_methods
 import libra_py.units as units
@@ -2717,7 +2717,8 @@ def compute_sd_overlaps_in_parallel(step, params):
 
     start_time = params['start_time']
     res_dir_1 = params['path_to_npz_files']
-    active_space = params['active_space']
+    active_space = params.get('active_space', None)
+    active_space2 = ( np.array(active_space) + 1 ).tolist() # conversion to a different convention
     sd_states_reindexed_sorted = params['sd_states_reindexed_sorted']
 
     t2 = time.time()
@@ -2739,24 +2740,27 @@ def compute_sd_overlaps_in_parallel(step, params):
     # Computing the overlaps for SDs
     t2 = time.time()
     if params['apply_orthonormalization']:
-        s_sd_1 = mapping.ovlp_mat_arb(
+        s_sd_1 = mapping3.ovlp_mat_arb(
             sd_states_reindexed_sorted[step],
             sd_states_reindexed_sorted[step],
             s_ks_1,
-            use_minimal=False,
-            use_mo_approach=False).real
-        s_sd_2 = mapping.ovlp_mat_arb(
+            active_space2)
+            #use_minimal=False,
+            #use_mo_approach=False).real
+        s_sd_2 = mapping3.ovlp_mat_arb(
             sd_states_reindexed_sorted[step + 1], 
             sd_states_reindexed_sorted[step + 1], 
             s_ks_2, 
-            use_minimal=False, 
-            use_mo_approach=False).real
+            active_space2)
+            #use_minimal=False, 
+            #use_mo_approach=False).real
 
-    st_sd = mapping.ovlp_mat_arb(sd_states_reindexed_sorted[step],
+    st_sd = mapping3.ovlp_mat_arb(sd_states_reindexed_sorted[step],
                                  sd_states_reindexed_sorted[step + 1],
                                  st_ks,
-                                 use_minimal=False,
-                                 use_mo_approach=False).real
+                                 active_space2)
+                                 #use_minimal=False,
+                                 #use_mo_approach=False).real
     # s_sd_1 = data_conv.MATRIX2nparray(s_sd_1)
     # s_sd_2 = data_conv.MATRIX2nparray(s_sd_2)
     # st_sd = data_conv.MATRIX2nparray(st_sd)