I need to calculate the average energy of a polymer chain confined inside a cylindrical pore as a part of the Thermodynamic Integration. Since accuracy is important for this purpose, I do not want the particles to ‘penetrate’ through the wall due to finite time steps. The way the velocity Verlet algorithm works allows particles to penetrate the wall and then immediately get repelled out in the subsequent time step. This artificially raises the average potential energy of the system rendering the thermodynamic integration highly inaccurate, numerically.
Are there any ways around this, other than using very small time steps? I am thinking of using a predictor-corrector MD algorithm to integrate the equations of motion, but this requires modifying LAMMP’s source code. How about doing MC simulation with LAMMPS? My polymer chain has a complex topology with side branches and rings attached to the backbone. I would greatly appreciate any suggestions.
Thank you.
All the “soft” walls (i.e. using harmonic, lj, morse or tabulated potential) do not suffer from this, since the wall’s repulsion starts way before the actual wall location and for as long as the displacement is small enough to not go “through” that barrier in a single MD step, the particles will convert kinetic energy into potential energy and then back to kinetic energy once they turn around.
The situation you describe would only happen for a very badly written MD code.
Thank you for your quick response. A reflective wall is not useful for my case because I NEED to have ‘walls’ with varying ‘softness’. Maybe I should not call them walls.
Imagine a cylindrical region (characterized by r<r0, axis parallel to x) inside my cubic simulation box. There is an external force field, characterized by its magnitude f0, only outside the cylindrical region. The direction of the force always points inward, toward the cylinder axis, so the effect of the force field is to push the particles into the cylindrical force-free region. If f0 is small, particles will explore both regions, when f0 is too large, the particles will be effectively confined inside the cylindrical region. In my simulations, I ‘confine’ particles inside a cylindrical region with varying degrees of confinement, corresponding to the lambda parameter of the thermodynamic integration method. At large f0, since all the particles should stay inside the cylindrical region, the average potential energy due to the external force should approach zero. But I do not see that because of the finite time step and my discontinuous (step function) force field. Making the force field smooth (eg tanh[kr]) is not a solution because the particular thermodynamic integration method we are implementing needs the force field to be step-like. Something like a predictor-corrector would prevent particles from going outside the cylindrical region when f0 is large, would you advise using that?
I know the problem is inherent to my non-smooth force field and the discrete nature of any numerical computation, and I am trying to mitigate it. A high degree of numerical accuracy is required for thermodynamic integration. Thank you so much again for your time and advice.
If you want to use a different time stepping algorithm than velocity Verlet, you have to use a different MD code.
If you want such an unusually “boundary potential”, you will have to write your own fix and then apply the same processing as fix wall/reflect does to mitigate the effect of the finite timestep:
When the particle has moved beyond the boundary surface, you can reconstruct the point where this happened from the velocity since it moved linearly at the position update. You can then determine for which fraction of the timestep the particle was inside the bounding potential and then add the bounding potential for that fraction of the time. This fix will have to run during the post-force step of time integration, so the force is added for the following two half velocity updates in final-integrate and initial-integrate of the next step.