Merge pull request #272 from suarez-va/devel

Fixed KC-RPMD compilation issues. new HEOM features documented. HEOM gen_heirarchy() now very efficient, matsubara frequencies switched to complex typing.

Author: alexvakimov

src/dyn/Dynamics.cpp

@@ -1227,7 +1227,6 @@ void propagate_electronic_kcrpmd(dyn_variables &dyn_var, nHamiltonian &Ham,
   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;
@@ -1245,7 +1244,7 @@ void propagate_electronic_kcrpmd(dyn_variables &dyn_var, nHamiltonian &Ham,
       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];
+                                                           a, b, c))[0];
       p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
     } else {
       sigma = sqrt(2.0 * gamKP / (beta * m_aux_var[0]));
@@ -1260,7 +1259,7 @@ void propagate_electronic_kcrpmd(dyn_variables &dyn_var, nHamiltonian &Ham,
       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];
+                                                           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 -
@@ -1273,7 +1272,7 @@ void propagate_electronic_kcrpmd(dyn_variables &dyn_var, nHamiltonian &Ham,
       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];
+                                                           a, b, c))[0];
       p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
     } else {
       sigma = sqrt(2.0 * gam / (beta * m_aux_var[0]));
@@ -1288,7 +1287,7 @@ void propagate_electronic_kcrpmd(dyn_variables &dyn_var, nHamiltonian &Ham,
       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];
+                                                           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 -

src/dyn/Energy_and_Forces.cpp

@@ -287,8 +287,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;
@@ -437,8 +436,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.h

@@ -531,8 +531,6 @@ class dyn_control_params{
   */
   double kcrpmd_c;
 
-  double kcrpmd_d;
-
   /**
     KC-RPMD Langevin frictional coefficient within donor-acceptor basin 
     Default: 0.0 a.u. (inverse time)