Merge branch 'devel' of https://github.com/Quantum-Dynamics-Hub/libra-code into devel

Author: alexvakimov

src/dyn/Dynamics.cpp

@@ -1129,7 +1129,6 @@ void propagate_electronic_kcrpmd(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_
   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;
@@ -1146,7 +1145,7 @@ void propagate_electronic_kcrpmd(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_
     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];
+      f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c))[0];
       p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
     }
     else{
@@ -1158,7 +1157,7 @@ void propagate_electronic_kcrpmd(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_
 	  	            - 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];
+      f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c))[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))
@@ -1169,7 +1168,7 @@ void propagate_electronic_kcrpmd(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_
     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];
+      f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c))[0];
       p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
     }
     else{
@@ -1181,7 +1180,7 @@ void propagate_electronic_kcrpmd(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_
 	  	            - 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];
+      f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c))[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))

src/dyn/Energy_and_Forces.cpp

@@ -288,8 +288,7 @@ vector<double> potential_energies(dyn_control_params& prms, dyn_variables& dyn_v
       double a = prms.kcrpmd_a; 
       double b = prms.kcrpmd_b; 
       double c = prms.kcrpmd_c; 
-      double d = prms.kcrpmd_d; 
-      Veff = ham.kcrpmd_effective_potential(dyn_vars.y_aux_var, *dyn_vars.q, *dyn_vars.iM, beta, eta, a, b, c, d);
+      Veff = ham.kcrpmd_effective_potential(dyn_vars.y_aux_var, *dyn_vars.q, *dyn_vars.iM, beta, eta, a, b, c);
       for(itraj=0; itraj<ntraj; itraj++){
         id[1] = itraj;
         res[itraj] = Veff;
@@ -438,8 +437,7 @@ void update_forces(dyn_control_params& prms, dyn_variables& dyn_vars, nHamiltoni
      double a = prms.kcrpmd_a; 
      double b = prms.kcrpmd_b; 
      double c = prms.kcrpmd_c; 
-     double d = prms.kcrpmd_d; 
-     *dyn_vars.f = ham.kcrpmd_effective_force(dyn_vars.y_aux_var, *dyn_vars.q, *dyn_vars.iM, beta, eta, a, b, c, d);
+     *dyn_vars.f = ham.kcrpmd_effective_force(dyn_vars.y_aux_var, *dyn_vars.q, *dyn_vars.iM, beta, eta, a, b, c);
     }
 
     // Adiabatic 

src/dyn/dyn_control_params.cpp

@@ -89,8 +89,7 @@ dyn_control_params::dyn_control_params(){
   kcrpmd_eta = 6.28;
   kcrpmd_a = 0.1;
   kcrpmd_b = 1000.;
-  kcrpmd_c = 0.5;
-  kcrpmd_d = 3.0;
+  kcrpmd_c = 1.0;
   kcrpmd_gamma = 0.0;
   kcrpmd_gammaKP = 0.0;
 
@@ -210,7 +209,6 @@ dyn_control_params::dyn_control_params(const dyn_control_params& x){
   kcrpmd_a = x.kcrpmd_a;
   kcrpmd_b = x.kcrpmd_b;
   kcrpmd_c = x.kcrpmd_c;
-  kcrpmd_d = x.kcrpmd_d;
   kcrpmd_gamma = x.kcrpmd_gamma;
   kcrpmd_gammaKP = x.kcrpmd_gammaKP;
 
@@ -422,7 +420,6 @@ void dyn_control_params::set_parameters(bp::dict params){
     else if(key=="kcrpmd_a"){ kcrpmd_a = bp::extract<double>(params.values()[i]);  }
     else if(key=="kcrpmd_b"){ kcrpmd_b = bp::extract<double>(params.values()[i]);  }
     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]);  }
 

src/dyn/dyn_control_params.h

@@ -509,43 +509,37 @@ class dyn_control_params{
 
   /**
     KC-RPMD free energy parameter eta
-    Default: 6.28
+    Default: 6.28 (dimensionless)
   */
   double kcrpmd_eta;
 
   /**
     KC-RPMD kinetic constraint parameter a
-    Default: 0.1
+    Default: 0.1 (dimensionless)
   */
   double kcrpmd_a;
 
   /**
     KC-RPMD heavy-side parameter b
-    Default: 1000.
+    Default: 1000. (dimensionless)
   */
   double kcrpmd_b;
 
   /**
     KC-RPMD kinetic constraint switching parameter c
-    Default: 0.5
+    Default: 1.0 (dimensionless)
   */
   double kcrpmd_c;
 
-  /**
-    KC-RPMD free energy switching parameter d
-    Default: 3.0
-  */
-  double kcrpmd_d;
-
   /**
     KC-RPMD Langevin frictional coefficient within donor-acceptor basin 
-    Default: 0.0
+    Default: 0.0 a.u. (inverse time)
   */
   double kcrpmd_gamma;
 
   /**
     KC-RPMD Langevin frictional coefficient within kinked-pair regime 
-    Default: 0.0
+    Default: 0.0 a.u. (inverse time)
   */
   double kcrpmd_gammaKP;