Merge branch 'devel'

src/dyn/Energy_and_Forces.cpp

@@ -316,7 +316,8 @@ void update_forces(dyn_control_params& prms, dyn_variables& dyn_vars, nHamiltoni
          prms.electronic_integrator==12 ){ option = 1; } 
 
       option = 0;
-      *dyn_vars.f = ham.Ehrenfest_forces_adi(*dyn_vars.ampl_adi, 1, option, dyn_vars.proj_adi).real();
+      //*dyn_vars.f = ham.Ehrenfest_forces_adi(*dyn_vars.ampl_adi, 1, option, dyn_vars.proj_adi).real();
+      *dyn_vars.f = ham.Ehrenfest_forces_adi(*dyn_vars.ampl_adi, 1, option).real();
 /*
       CMATRIX& U = *ham.basis_transform;
       CMATRIX& C = *dyn_vars.ampl_adi;

src/dyn/dyn_control_params.cpp

@@ -52,7 +52,7 @@ dyn_control_params::dyn_control_params(){
   do_ssy = 0;
   do_phase_correction = 1;
   phase_correction_tol = 1e-3; 
-  state_tracking_algo = 2;
+  state_tracking_algo = -1;
   MK_alpha = 0.0;
   MK_verbosity = 0;
   convergence = 0;
@@ -79,10 +79,14 @@ dyn_control_params::dyn_control_params(){
   collapse_option = 0;
   decoherence_rates = NULL;
   ave_gaps = NULL;
-  wp_width = 0.1;
+  wp_width = NULL;
+  wp_v = NULL;
   coherence_threshold = 0.01;
+  e_mask = 0.0001;
   use_xf_force = 0;
   project_out_aux = 0;
+  tp_algo = 1;
+  use_td_width = 0;
 
   ///================= Entanglement of trajectories ================================
   entanglement_opt = 0;
@@ -104,6 +108,7 @@ dyn_control_params::dyn_control_params(){
   dt = 41.0;
   num_electronic_substeps = 1;
   electronic_integrator = 0;
+  ampl_transformation_method = 1; 
   assume_always_consistent = 0;
 
   thermally_corrected_nbra = 0;
@@ -158,10 +163,17 @@ dyn_control_params::dyn_control_params(const dyn_control_params& x){
   dephasing_informed = x.dephasing_informed;
   instantaneous_decoherence_variant = x.instantaneous_decoherence_variant; 
   collapse_option = x.collapse_option;
-  wp_width = x.wp_width;
+
+  wp_width = new MATRIX( x.wp_width->n_rows, x.wp_width->n_cols );
+  *wp_width = *x.wp_width;
+  wp_v = new MATRIX( x.wp_v->n_rows, x.wp_v->n_cols );
+  *wp_v = *x.wp_v;
   coherence_threshold = x.coherence_threshold;
+  e_mask = x.e_mask;
   use_xf_force = x.use_xf_force;
   project_out_aux = x.project_out_aux;
+  tp_algo = x.tp_algo;
+  use_td_width = x.use_td_width;
 
   ///================= Entanglement of trajectories ================================
   entanglement_opt = x.entanglement_opt;
@@ -184,6 +196,7 @@ dyn_control_params::dyn_control_params(const dyn_control_params& x){
   dt = x.dt;
   num_electronic_substeps = x.num_electronic_substeps;
   electronic_integrator = x.electronic_integrator;
+  ampl_transformation_method = x.ampl_transformation_method;
   assume_always_consistent = x. assume_always_consistent;
 
   decoherence_rates = new MATRIX(x.decoherence_rates->n_rows, x.decoherence_rates->n_cols);  
@@ -206,6 +219,8 @@ dyn_control_params::~dyn_control_params() {
 
   //cout<<"dyn_control_params destructor\n";
 
+  delete wp_width;
+  delete wp_v;
   delete decoherence_rates;  
   delete ave_gaps;
   delete schwartz_decoherence_inv_alpha;
@@ -215,7 +230,8 @@ dyn_control_params::~dyn_control_params() {
 void dyn_control_params::sanity_check(){
 
   ///=================== Options for state tracking ======================
-  if(state_tracking_algo==0 || state_tracking_algo==1 ||
+  if(state_tracking_algo==-1 ||
+     state_tracking_algo==0 || state_tracking_algo==1 ||
      state_tracking_algo==2 || state_tracking_algo==3 ||
      state_tracking_algo==32 || state_tracking_algo==33){ ; ; }
   else{
@@ -343,10 +359,26 @@ void dyn_control_params::set_parameters(bp::dict params){
         for(int b=0;b<x.n_cols;b++){ ave_gaps->set(a, b, x.get(a,b));   }
       } 
     }
-    else if(key=="wp_width"){ wp_width = bp::extract<double>(params.values()[i]); }
+    else if(key=="wp_width"){ 
+      MATRIX x( bp::extract<MATRIX>(params.values()[i]) );
+      wp_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++){ wp_width->set(a, b, x.get(a,b));   }
+      } 
+    }
+    else if(key=="wp_v"){ 
+      MATRIX x( bp::extract<MATRIX>(params.values()[i]) );
+      wp_v = 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++){ wp_v->set(a, b, x.get(a,b));   }
+      } 
+    }
     else if(key=="coherence_threshold"){ coherence_threshold = bp::extract<double>(params.values()[i]); }
+    else if(key=="e_mask"){ e_mask = bp::extract<double>(params.values()[i]); }
     else if(key=="use_xf_force"){ use_xf_force = bp::extract<int>(params.values()[i]); }
     else if(key=="project_out_aux"){ project_out_aux = bp::extract<int>(params.values()[i]); }
+    else if(key=="tp_algo"){ tp_algo = bp::extract<int>(params.values()[i]); }
+    else if(key=="use_td_width"){ use_td_width = bp::extract<int>(params.values()[i]); }
 
     ///================= Entanglement of trajectories ================================
     else if(key=="entanglement_opt"){ entanglement_opt = bp::extract<int>(params.values()[i]); }
@@ -378,6 +410,7 @@ void dyn_control_params::set_parameters(bp::dict params){
     else if(key=="dt") { dt = bp::extract<double>(params.values()[i]);  }
     else if(key=="num_electronic_substeps") { num_electronic_substeps = bp::extract<int>(params.values()[i]);  }
     else if(key=="electronic_integrator"){ electronic_integrator = bp::extract<int>(params.values()[i]); }
+    else if(key=="ampl_transformation_method"){ ampl_transformation_method = bp::extract<int>(params.values()[i]); }
     else if(key=="assume_always_consistent"){  assume_always_consistent = bp::extract<int>(params.values()[i]); }
 
     else if(key=="thermally_corrected_nbra"){ thermally_corrected_nbra = bp::extract<int>(params.values()[i]); }

src/dyn/dyn_control_params.h

@@ -244,9 +244,10 @@ class dyn_control_params{
 
   /** 
     State tracking algorithm:
+      - -1: use LD approach, it includes phase correction too [ default ]
       - 0: no state tracking
       - 1: method of Kosuke Sato (may fail by getting trapped into an infinite loop)
-      - 2: Munkres-Kuhn (Hungarian) algorithm [ default ]
+      - 2: Munkres-Kuhn (Hungarian) algorithm 
       - 3: experimental stochastic algorithm, the original version with elimination (known problems)
       - 32: experimental stochastic algorithms with all permutations (too expensive)
       - 33: the improved stochastic algorithm with good scaling and performance, on par with the mincost
@@ -389,6 +390,7 @@ class dyn_control_params{
       - 5: SHXF of Min
       - 6: MQCXF
       - 7: DISH, rev2023
+      - 8: diabatic IDA, experimental
 
   */
   double decoherence_algo;
@@ -464,8 +466,13 @@ class dyn_control_params{
     only used with decoherence_algo == 1
 
       - 0: ID-S
-      - 1: ID-A [default]
+      - 1: ID-A [default] - if the proposed hop is not successful, we project back to the initial state
+                            if the proposed hop is accepted - we project onto that state
       - 2: ID-C - consistent ID - an experimental algorithm
+      - 3: ID-A, new: if the proposed hop is not successful, we project out the proposed states
+                      if the proposed hop is accepted - we project onto that state
+      - 4: ID-F, new: if the proposed hop is not successful, we project out the proposed states
+                      but we don't do anything if the hop is successful
   */
   int instantaneous_decoherence_variant;
 
@@ -496,31 +503,71 @@ class dyn_control_params{
 
 
   /**
-  The width of frozen Gaussian for the decoherence from SHXF & MQCXF
-  [ default : 0.1 ]
+  MATRIX(ndof, 1) of (initial) wave packet widths for the decoherence from SHXF & MQCXF
+  [ default : NULL ]
   */
-  double wp_width;
+  MATRIX* wp_width;
+
+
+  /**
+  MATRIX(ndof, 1) of wave packet velocities for the decoherence from SHXF & MQCXF
+  This value is applied when use_td_width = 1
+  [ default : NULL ]
+  */
+  MATRIX* wp_v;
 
 
   /**
   The criterion whether the electronic state is in a coherence
   [ default : 0.01 ]
   */
   double coherence_threshold;
-
+
+
+  /**
+  The masking parameter for computing nabla phase vectors in the MQCXF
+  [ default : 0.0001 ]
+  */
+  double e_mask;
+
+
   /**
   Whether to use the decoherence force in MQCXF
   The corresponding electronic propagation is adjusted for the energy conservation
   [ default : 0 ]
   */
   int use_xf_force;
 
+
   /**
   Whether to project out the density on an auxiliary trajectory when its motion is classically forbidden
   [ default : 0 ]
   */
   int project_out_aux;
 
+
+  /**
+  Turning-point algorithm for auxiliary trajectories
+
+  Options:
+      - 0: no treatment of a turning point
+      - 1: collapse to the active state [default]
+      - 2: fix auxiliary positions of adiabatic states except for the active state
+      - 3: keep auxiliary momenta of adiabatic states except for the active state
+  */
+  int tp_algo;
+
+
+  /**
+  Whether to use the td Gaussian width for the nuclear wave packet approximation
+  This option can be considered when it comes to unbounded systems.
+  This approximation is based on a nuclear wave packet on a free surface:
+    \sigma_x(t)=\sqrt[\sigma_x(0)^2 + (wp_v * t)^2]
+  [ default : 0 ]
+  */
+  int use_td_width;
+
+
   /**
   A flag for NBRA calculations. Since in NBRA, the Hamiltonian is the same for all the trajectories
   we can only compute the Hamiltonian related properties once for one trajectory and increase the speed of calculations.
@@ -695,6 +742,14 @@ class dyn_control_params{
   */
   int electronic_integrator;
 
+
+  /**
+    Whether transform the amplitudes by the T transformation matrix
+         0              - do not transform by the T matrix (naive, but potentially correct approach) 
+         1              - do transform it (as in LD, but maybe not needed if we directly transform basis)
+  */
+  int ampl_transformation_method;
+
 
   /**
     If set to True (1), we will force the reprojection matrix T_new to be the identity matrix. This effectively

src/dyn/dyn_decoherence.h

@@ -54,6 +54,9 @@ void instantaneous_decoherence(CMATRIX& Coeff,
    vector<int>& accepted_states, vector<int>& proposed_states, vector<int>& initial_states,
    int instantaneous_decoherence_variant, int collapse_option);
 
+void instantaneous_decoherence_dia(CMATRIX& Coeff, nHamiltonian& ham,
+   vector<int>& accepted_states, vector<int>& proposed_states, vector<int>& initial_states,
+   int instantaneous_decoherence_variant, int collapse_option);
 
 CMATRIX afssh_dzdt(CMATRIX& dz, CMATRIX& Hvib, CMATRIX& F, CMATRIX& C, double mass, int act_state);
 void integrate_afssh_moments(CMATRIX& dR, CMATRIX& dP, CMATRIX& Hvib, CMATRIX& F, CMATRIX& C, double mass, int act_state, double dt, int nsteps);
@@ -72,13 +75,15 @@ CMATRIX mfsd(MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector<MATRIX>&
 
 // Independent-trajectory XF methods
 void xf_hop_reset(dyn_variables& dyn_var, vector<int>& accepted_states, vector<int>& initial_states);
+void update_ham_xf(dyn_variables& dyn_var);
 
 void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn_control_params& prms); // For SHXF
 void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn_control_params& prms); // For MQCXF
 
 // XF propagation
+void rotate_nab_phase(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms);
 void update_forces_xf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev);
-void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms, int do_rotation);
+void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms);
 void XF_correction(CMATRIX& Ham, dyn_variables& dyn_var, CMATRIX& C, double wp_width, CMATRIX& T, int traj);
 
 ///================  In dyn_decoherence_time.cpp  ===================================