Merge branch 'devel'

Author: alexvakimov

README.md

@@ -138,6 +138,11 @@ Next, all what we actually need:
     conda install -y -c psi4/label/dev libint2=2.7.1
     conda install -y -c anaconda h5py gmp
 
+
+Installation instruction of Scikit-learn from its official website:
+
+    pip install -U scikit-learn
+
 
     >
     >  YES - IT GOT SMALLER AND MORE COMPACT !

src/dyn/Dynamics.cpp

@@ -1295,6 +1295,9 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
   // Recompute diabatic/adiabatic states, time-overlaps, NAC, Hvib, etc. in response to change of q
   update_Hamiltonian_variables(prms, dyn_var, ham, ham_aux, py_funct, params, 0);
 
+  // Copy diabatic-to-adiabatic basis transformation to the dynamical variable
+  dyn_var.update_basis_transform(ham);
+
   // Recompute the orthogonalized reprojection matrices, stored in dyn_var.proj_adi
   // this calculaitons used ham.children[i].time_overlap matrix, updated in the previous step
   update_proj_adi(prms, dyn_var, ham); 
@@ -1326,15 +1329,20 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
   }
 
   // Recompute density matrices in response to the updated amplitudes  
-  dyn_var.update_amplitudes(prms, ham);
-  dyn_var.update_density_matrix(prms, ham, 1); 
+//  dyn_var.update_amplitudes(prms, ham);
+  dyn_var.update_amplitudes(prms);
+//  dyn_var.update_density_matrix(prms, ham, 1); 
+  dyn_var.update_density_matrix(prms);
 
   vector<int> old_states( dyn_var.act_states);
 
   // In the interval [t, t + dt], we may have experienced the basis reordering, so we need to 
   // change the active adiabatic state
-  if(prms.tsh_method != 3 && prms.tsh_method != 4 ){  // Don't update states based on amplitudes, in the LZ method
-    dyn_var.update_active_states(1, 0); // 1 - forward; 0 - only active state
+  if(prms.tsh_method == 3 or prms.tsh_method == 4 ){  
+  // Don't update states based on amplitudes, in the LZ method
+    ;;
+  }
+  else{  dyn_var.update_active_states(1, 0); // 1 - forward; 0 - only active state
   }
 
   // For now, this function also accounts for the kinetic energy adjustments to reflect the adiabatic evolution
@@ -1406,6 +1414,10 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
     decoherence_rates = schwartz_2(prms, ham, *prms.schwartz_decoherence_inv_alpha); 
   }
 
+  else if(prms.decoherence_times_type==4){
+    decoherence_rates = schwartz_1(prms, *dyn_var.ampl_adi, *dyn_var.p, ham, *prms.schwartz_interaction_width); 
+  }
+
   ///== Optionally, apply the dephasing-informed correction ==
   if(prms.dephasing_informed==1){
     Eadi = get_Eadi(ham); 
@@ -1438,7 +1450,7 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
     if(prms.rep_tdse==1){
       p = *dyn_var.p;
       //cout<<"p before mfsd\n"; dyn_var.p->show_matrix();
-      *dyn_var.ampl_adi = mfsd(p, *dyn_var.ampl_adi, *dyn_var.iM, prms.dt, decoherence_rates, ham, rnd, prms.isNBRA);
+      *dyn_var.ampl_adi = mfsd(p, *dyn_var.ampl_adi, *dyn_var.iM, prms.dt, dyn_var.act_states, decoherence_rates, ham, rnd, prms.isNBRA);
       *dyn_var.p = p;
       //cout<<"p after mfsd\n"; dyn_var.p->show_matrix();
 
@@ -1464,25 +1476,28 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
 
 
   // Update amplitudes and density matrices in response to decoherence corrections
-  dyn_var.update_amplitudes(prms, ham);
-  dyn_var.update_density_matrix(prms, ham, 1);
+//  dyn_var.update_amplitudes(prms, ham);
+  dyn_var.update_amplitudes(prms);
+//  dyn_var.update_density_matrix(prms, ham, 1);
+  dyn_var.update_density_matrix(prms);
 
 
   //************************************ TSH options ****************************************
   // Adiabatic dynamics
   if(prms.tsh_method==-1){ ;; } 
 
-  // FSSH, GFSH, MSSH, LZ, ZN, DISH, MASH, FSSH2
+  // 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 == 4 || prms.tsh_method == 5 || prms.tsh_method == 6 || prms.tsh_method == 7
+       || prms.tsh_method == 8 ){
 
 
     vector<int> old_states(dyn_var.act_states); 
     //========================== Hop proposal and acceptance ================================
 
-    // FSSH (0), GFSH (1), MSSH (2), LZ(3), ZN (4), MASH(6), FSSH2(7)
+    // FSSH (0), GFSH (1), MSSH (2), LZ(3), ZN (4), MASH(6), FSSH2(7), FSSH3(8)
     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 == 4 || prms.tsh_method == 6 || prms.tsh_method == 7 || prms.tsh_method == 8){
 
       /// Compute hop proposal probabilities from the active state of each trajectory to all other states 
       /// of that trajectory
@@ -1549,15 +1564,17 @@ 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
+  }// tsh_method == 0, 1, 2, 3, 4, 5, 6, 7, 8
 
   else{   cout<<"tsh_method == "<<prms.tsh_method<<" is undefined.\nExiting...\n"; exit(0);  }
 
 
   // Update the amplitudes and DM, in response to state hopping and other changes in the TSH part
   // so that we have them consistent in the output
-  dyn_var.update_amplitudes(prms, ham); 
-  dyn_var.update_density_matrix(prms, ham, 1);
+//  dyn_var.update_amplitudes(prms, ham); 
+  dyn_var.update_amplitudes(prms);
+//  dyn_var.update_density_matrix(prms, ham, 1);
+  dyn_var.update_density_matrix(prms);
 
 
   // Saves the current density matrix into the previous - needed for FSSH2

src/dyn/dyn_control_params.cpp

@@ -1,5 +1,5 @@
 /*********************************************************************************
-* Copyright (C) 2019-2022 Alexey V. Akimov
+* Copyright (C) 2019-2024 Alexey V. Akimov
 *
 * This file is distributed under the terms of the GNU General Public License
 * as published by the Free Software Foundation, either version 3 of
@@ -64,14 +64,27 @@ dyn_control_params::dyn_control_params(){
   tsh_method = 0;
   hop_acceptance_algo = 0;
   momenta_rescaling_algo = 0;
+  use_Jasper_Truhlar_criterion = 1;
   use_boltz_factor = 0;
 
+  ///================= FSSH2 specific ====================
+  fssh2_revision = 0;
+
+  ///================= FSSH3 specific ====================
+  fssh3_size_option = 1;
+  fssh3_approach_option = 0;
+  fssh3_decomp_option = 3;
+  fssh3_dt = 0.001;
+  fssh3_max_steps = 1000;
+  fssh3_err_tol = 1e-7;
+
   ///================= Decoherence options =========================================
   decoherence_algo = -1; 
   sdm_norm_tolerance = 0.0;
   dish_decoherence_event_option = 1;
   decoherence_times_type = -1;
   schwartz_decoherence_inv_alpha = NULL;
+  schwartz_interaction_width = NULL;
   decoherence_C_param = 1.0;
   decoherence_eps_param = 0.1;
   dephasing_informed = 0;
@@ -151,8 +164,20 @@ dyn_control_params::dyn_control_params(const dyn_control_params& x){
   tsh_method = x.tsh_method;
   hop_acceptance_algo = x.hop_acceptance_algo;
   momenta_rescaling_algo = x.momenta_rescaling_algo;
+  use_Jasper_Truhlar_criterion = x.use_Jasper_Truhlar_criterion;
   use_boltz_factor = x.use_boltz_factor;
 
+  ///================= FSSH2 specific ====================
+  fssh2_revision = x.fssh2_revision;
+ 
+  ///================= FSSH3 specific ====================
+  fssh3_size_option = x.fssh3_size_option;
+  fssh3_approach_option = x.fssh3_approach_option;
+  fssh3_decomp_option = x.fssh3_decomp_option;
+  fssh3_dt = x.fssh3_dt;
+  fssh3_max_steps = x.fssh3_max_steps;
+  fssh3_err_tol = x.fssh3_err_tol;
+
   ///================= Decoherence options =========================================
   decoherence_algo = x.decoherence_algo; 
   sdm_norm_tolerance = x.sdm_norm_tolerance;
@@ -207,6 +232,9 @@ dyn_control_params::dyn_control_params(const dyn_control_params& x){
 
   schwartz_decoherence_inv_alpha = new MATRIX( x.schwartz_decoherence_inv_alpha->n_rows, x.schwartz_decoherence_inv_alpha->n_cols );
   *schwartz_decoherence_inv_alpha = *x.schwartz_decoherence_inv_alpha;
+  
+  schwartz_interaction_width = new MATRIX( x.schwartz_interaction_width->n_rows, x.schwartz_interaction_width->n_cols );
+  *schwartz_interaction_width = *x.schwartz_interaction_width;
 
   thermally_corrected_nbra = x.thermally_corrected_nbra;
   total_energy = x.total_energy;
@@ -224,6 +252,7 @@ dyn_control_params::~dyn_control_params() {
   delete decoherence_rates;  
   delete ave_gaps;
   delete schwartz_decoherence_inv_alpha;
+  delete schwartz_interaction_width;
 }
 
 
@@ -326,8 +355,20 @@ void dyn_control_params::set_parameters(bp::dict params){
     else if(key=="tsh_method") { tsh_method = bp::extract<int>(params.values()[i]);  }
     else if(key=="hop_acceptance_algo") { hop_acceptance_algo = bp::extract<int>(params.values()[i]);  }
     else if(key=="momenta_rescaling_algo"){ momenta_rescaling_algo = bp::extract<int>(params.values()[i]);  }
+    else if(key=="use_Jasper_Truhlar_criterion"){ use_Jasper_Truhlar_criterion = bp::extract<int>(params.values()[i]); }
     else if(key=="use_boltz_factor"){ use_boltz_factor = bp::extract<int>(params.values()[i]);  }
 
+    ///================= FSSH2 specific ====================
+    else if(key=="fssh2_revision"){ fssh2_revision = bp::extract<int>(params.values()[i]);  }
+
+    ///================= FSSH3 specific ====================
+    else if(key=="fssh3_size_option"){ fssh3_size_option = bp::extract<int>(params.values()[i]);  }
+    else if(key=="fssh3_approach_option"){ fssh3_approach_option = bp::extract<int>(params.values()[i]);  }
+    else if(key=="fssh3_decomp_option"){ fssh3_decomp_option = bp::extract<int>(params.values()[i]);  }
+    else if(key=="fssh3_dt"){ fssh3_dt = bp::extract<double>(params.values()[i]);  }
+    else if(key=="fssh3_max_steps"){ fssh3_max_steps = bp::extract<int>(params.values()[i]);  }
+    else if(key=="fssh3_err_tol"){ fssh3_err_tol = bp::extract<double>(params.values()[i]);  }
+
     ///================= Decoherence options =========================================
     else if(key=="decoherence_algo"){ decoherence_algo = bp::extract<int>(params.values()[i]); }
     else if(key=="sdm_norm_tolerance"){ sdm_norm_tolerance = bp::extract<double>(params.values()[i]); }
@@ -340,6 +381,13 @@ void dyn_control_params::set_parameters(bp::dict params){
         for(int b=0;b<x.n_cols;b++){ schwartz_decoherence_inv_alpha->set(a, b, x.get(a,b));   }
       } 
     }
+    else if(key=="schwartz_interaction_width"){ 
+      MATRIX x( bp::extract<MATRIX>(params.values()[i]) );
+      schwartz_interaction_width = new MATRIX(x.n_rows, x.n_cols);      
+      for(int a=0;a<x.n_rows;a++){
+        for(int b=0;b<x.n_cols;b++){ schwartz_interaction_width->set(a, b, x.get(a,b));   }
+      } 
+    }
     else if(key=="decoherence_C_param"){ decoherence_C_param = bp::extract<double>(params.values()[i]); }
     else if(key=="decoherence_eps_param"){ decoherence_eps_param = bp::extract<double>(params.values()[i]); }
     else if(key=="dephasing_informed"){ dephasing_informed = bp::extract<int>(params.values()[i]); }

src/dyn/dyn_control_params.h

@@ -309,6 +309,7 @@ class dyn_control_params{
       - 5: DISH
       - 6: MASH
       - 7: FSSH2
+      - 8: FSSH3 - experimental
   */
   int tsh_method;
 
@@ -358,6 +359,23 @@ class dyn_control_params{
   */
   int momenta_rescaling_algo;
 
+
+  /**
+    A flag to turn on/off the Jasper-Truhlar criterion for reversal of velocities on frustrated hops.
+    According to: Jasper, A. W.; Truhlar, D. G. Chem. Phys. Lett. 2003, 369, 60− 67
+
+    Options:
+    
+      - 0: don't use this criterion (naive handling)
+      - 1: use it [ default ] - the velocities are reversed along the direction d_{a,j} if
+        a) (F_a * d_{a,j}) * (F_j * d_{a,j}) < 0 and b) (v * d_{a,j}) * (F_j * d_{a,j}) < 0 
+        where a - is the active state index;  Only in effect, if `momenta_rescaling_algo == 201`
+        
+
+  */
+  int use_Jasper_Truhlar_criterion;
+
+
   /**
     A flag to scale the proposed hopping probabilities by the
     Boltzmann factor. This is needed for the libra_py/workflows/nbra calculations, where the hop proposal 
@@ -371,6 +389,72 @@ class dyn_control_params{
    int use_boltz_factor;
 
 
+  //=========== FSSH2 options ==========
+  /**
+    Whether to use the revised FSSH2
+
+    Options:
+
+      - 0: use the FSSH2 as it was formulated in the original paper [ default ]
+      - 1: apply the revised version 
+  */
+  int fssh2_revision;
+
+
+  //=========== FSSH3 options ==========
+  /**
+    The size of the vectorized density matrix in equations to determine hopping probabilites/fluxes
+
+    Options:
+
+      - 0: N elements - only populations; the matrices are overdetermined
+      - 1: N^2 elements - first N elements are populations, then Re and Im parts of upper-triangular coherences
+                         that is rho_{0,1}, rho_{0,2}, ... rho_{0,N-1}, rho_{1,2}, ... rho_{1,N-1}, ... rho_{N-2,N-1}  [ default ]
+  */
+  int fssh3_size_option;
+
+  /**
+    The approach to determine the hopping probabilities:
+
+    Options:
+
+      - 0: based on master equation, rho(t+dt) = J * rho(t);  J matrix contains hopping probabilities directly  [ defualt ]
+      - 1: based on kinetic approach, drho/dt = J * rho; J matrix contains fluxes
+  */
+  int fssh3_approach_option;
+
+  /** 
+    The matrix decomposition method for solving the least-squares problem
+
+    Options:
+      - 0: bdcSvd
+      - 1: fullPivLu
+      - 2: fullPivHouseholderQr
+      - 3: completeOrthogonalDecomposition [ default ]
+  */
+  int fssh3_decomp_option;
+
+  
+  /**
+    The time-step of the optimization procedure in the FSSH3 calculations
+    Default: 0.001
+  */
+  double fssh3_dt;
+
+
+  /**
+    The maximal number of steps in the FSSH3 optimization step
+    Default: 1000
+  */
+  int fssh3_max_steps;
+
+
+  /**
+    FSSH3 error tolerance
+    Default: 1e-7
+  */
+  double fssh3_err_tol;
+
 
   ///===============================================================================
   ///================= Decoherence options =========================================
@@ -437,6 +521,12 @@ class dyn_control_params{
     computing decoherence rates [ default: NULL ]
   */
   MATRIX* schwartz_decoherence_inv_alpha;
+  
+  /**
+    MATRIX(ndof, 1) - the parameters for the spatial extent of NAC in
+    computing decoherence rates [ default: NULL ]
+  */
+  MATRIX* schwartz_interaction_width;
 
 
   /**

src/dyn/dyn_decoherence.h

@@ -70,8 +70,8 @@ CMATRIX bcsh(CMATRIX& Coeff, double dt, vector<int>& act_states, MATRIX& reversa
 
 
 // For MF-SD of Schwartz
-CMATRIX mfsd(MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector<MATRIX>& decoherence_rates, nHamiltonian& ham, Random& rnd, int isNBRA);
-CMATRIX mfsd(MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector<MATRIX>& decoherence_rates, nHamiltonian& ham, Random& rnd);
+CMATRIX mfsd(MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector<int>& act_st, vector<MATRIX>& decoherence_rates, nHamiltonian& ham, Random& rnd, int isNBRA);
+CMATRIX mfsd(MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector<int>& act_st, vector<MATRIX>& decoherence_rates, nHamiltonian& ham, Random& rnd);
 
 // Independent-trajectory XF methods
 void xf_hop_reset(dyn_variables& dyn_var, vector<int>& accepted_states, vector<int>& initial_states);
@@ -107,6 +107,7 @@ MATRIX coherence_intervals(CMATRIX& Coeff, MATRIX& rates);
 MATRIX coherence_intervals(CMATRIX& Coeff, vector<MATRIX>& rates);
 
 vector<MATRIX> schwartz_1(dyn_control_params& prms, CMATRIX& amplitudes, nHamiltonian& ham, MATRIX& inv_alp);
+vector<MATRIX> schwartz_1(dyn_control_params& prms, CMATRIX& amplitudes, MATRIX& p, nHamiltonian& ham, MATRIX& inv_alp);
 vector<MATRIX> schwartz_2(dyn_control_params& prms, nHamiltonian& ham, MATRIX& inv_alp);