added kcrpmd langevin thermostat parameter for kinked-pair regime

Author: suarez-va

src/dyn/Dynamics.cpp

@@ -1126,12 +1126,13 @@ void propagate_electronic_kcrpmd(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_
   int num_el = prms.num_electronic_substeps;
   double dt = prms.dt / num_el;
   double beta = (hartree/boltzmann) / prms.Temperature; 
-  double eta = prms.kcrpmd_eta; 
-  double gam = prms.kcrpmd_gamma;
   double a = prms.kcrpmd_a; 
   double b = prms.kcrpmd_b; 
   double c = prms.kcrpmd_c; 
   double d = prms.kcrpmd_d; 
+  double eta = prms.kcrpmd_eta; 
+  double gam = prms.kcrpmd_gamma;
+  double gamKP = prms.kcrpmd_gammaKP;
   double sigma;
   double xi;
   double theta;
@@ -1141,26 +1142,51 @@ void propagate_electronic_kcrpmd(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_
   vector<double>& p_aux_var = dyn_var.p_aux_var;
   vector<double>& f_aux_var = dyn_var.f_aux_var;
 
-  if(gam == 0.0){
-    p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
-    y_aux_var[0] += p_aux_var[0] / m_aux_var[0] * dt;
-    f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c, d))[0];
-    p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
+  if(y_aux_var[0] >= -0.5 && y_aux_var[0] <= 0.5){
+    if(gamKP == 0.0){
+      p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
+      y_aux_var[0] += p_aux_var[0] / m_aux_var[0] * dt;
+      f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c, d))[0];
+      p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
+    }
+    else{
+      sigma = sqrt(2.0 * gamKP / (beta * m_aux_var[0]));
+      xi = rnd.normal();
+      theta = rnd.normal();
+      p_aux_var[0] += 0.5 * xi * sigma * m_aux_var[0] * sqrt(dt)
+                    + 0.5 * (f_aux_var[0] - gamKP * p_aux_var[0]) * dt
+	  	            - 0.25 * (xi / 2.0 + theta / sqrt(3.0)) * gamKP * sigma * m_aux_var[0] * sqrt(pow(dt,3))
+	  	            - 0.125 * gamKP * (f_aux_var[0] - gamKP * p_aux_var[0]) * pow(dt,2);
+      y_aux_var[0] += p_aux_var[0] / m_aux_var[0] * dt + 0.5 * theta / sqrt(3.0) * sigma * sqrt(pow(dt,3));
+      f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c, d))[0];
+      p_aux_var[0] += 0.5 * xi * sigma * m_aux_var[0] * sqrt(dt)
+	                 + 0.5 * (f_aux_var[0] - gamKP * p_aux_var[0]) * dt
+		             - 0.25 * (xi / 2.0 + theta / sqrt(3.0)) * gamKP * sigma * m_aux_var[0] * sqrt(pow(dt,3))
+		             - 0.125 * gamKP * (f_aux_var[0] - gamKP * p_aux_var[0]) * pow(dt,2);
+    }
   }
   else{
-    sigma = sqrt(2.0 * gam / (beta * m_aux_var[0]));
-    xi = rnd.normal();
-    theta = rnd.normal();
-    p_aux_var[0] += 0.5 * xi * sigma * m_aux_var[0] * sqrt(dt)
-	          + 0.5 * (f_aux_var[0] - gam * p_aux_var[0]) * dt
-		  - 0.25 * (xi / 2.0 + theta / sqrt(3.0)) * gam * sigma * m_aux_var[0] * sqrt(pow(dt,3))
-		  - 0.125 * gam * (f_aux_var[0] - gam * p_aux_var[0]) * pow(dt,2);
-    y_aux_var[0] += p_aux_var[0] / m_aux_var[0] * dt + 0.5 * theta / sqrt(3.0) * sigma * sqrt(pow(dt,3));
-    f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c, d))[0];
-    p_aux_var[0] += 0.5 * xi * sigma * m_aux_var[0] * sqrt(dt)
-	          + 0.5 * (f_aux_var[0] - gam * p_aux_var[0]) * dt
-		  - 0.25 * (xi / 2.0 + theta / sqrt(3.0)) * gam * sigma * m_aux_var[0] * sqrt(pow(dt,3))
-		  - 0.125 * gam * (f_aux_var[0] - gam * p_aux_var[0]) * pow(dt,2);
+    if(gam == 0.0){
+      p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
+      y_aux_var[0] += p_aux_var[0] / m_aux_var[0] * dt;
+      f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c, d))[0];
+      p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
+    }
+    else{
+      sigma = sqrt(2.0 * gam / (beta * m_aux_var[0]));
+      xi = rnd.normal();
+      theta = rnd.normal();
+      p_aux_var[0] += 0.5 * xi * sigma * m_aux_var[0] * sqrt(dt)
+                    + 0.5 * (f_aux_var[0] - gam * p_aux_var[0]) * dt
+	  	            - 0.25 * (xi / 2.0 + theta / sqrt(3.0)) * gam * sigma * m_aux_var[0] * sqrt(pow(dt,3))
+	  	            - 0.125 * gam * (f_aux_var[0] - gam * p_aux_var[0]) * pow(dt,2);
+      y_aux_var[0] += p_aux_var[0] / m_aux_var[0] * dt + 0.5 * theta / sqrt(3.0) * sigma * sqrt(pow(dt,3));
+      f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c, d))[0];
+      p_aux_var[0] += 0.5 * xi * sigma * m_aux_var[0] * sqrt(dt)
+	                 + 0.5 * (f_aux_var[0] - gam * p_aux_var[0]) * dt
+		             - 0.25 * (xi / 2.0 + theta / sqrt(3.0)) * gam * sigma * m_aux_var[0] * sqrt(pow(dt,3))
+		             - 0.125 * gam * (f_aux_var[0] - gam * p_aux_var[0]) * pow(dt,2);
+    }
   }
 }
 

src/dyn/dyn_control_params.cpp

@@ -92,6 +92,7 @@ dyn_control_params::dyn_control_params(){
   kcrpmd_c = 0.5;
   kcrpmd_d = 3.0;
   kcrpmd_gamma = 0.0;
+  kcrpmd_gammaKP = 0.0;
 
   ///================= Decoherence options =========================================
   decoherence_algo = -1; 
@@ -211,6 +212,7 @@ dyn_control_params::dyn_control_params(const dyn_control_params& x){
   kcrpmd_c = x.kcrpmd_c;
   kcrpmd_d = x.kcrpmd_d;
   kcrpmd_gamma = x.kcrpmd_gamma;
+  kcrpmd_gammaKP = x.kcrpmd_gammaKP;
 
   ///================= Decoherence options =========================================
   decoherence_algo = x.decoherence_algo; 
@@ -422,6 +424,7 @@ void dyn_control_params::set_parameters(bp::dict params){
     else if(key=="kcrpmd_c"){ kcrpmd_c = bp::extract<double>(params.values()[i]);  }
     else if(key=="kcrpmd_d"){ kcrpmd_d = bp::extract<double>(params.values()[i]);  }
     else if(key=="kcrpmd_gamma"){ kcrpmd_gamma = bp::extract<double>(params.values()[i]);  }
+    else if(key=="kcrpmd_gammaKP"){ kcrpmd_gammaKP = bp::extract<double>(params.values()[i]);  }
 
     ///================= Decoherence options =========================================
     else if(key=="decoherence_algo"){ decoherence_algo = bp::extract<int>(params.values()[i]); }

src/dyn/dyn_control_params.h

@@ -538,11 +538,17 @@ class dyn_control_params{
   double kcrpmd_d;
 
   /**
-    KC-RPMD Langevin frictional coefficient 
+    KC-RPMD Langevin frictional coefficient within donor-acceptor basin 
     Default: 0.0
   */
   double kcrpmd_gamma;
 
+  /**
+    KC-RPMD Langevin frictional coefficient within kinked-pair regime 
+    Default: 0.0
+  */
+  double kcrpmd_gammaKP;
+
 
   ///===============================================================================
   ///================= Decoherence options =========================================