Maintaining positional uncertainty for changing temperatures

Hello Dear LAMMPS users and developers!

I am using the LAMMPS installation from 29 September 2021 on a Linux machine.

In addition to a standard eam potential (e.g. Al99.eam.alloy), I am interested in applying a spring force unique to each atom that is proportional both to initial atomic position and to the temperature (i.e. Fspring = -K*(r - r0) , where K is proportional to temperature T.)

I would like to use this spring force to impose a restriction on a subset of atoms so that they do not deviate from their initial positions too greatly, with some specified anisotropic standard deviation (e.g. 2 nm allowed deviation in the xz plane, 10 nm allowed deviation in the y direction.)

The most sensible fix for this seems to be fix spring/self.

The issue: If I change the temperature of the system over the duration of the simulation, I need to be able to adjust the spring constant such that the standard deviation remains constant (i.e. increase in T requires increase in K_spring to maintain std).
fix spring/self uses the position at the time of the call as the initial position (r0) for calculating the force [F = -k * (r - r0)].
Therefore, updating the spring constants and re-calling fix spring/self causes the r0 to wander away from the true initial position (undesirable).

The question: Is there an alternative way to implement the described system that I have overlooked? Can the spring constant be dynamically updated without re-calling the fix? Do I need to expand the fix spring/self code s.t. I can read in the r0 values instead of writing the present configuration to a binary restart file (a bit daunting since I am new to LAMMPS)?

Any insight would be greatly appreciated,

Thank you for your thoughts,
Dylan

The most straightforward solution to your issue would be to modify fix spring/self to allow specifying K as an equal style (or even atom style) variable. Similar things exist for other fixes.

You could also see if you can implement your position restraint with fix store/state (to record r0) and then a bunch atom style variables for Fx, Fy, Fz for each particle (you need variables for r, dx, dy, and dz for that) and have K as an equal style variable which can be made dependent on the timestep with the β€œramp” function or similar. Then you could use fix addforce to apply those forces. This can be a bit tricky with periodic boundaries, so you probably also would need to use compute property/atom to get access to unwrapped positions.

1 Like

@dcmiley
FYI, I have just implemented support for entering the spring constant K as either an equal style or atom style variable, which allows to adjust the spring constant over time and through space.

1 Like

I greatly appreciate you implementing this feature - you surely saved me a week of headache and this will be a good example for learning how to make modifications for myself in the future.