Added original GFSH

Author: alexvakimov

README.md

@@ -225,7 +225,32 @@ FIND_PACKAGE(Python3 3.6 REQUIRED COMPONENTS Development)
 
 This will force the make to search within the libra environment in the specified locations where you know the files exist. 
 
- 
+
+### 2. 4/17/2025 (Alexey Akimov)
+
+A good way to setup the conda environment to have Boost and Python version consistent is this:
+
+
+    conda install -c conda-forge boost=1.82 python=3.10
+
+
+### 3. 4/17/2025 (Alexey Akimov)
+
+Another useful recipe for setting up jupyter notebook specific to a selected Conda environment:
+
+Step 1: Activate the environment
+
+    conda activate libra
+
+
+Step 2: Install ipykernel and register the kernel
+
+
+    conda install ipykernel
+    python -m ipykernel install --user --name=libra --display-name "Python (libra)"
+
+
+Now, in Jupyter, you'll see a new kernel called "Python (libra)". Select that in your notebook.
 
 
 

src/dyn/Dynamics.cpp

@@ -1441,16 +1441,17 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
   // FSSH, GFSH, MSSH, LZ, ZN, DISH, MASH, FSSH2, FSSH3
   else if(prms.tsh_method == 0 || prms.tsh_method == 1 || prms.tsh_method == 2 || prms.tsh_method == 3 
        || prms.tsh_method == 4 || prms.tsh_method == 5 || prms.tsh_method == 6 || prms.tsh_method == 7
-       || prms.tsh_method == 8 ){
+       || prms.tsh_method == 8 || prms.tsh_method == 9 ){
 
 
     vector<int> old_states(dyn_var.act_states); 
     vector<int> old_states_dia(dyn_var.act_states_dia); 
     //========================== Hop proposal and acceptance ================================
 
-    // FSSH (0), GFSH (1), MSSH (2), LZ(3), ZN (4), MASH(6), FSSH2(7), FSSH3(8)
+    // FSSH (0), GFSH (1), MSSH (2), LZ(3), ZN (4), MASH(6), FSSH2(7), FSSH3(8), GFSH original (9)
     if(prms.tsh_method == 0 || prms.tsh_method == 1 || prms.tsh_method == 2 || prms.tsh_method == 3  
-    || prms.tsh_method == 4 || prms.tsh_method == 6 || prms.tsh_method == 7 || prms.tsh_method == 8){
+    || prms.tsh_method == 4 || prms.tsh_method == 6 || prms.tsh_method == 7 || prms.tsh_method == 8
+    || prms.tsh_method == 9){
 
       /// Compute hop proposal probabilities from the active state of each trajectory to all other states 
       /// of that trajectory
@@ -1532,7 +1533,7 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
     // Update vib Hamiltonian to reflect the change of the momentum
     update_Hamiltonian_variables(prms, dyn_var, ham, ham_aux, py_funct, params, 1); 
 
-  }// tsh_method == 0, 1, 2, 3, 4, 5, 6, 7, 8
+  }// tsh_method == 0, 1, 2, 3, 4, 5, 6, 7, 8, 9
 
   else{   cout<<"tsh_method == "<<prms.tsh_method<<" is undefined.\nExiting...\n"; exit(0);  }
 
@@ -1543,7 +1544,7 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
   dyn_var.update_density_matrix(prms);
 
 
-  // Saves the current density matrix into the previous - needed for FSSH2
+  // Saves the current density matrix into the previous - needed for FSSH2 and GFSH (original)
   dyn_var.save_curr_dm_into_prev();
 
 

src/dyn/dyn_control_params.h

@@ -323,6 +323,7 @@ class dyn_control_params{
       - 6: MASH
       - 7: FSSH2
       - 8: FSSH3
+      - 9: GFSH (original)
   */
   int tsh_method;
 

src/dyn/dyn_hop_proposal.cpp

@@ -370,6 +370,10 @@ MATRIX hopping_probabilities_gfsh(dyn_control_params& prms, CMATRIX& Coeff, CMAT
 
 
 
+
+
+
+
 vector<double> hopping_probabilities_gfsh(dyn_control_params& prms, CMATRIX& denmat, CMATRIX& Hvib, int act_state_indx){
 /**
   \brief Compute the GFSH surface hopping probabilities for a single trajectory
@@ -468,6 +472,100 @@ vector<double> hopping_probabilities_gfsh(dyn_control_params& prms, CMATRIX& den
 
 
 
+
+vector<double> hopping_probabilities_gfsh_orig(dyn_control_params& prms, CMATRIX& denmat, CMATRIX& denmat_old, int act_state_indx){
+/**
+  \brief Compute the GFSH surface hopping probabilities for a single trajectory
+
+  Abbreviation: GFSH - global flux surface hopping
+  References: 
+  (1) Wang, L.; Trivedi, D.; Prezhdo, O. V. Global Flux Surface Hopping Approach for Mixed Quantum-Classical Dynamics. J. Chem. Theory Comput. 2014, 10, 3598–3605.
+  (2) Akimov, A. V. Libra: An Open-Source “Methodology Discovery” Library for Quantum and Classical Dynamics Simulations.
+  J. Comput. Chem. 2016, 37, 1626–1649.
+
+  \param[in] key parameters needed for this type of calculations
+    - dt - integration timestep [a.u.]
+    - Temperature - temperature [ K ]
+    - use_boltz_factor - whether to scale the computed probabilities by a Boltzmann factor
+  \param[in] denmat - [nstates x nstates] - current density matrix
+  \param[in] denmat_old - [nstates x nstates] - previous density matrix
+  \param[in] act_state_indx - index of the initial state
+
+  Returns: A nstates-vector of hopping probabilities to all states from the current active state
+
+
+*/
+
+  const double kb = 3.166811429e-6; // Hartree/K
+  int i,j,k;
+  double sum,g_ij,argg;
+
+  double dt = prms.dt;
+  double T = prms.Temperature;
+  int use_boltz_factor = prms.use_boltz_factor;
+
+
+  int nstates = denmat.n_rows;
+  vector<double> g(nstates, 0.0);
+
+  CMATRIX denmat_dot(nstates, nstates);
+  denmat_dot = denmat - denmat_old;
+
+
+  // compute a_kk and a_dot_kk
+  vector<double> a(nstates,0.0);
+  vector<double> a_dot(nstates,0.0);
+  double norm = 0.0; // normalization factor
+
+  for(i=0;i<nstates;i++){
+    a[i] = denmat.get(i,i).real();
+    a_dot[i] = denmat_dot.get(i,i).real();
+
+    if(a_dot[i]<0.0){ norm += a_dot[i]; } // total rate of population decrease in all decaying states
+
+  }// for i
+
+
+  // Now calculate the hopping probabilities
+  i = act_state_indx;
+  double sumg = 0.0;
+
+  for(j=0;j<nstates;j++){
+
+    if(j!=i){  // off-diagonal = probabilities to hop to other states
+
+      if(a[i]<1e-12){  g[j] = 0.0; }  // since the initial population is almost zero, so no need for hops
+      else{
+
+        if( fabs(norm) > 1e-12 ){ g[j] = (a_dot[j]/a[i]) * a_dot[i] / norm;  }
+
+        // since norm is negative, than this condition means that a_dot[i] and a_dot[j] have same signs
+        // which is bad - so no transitions are assigned
+        if(g[j]<0.0){  g[j] = 0.0;   }
+
+        // here we have opposite signs of a_dot[i] and a_dot[j], but this is not enough yet
+        else{  
+          if(a_dot[i]<0.0 && a_dot[j]>0.0){ ;; } // this is out transition probability, but it is already computed
+          else{  g[j] = 0.0; } // wrong transition
+        }
+      }// a[i]>1e-12
+
+      sumg += g[j];
+    }
+  }// for j
+
+  g[i] = 1.0 - sumg;  // probability to stay in state i
+
+  if(g[i]<0.0){  g[i] = 0.0; }
+
+  return g;
+
+}// gfsh - original version 
+
+
+
+
+
 vector<double> hopping_probabilities_fssh2(dyn_control_params& prms, CMATRIX& denmat, CMATRIX& denmat_old, int act_state_indx){
 /**
   \brief This function computes the surface hopping probabilities according to Leonardo Araujo's hopping probability 
@@ -1262,10 +1360,24 @@ nHamiltonian& ham, nHamiltonian& ham_prev){
         CMATRIX dm_prev_trans(*dyn_var.dm_dia_prev[traj]);
         g[traj] = hopping_probabilities_fssh3(prms, dm, dm_prev_trans, dyn_var.act_states_dia[traj], dyn_var.fssh3_errors[traj]);
       }
-    }
+    } // FSSH-3
+	
+	else if(prms.tsh_method == 9){ // GFSH - original
+ 
+      if(prms.rep_sh==1){
+        CMATRIX& dm_prev = *dyn_var.dm_adi_prev[traj];
+        if(prms.rep_tdse==2){ dm = *dyn_var.dm_dia_prev[traj]; }
+        g[traj] = hopping_probabilities_gfsh_orig(prms, dm, dm_prev, dyn_var.act_states[traj]);
+      }
+      else if(prms.rep_sh==0){
+        CMATRIX& dm_prev = *dyn_var.dm_dia_prev[traj];
+        g[traj] = hopping_probabilities_gfsh_orig(prms, dm, dm_prev, dyn_var.act_states_dia[traj]);
+      }
+    }// GFSH original
+		
 
     else{
-      cout<<"Error in tsh1: tsh_method can be -1, 0, 1, 2, 3, 4, 5, 6, 7, or 8. Other values are not defined\n";
+      cout<<"Error in tsh1: tsh_method can be -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9. Other values are not defined\n";
       cout<<"Exiting...\n";
       exit(0);
     }

src/dyn/dyn_hop_proposal.h

@@ -44,6 +44,8 @@ vector<double> hopping_probabilities_fssh(dyn_control_params& prms, CMATRIX& den
 MATRIX hopping_probabilities_gfsh(dyn_control_params& prms, CMATRIX& Coeff, CMATRIX& Hvib);
 vector<double> hopping_probabilities_gfsh(dyn_control_params& prms, CMATRIX& denmat, CMATRIX& Hvib, int atc_state_indx);
 
+vector<double> hopping_probabilities_gfsh_orig(dyn_control_params& prms, CMATRIX& denmat, CMATRIX& denmat_old, int act_state_indx);
+
 vector<double> hopping_probabilities_fssh2(dyn_control_params& prms, CMATRIX& denmat, CMATRIX& denmat_old, int act_state_indx);
 
 MATRIX hopping_probabilities_mssh(dyn_control_params& prms, CMATRIX& Coeff, CMATRIX& Hvib);

src/dyn/libdyn.cpp

@@ -639,6 +639,9 @@ void export_dyn_hop_proposal_objects(){
   (dyn_control_params& prms, CMATRIX& denmat, CMATRIX& Hvib, int act_state_indx) = &hopping_probabilities_gfsh;
   def("hopping_probabilities_gfsh", expt_hopping_probabilities_gfsh_v2);
 
+  vector<double> (*expt_hopping_probabilities_gfsh_orig_v1)
+  (dyn_control_params& prms, CMATRIX& denmat, CMATRIX& denmat_old, int act_state_indx) = &hopping_probabilities_gfsh_orig;
+  def("hopping_probabilities_gfsh_orig", expt_hopping_probabilities_gfsh_orig_v1);
 
   vector<double> (*expt_hopping_probabilities_fssh2_v1)
   (dyn_control_params& prms, CMATRIX& denmat, CMATRIX& denmat_old, int act_state_indx) = &hopping_probabilities_fssh2;

src/libra_py/dynamics/tsh/compute.py

@@ -990,7 +990,7 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd):
         dyn_var.allocate_mqcxf()
     if tsh_method == 5 or decoherence_algo == 7:  # DISH or DISH_rev2023
         dyn_var.allocate_dish()
-    if tsh_method == 7 or tsh_method == 8:  # FSSH2 or FSSH3
+    if tsh_method == 7 or tsh_method == 8 or tsh_method == 9:  # FSSH2 or FSSH3 or original GFSH
         dyn_var.allocate_fssh2()
         dyn_var.save_curr_dm_into_prev()