Added the new experimental decoherence scheme

Author: alexvakimov

src/dyn/Dynamics.cpp

@@ -1120,6 +1120,63 @@ void propagate_electronic(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonia
 
 }
 
+void renormalize_hopping_probabilities(vector< vector<double> >& g, vector< vector<vector<double> > >& coherence_factors, vector<int>& act_states){
+/**
+  Scale the hopping probabilities by the factors and ensure the sum of all 
+  hopping probabilities is still 1
+
+  g[itraj][i] - hopping probability from the current active state for the trajectory itraj to the state i
+
+  act_states[itraj] - the active state for the trajectory itraj
+
+  coherence_factors[itraj][i][j] - the scaling of the hopping probability from state i to state j for the trajectory itraj
+
+
+  Result: the g vector is modified
+*/
+  int ntraj = g.size();
+  int nstates = coherence_factors[0].size();
+
+  for(int itraj=0; itraj<ntraj; itraj++){
+    int a = act_states[itraj];
+    double sum = 0.0;
+    for(int j=0;j<nstates; j++){
+      if(j!=a){ 
+        double val = g[itraj][j] * coherence_factors[itraj][a][j];
+        sum += val;
+        g[itraj][j] =  val;
+      } 
+    }// for j
+    // The self-element: 
+    g[itraj][a] = 1.0 - sum;
+    if(g[itraj][a]<0){  g[itraj][a] = 0.0; }
+
+  }// for itraj  
+
+}
+
+void reset_coherence_factors(vector< vector<vector<double> > >& coherence_factors, vector<int>& act_states, vector<int>& old_states){
+
+  int ntraj = coherence_factors.size();
+  int nstates = coherence_factors[0].size();
+
+  for(int itraj=0; itraj<ntraj; itraj++){
+    if( act_states[itraj] != old_states[itraj]){
+      // If the hop has happened for this trajectory, we start 
+      // new evolution of the wavepackets on all the surfaces
+      // so we have to reset all the coherence factors to 1.0
+      for(int i=0; i<nstates; i++){
+        for(int j=0; j<nstates; j++){
+          coherence_factors[itraj][i][j] = 1.0;
+        }// for j
+      }// for i
+
+    }// if
+  }// for itraj
+
+}
+
+
 void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
               nHamiltonian& ham, nHamiltonian& ham_aux, bp::object py_funct, bp::dict params,  Random& rnd,
               vector<Thermostat>& therm){
@@ -1151,7 +1208,7 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
 
 //  cout<<"In compute_dynamics\n";
   //======== General variables =======================
-  int i, cdof, traj, dof, idof, ntraj1;
+  int i, j, cdof, traj, dof, idof, ntraj1;
 
   //========= Control parameters variables ===========
   dyn_control_params prms;
@@ -1378,7 +1435,7 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
   }
 
   else if(prms.decoherence_times_type==5){
-    decoherence_rates = Gu_Franco(prms, *dyn_var.ampl_adi);
+    decoherence_rates = Gu_Franco(prms, *dyn_var.ampl_adi); 
   }
 
   ///== Optionally, apply the dephasing-informed correction ==
@@ -1435,7 +1492,19 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
   }
 
   // DISH, rev2023
-  if(prms.decoherence_algo==7){  dish_rev2023(dyn_var, ham,  decoherence_rates, prms, rnd);   }
+  else if(prms.decoherence_algo==7){  dish_rev2023(dyn_var, ham,  decoherence_rates, prms, rnd);   }
+
+  // Simple decoherence
+  else if(prms.decoherence_algo==9){
+    for(traj=0; traj<ntraj; traj++){
+      for(i=0; i<nadi; i++){
+        for(j=0; j<nadi; j++){
+          double argg = decoherence_rates[traj].get(i,j) * prms.dt;
+          dyn_var.coherence_factors[traj][i][j] *= exp( - argg * argg );
+        }// for j
+      }// for i
+    }// for traj
+  }// simple decoherence
 
 
   // Update amplitudes and density matrices in response to decoherence corrections
@@ -1467,6 +1536,11 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
       vector< vector<double> > g;
       g = hop_proposal_probabilities(prms, dyn_var, ham, ham_aux);
 
+      // simple decoherence correction
+      if(prms.decoherence_algo==9){   // intended only for adiabatic states for now 
+        renormalize_hopping_probabilities(g, dyn_var.coherence_factors, dyn_var.act_states); 
+      }
+
       // Propose new discrete states for all trajectories
       vector<int> prop_states(ntraj, 0);
       if(prms.rep_sh==1){
@@ -1511,7 +1585,11 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
                                     prms.instantaneous_decoherence_variant, prms.collapse_option);
         }
         else{ cout<<"ERROR: Instantaneous Decoherence requires rep_tdse = 1\nExiting now...\n"; exit(0); }
-      }// algo == 6
+      }// algo == 8
+
+      else if(prms.decoherence_algo==9){ // simple decoherence method
+        reset_coherence_factors(dyn_var.coherence_factors, act_states, old_states);   
+      }// algo == 9 
 
     }
     // DISH - the old one
@@ -1522,6 +1600,7 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
     }// DISH
 
 
+
     //====================== Momenta adjustment after successful/frustrated hops ===================
     // Velocity rescaling: however here we may be changing velocities
     if(prms.rep_sh==1){

src/dyn/dyn_control_params.h

@@ -292,15 +292,15 @@ class dyn_control_params{
   int convergence;
 
   /**
-  The maximum number of hops that an be attempted before either choosing the identity or exiting in stochastic reordering algorithm 3. 
+    The maximum number of hops that an be attempted before either choosing the identity or exiting in stochastic reordering algorithm 3. 
   */
   int max_number_attempts;
 
   /**
-  The probability threshold for stochastic state reordering algorithm. 
-  If a probability for a multi-state stransition is below this value, it will be disregarded and set to 0
-  The rest of the probabilities will be renormalized
-  Default: 0.0 
+    The probability threshold for stochastic state reordering algorithm. 
+    If a probability for a multi-state stransition is below this value, it will be disregarded and set to 0
+    The rest of the probabilities will be renormalized
+    Default: 0.0 
   */
   double min_probability_reordering; 
 
@@ -513,6 +513,7 @@ class dyn_control_params{
       - 6: MQCXF
       - 7: DISH, rev2023
       - 8: diabatic IDA, experimental
+      - 9: simple decoherence, experimental
 
   */
   double decoherence_algo;

src/dyn/dyn_variables.cpp

@@ -119,6 +119,9 @@ dyn_variables::dyn_variables(int _ndia, int _nadi, int _ndof, int _ntraj){
   ///================= KC-RPMD ====================
   kcrpmd_vars_status = 0;
 
+  ///================= simple decoherence ====================
+  simple_decoherence_vars_status = 0;
+
 }
 
 
@@ -313,6 +316,18 @@ void dyn_variables::allocate_kcrpmd(){
 }// allocate_kcrpmd
 
 
+void dyn_variables::allocate_simple_decoherence(){
+
+  if(simple_decoherence_vars_status==0){
+
+    coherence_factors = std::vector< std::vector< std::vector<double> >  >(ntraj, std::vector< std::vector<double> >(nadi, std::vector<double>(nadi, 1.0) ) );
+    simple_decoherence_vars_status = 1;
+  }
+
+}
+
+
+
 dyn_variables::dyn_variables(const dyn_variables& x){     
   //cout<<"dyn_variables copy constructor\n";
   int itraj, idof;
@@ -470,6 +485,10 @@ dyn_variables::dyn_variables(const dyn_variables& x){
 
   }// if KCRPMD vars
 
+  if(x.simple_decoherence_vars_status == 1 ){
+    coherence_factors = x.coherence_factors;
+  }
+
 }// dyn_variables cctor
 
 
@@ -629,6 +648,12 @@ dyn_variables::~dyn_variables(){
     kcrpmd_vars_status = 0;
   }
 
+  if(simple_decoherence_vars_status==1){
+    coherence_factors.clear();
+
+    simple_decoherence_vars_status = 0;
+  }
+
 }
 
 

src/dyn/dyn_variables.h

@@ -552,12 +552,20 @@ class dyn_variables{
   int timestep; 
 
 
-  ///=============== For new decoherence method ==============
+  ///=============== For new decoherence method (let's call it simple_decoherence ) ==============
+  /**
+    Status of the new decoherence method vars
+
+    0 - not allocated;
+    1 - allocated 
+  */
+  int simple_decoherence_vars_status;
+
   /**
    Integrated exp( -(dt/tau(t))**2 ) over many timesteps 
-   Reset to 1 at every accepted hop
+   Reset to 1 at every accepted hop, separate for each trajectory
   */
-  vector< vector<double> > coherence_factors;
+  vector< vector< vector<double> > > coherence_factors;
 
 
   ///====================== In dyn_variables.cpp =====================
@@ -573,6 +581,8 @@ class dyn_variables{
   void allocate_mqcxf();
   void allocate_qtsh();
   void allocate_kcrpmd();
+  void allocate_simple_decoherence();
+
 
   dyn_variables(int _ndia, int _nadi, int _ndof, int _ntraj);
   dyn_variables(const dyn_variables& x); 

src/dyn/libdyn.cpp

@@ -248,6 +248,7 @@ void export_dyn_variables_objects(){
       .def_readwrite("tcnbra_ekin", &dyn_variables::tcnbra_ekin)
       .def_readwrite("qtsh_vars_status", &dyn_variables::qtsh_vars_status)
       .def_readwrite("kcrpmd_vars_status", &dyn_variables::kcrpmd_vars_status)
+      .def_readwrite("simple_decoherence_vars_status", &dyn_variables::simple_decoherence_vars_status)
 
       .def("set_parameters", expt_set_parameters_v1)
 
@@ -262,6 +263,7 @@ void export_dyn_variables_objects(){
       .def("allocate_mqcxf", &dyn_variables::allocate_mqcxf)
       .def("allocate_qtsh", &dyn_variables::allocate_qtsh)
       .def("allocate_kcrpmd", &dyn_variables::allocate_kcrpmd)
+      .def("allocate_simple_decoherence", &dyn_variables::allocate_simple_decoherence)
 
       .def("set_q", &dyn_variables::set_q)
       .def("set_p", &dyn_variables::set_p)

src/libra_py/dynamics/tsh/compute.py

@@ -988,6 +988,9 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd):
         dyn_var.allocate_shxf()
     elif decoherence_algo == 6:  # MQCXF
         dyn_var.allocate_mqcxf()
+    elif decoherence_algo == 9:  # simple decoherence
+        dyn_var.allocate_simple_decoherence()
+
     if tsh_method == 5 or decoherence_algo == 7:  # DISH or DISH_rev2023
         dyn_var.allocate_dish()
     if tsh_method == 7 or tsh_method == 8 or tsh_method == 9:  # FSSH2 or FSSH3 or original GFSH