Implementation of MTK + Langevin propagation

[RudolfWeeber]

Goal

• add a propagation with barostat based on MTK equation and Langevin thermostat.
• based on implementation notes by Cristiano.
• replaces the present npt implementation

steps

• the npt state data structure should contain:
  • instantaneous pressure (p_at in notes),
  • degree of freedom of barostat (p_eps in notes),
  • Degrees of freedom (with center_of mass conservation) 3n -3
  • maybe sum of kinetic energy of all particles (K2 in c++ code in notes)
• the npt parameter data struct shoudl contain
  • external pressure
  • piston mass
  • langevin friction coeff (gamma) and thermal energy (kT)
• in src/core/integrators/velocity_verlet_npt.* implement single particle functions for the propagation steps in the notes. Keep in mind that the implementation notes refer to momenta (p), whereas Espresso uses velocities v.
  • Utils::Vector3d velocity_mtk_half_step(double Nf, double p_eps, double W, double dt)
  • double kinetic_energy(Vector3d velocity, double mass)
  • double p_eps_half_step(p_eps, p_inst, p_ext,sum_kinetic_energy,dt)
  • velocity_force_half_step(velocity, force, mass, dt)
  • postion_mtk_hafl_step(position, p_eps, W, dt)
  • volume_full_step(volume, p_eps, W, dt)
  • position_velocity_full_step(position, velocity, dt)
• Following the c++ sample code, deploy above propagation functions in the correct order. Espresso uses two enty points (step 1 before force calculation, step 2 after force calculation). For npt, these are velocity_verlet_npt_step_1 and 2. This functions will loop over all particles (seeral times) and execute the individual steps introduced above. Note that some steps can be executed in the same loop, e.g. velocity_mtk_hafl_step and the kinetic energy summation. Furthermore, some prefactors should probably be pre-calculated and passed into the respective steps, e.g., exp(0.5*dt*p_eps/W)

[RudolfWeeber]

@stekajack would you be willing to check whether that makes sense to you?

[stekajack]

@stekajack would you be willing to check whether that makes sense to you?

Based on a cursory reading, the descriptions align with my understanding of the problem and what I anticipated the implementation would look like. The implementation is based on the work of Cristiano de Michele (https://doi.org/10.1080/00268976.2024.2364834). This structure should make it straightforward to couple the MTK equations with another thermostat (i.e. the Noosé-Hoover thermostat, also discussed in the aforementioned publication )