lammps
Reversible integration of SLLOD
Summary
This PR updates fix nvt/sllod, fix deform, and compute temp/deform (and some other associated code) for stable integration of the equations of motion in a way which is reversible and energy-conserving (once work done to drive the flow is accounted for). With these changes, energy dissipation is now in good agreement with the expected analytical value ($-V\mathbf{P}^T:\nabla\mathbf{u}$), and it is now possible to do this:
Logo input script
# simulation of LAMMPS logo under SLLOD
units lj
atom_style atomic
dimension 2
# create 2d lattice, low density so diffuses
lattice sq2 0.5 origin 0.25 0.25 0.0
region box block 0 31 0 7 -0.5 0.5
create_box 2 box
create_atoms 1 box
# LAMMPS letters via regions, convert to type 2 atoms
region L1 block 2 3 1 6 -0.5 0.5
region L2 block 2 5 1 2 -0.5 0.5
region L union 2 L1 L2
region A1 block 6 7 1 6 -0.5 0.5
region A2 block 8 9 1 6 -0.5 0.5
region A3 block 6 9 3 4 -0.5 0.5
region A4 block 6 9 5 6 -0.5 0.5
region A union 4 A1 A2 A3 A4
region 1M1 block 10 11 1 6 -0.5 0.5
region 1M2 block 14 15 1 6 -0.5 0.5
region 1M3 block 10 15 5 6 -0.5 0.5
region 1M4 block 12 13 3 6 -0.5 0.5
region 1M union 4 1M1 1M2 1M3 1M4
region 2M1 block 16 17 1 6 -0.5 0.5
region 2M2 block 20 21 1 6 -0.5 0.5
region 2M3 block 16 21 5 6 -0.5 0.5
region 2M4 block 18 19 3 6 -0.5 0.5
region 2M union 4 2M1 2M2 2M3 2M4
region P1 block 22 23 1 6 -0.5 0.5
region P2 block 24 25 3 6 -0.5 0.5
region P3 block 22 25 3 4 -0.5 0.5
region P4 block 22 25 5 6 -0.5 0.5
region P union 4 P1 P2 P3 P4
region S1 block 26 29 5 6 -0.5 0.5
region S2 block 26 27 3 6 -0.5 0.5
region S3 block 26 29 3 4 -0.5 0.5
region S4 block 28 29 1 4 -0.5 0.5
region S5 block 26 29 1 2 -0.5 0.5
region S union 5 S1 S2 S3 S4 S5
region LAMMPS union 6 L A 1M 2M P S
set region LAMMPS type 2
# LJ and other system parameters
mass * 1.0
pair_style lj/cut 2.5
pair_coeff * * 1.0 1.0 2.5
timestep 0.005
neighbor 0.3 bin
neigh_modify delay 0 every 1 check yes
change_box all triclinic
velocity all create 2.0 87287 loop geom
# shear with compression in x, expansion in y
# for volume preservation, can set yy = -xx
variable xy index 0.7
variable xx index -0.05
variable yy index 0.2
# run long enough to get a box flip, but short enough
# to stay reversible with floating point error
variable runtime index 800
thermo 100
dump 1 all atom 20 logo.lammpstrj
# run without integrator to "pause" the visualization
run 100
# dissolve the lattice
fix def all deform 0 xy erate ${xy} y trate ${yy} x trate ${xx} remap none couple/erate yes
fix 1 all nvt/sllod temp 1 1 0.1 peculiar yes
compute off all temp/partial 0 0 0
compute t all temp/deform temp off
fix_modify 1 temp off
fix 2 all enforce2d
run ${runtime}
# run without integrator to "pause" the visualization
unfix 1
unfix 2
unfix def
run 200
# reverse the peculiar-frame velocities and continue run
# logo and lattice should reassemble within round-off errors
variable vxflip atom -vx
variable vyflip atom -vy
velocity all set v_vxflip v_vyflip NULL
# do reverse step with lab-frame velocities, because we can
fix def all deform 0 xy erate $(-v_xy) y trate $(-v_yy) x trate $(-v_xx) remap v couple/erate yes
fix 1 all nvt/sllod temp 1 1 0.1 kick yes
fix_modify 1 temp t
fix 2 all enforce2d
run ${runtime}
# run without integrator to "pause" the visualization
unfix 1
unfix 2
run 100
The initial motivation for these changes was some discrepancies/issues we noticed when looking at the transient response and trying to do response theory calculations, but in the end we have found that even viscosity calculations at steady state were slightly affected by the issues fixed in this PR (viscosity was slightly overestimated, becoming significant at high shear rates or under mixed flows). I'm currently in the final stages of preparing a paper with the details, which I'm happy to share privately if you want to see it before we put it on arXiv. I've already sent it to Pieter in 't Veld for his comments, and he hasn't flagged any concerns.
The main changes are:
- Account for the change in streaming velocity caused by updating particle positions. Without this, the thermal velocity is implicitly changed.
- Add an option to update the unit cell with
fix deformduringpost_integraterather thanend_of_stepto avoid an off-by-one error on the first step (SLLOD technically requires that the box be updated at the same time as particle positions to avoid an additional force calculation). - Apply the position/velocity updates due to SLLOD in a reversible way (important for mixed flows).
- Account for the effects of elongation (x/y/z TRATE) on shear (xy/xz/yz ERATE), as well as the effect of xy shear on xz tilt when there is yz tilt or shear (including handling box flips). Without this, mixed flows have some major issues, including the flow tensor (as calculated by
fix nvt/sllod) not staying constant. - Fix a bug where velocities were not re-mapped after a box flip or when unwrapping coordinates for calculations like angular momentum.
- Set
h_ratefrom the start for TRATE deform style so that the correct flow rate tensor is used byfix nvt/sllodin the first time step.
See the implementation details section below for additional notes on some of these changes (numbers match up).
Still to-do:
- Updates to documentation - I've left this for now until the user interface is finalised, but see below for an example input script.
- Accelerator styles - I've focused on the bare styles for now, but it should be easy enough to propagate the changes to the accelerator variants once any issues have been sorted out.
- ~~Fix merge conflicts with
develop- I based this on the29Aug2024_update4release to be able to reference it in the paper.~~
Example Input Script
# Examples:
# Simple shear:
# gamma_xy = 0.1
# others = 0
# Mixed shear:
# gamma_xy = 0.04
# gamma_xz = 0.1
# gamma_yz = 0.02
# others = 0
# Mixed shear v2 (only supported by this patch unless run for a shorter time):
# gamma_xy = 0.04
# gamma_xz = 0.06
# gamma_yz = 0.1
# others = 0
# Simple planar elongation:
# eps_xx = 0.05
# eps_yy = -0.05
# others = 0
# Mixed shear + elongation (only supported by this patch unless run for a shorter time):
# gamma_xy = 0.04
# gamma_xz = 0.06
# gamma_yz = 0.01
# eps_xx = 0.025
# eps_yy = -0.05
# eps_zz = 0.025 (calculated for volume preservation)
variable gamma_xy index 0.1
variable gamma_xz index 0.06
variable gamma_yz index 0.1
variable eps_xx index 0.0025
variable eps_yy index -0.005
variable eps_zz index $(-v_eps_xx-v_eps_yy)
variable tcouple index 0.1
variable wx equal omega(all,x)
variable wy equal omega(all,y)
variable wz equal omega(all,z)
# set no to turn off thermostat and have atoms with zero thermal velocity (useful for debugging)
variable thermostat index yes
# set yes to store velocity in peculiar frame (i.e. velocity relative to streaming velocity)
variable peculiar index no
units lj
atom_style atomic
lattice sc 0.7
# Check for box origin dependence.
region box block 0 5 0 20 0 10
# region box block -2.5 2.5 -10 10 -5 5
create_box 1 box
create_atoms 1 box
mass * 1.0
pair_style lj/cut 2.5
pair_coeff * * 1.0 1.0 1.0
pair_modify shift yes
change_box all triclinic
timestep 0.0005
neighbor 0.3 bin
neigh_modify delay 0 every 1
delete_atoms random fraction 0.5 yes all NULL 481913
#
if "${thermostat} == yes" then &
"velocity all create 1.0 902384 mom yes" &
else &
"velocity all set 0 0 0"
# Get current energy for later
thermo_modify norm no
run 0
variable Einitial index $(econserve)
if "${peculiar} == yes" then "variable remap index none" &
else "variable remap index v"
# Modified version of fix deform uses N = 0 to indicate
# post_force() box updates instead of end_of_step()
fix def all deform 0 &
x trate ${eps_xx} &
y trate ${eps_yy} &
z trate ${eps_zz} &
xy erate ${gamma_xy} &
xz erate ${gamma_xz} &
yz erate ${gamma_yz} &
remap ${remap} couple/erate yes
# Using default `kick yes` flag.
# Internal temperature compute is either temp/deform or temp depending on ${peculiar}
fix sllod all nvt/sllod temp 1.0 1.0 ${tcouple} tchain 1 peculiar ${peculiar}
if "${thermostat} == yes" then "jump SELF END_THERMOSTAT"
# Disable thermostat
compute off all temp/partial 0 0 0
if "${peculiar} == yes" then "fix_modify sllod temp off" "jump SELF ALIAS"
compute tdeform all temp/deform temp off
fix_modify sllod temp tdeform
label ALIAS
# Alias sllod_temp for pressure calculation
if "${peculiar} == yes" then &
"compute sllod_temp all temp" &
else &
"compute sllod_temp all temp/deform"
label END_THERMOSTAT
compute press all pressure sllod_temp
# dump 1 all custom 100 data.lammpstrj id mass x y z vx vy vz
# Calculate expected energy dissipation rate as dE/dt = -VP^T:∇u - 2Kα
# -2Kα term handled by econserve.
# LJ units so unit conversion factor is 1.
variable Edot_expected equal -vol*(${eps_xx}*c_press[1]+${eps_yy}*c_press[2]+${eps_zz}*c_press[3]+${gamma_xy}*c_press[4]+${gamma_xz}*c_press[5]+${gamma_yz}*c_press[6])
# Simple integration of Edot_expected
fix running all ave/time 1 1 1 v_Edot_expected ave running
variable Eexpected equal v_Einitial+f_running*step*dt
variable Ediff equal v_Eexpected-econserve
variable Ediff_percent equal abs(v_Ediff*100/(v_Eexpected+(v_Eexpected==0)*econserve+(v_Ediff==0)))
# Expect Ediff and Ediff_percent to stay near 0.
thermo_style custom step temp pxx pyy pzz pxy pxz pyz vol pe ke v_wx v_wy v_wz econserve v_Eexpected v_Ediff v_Ediff_percent
thermo_modify temp sllod_temp press press norm no
thermo 500
# run_style respa 2 4
# Test that kick is only applied once
run 100 post no
# Could manually force kick to apply again with
# fix_modify sllod kick yes
run 1000000
Related Issue(s)
Author(s)
Stephen Sanderson - The University of Queensland
Licensing
By submitting this pull request, I agree, that my contribution will be included in LAMMPS and redistributed under either the GNU General Public License version 2 (GPL v2) or the GNU Lesser General Public License version 2.1 (LGPL v2.1).
Artificial Intelligence (AI) Tools Usage
By submitting this pull request, I confirm that I did NOT use any AI tools to generate all or parts of the code and modifications in this pull request.
Backward Compatibility
Since the new kick flag of fix nvt/sllod defaults to yes when using lab-frame velocity, scripts which previously applied the velocity profile manually would need to either stop doing that or set kick no.
Previous scripts will also generate a warning about fix deform being applied at the end of the step, and suggest setting nevery = 0 in fix deform to correct the issue.
With fix nvt/sllod now having compute temp/deform as a hard requirement, scripts which used some other velocity bias will need to adapt by setting that bias via the new temp argument of compute temp/deform.
Implementation Notes
- The cleanest way I found to handle velocity bias changing with position updates was to modify
compute temp/deformto have an additional method for applying the velocity bias based on current particle positions, rather than just restoring it from whatever was previously calculated. This meant that use oftemp/deformbysllodbecame a hard requirement, so to preserve the ability to use other biases (and reduce duplicate code),temp/deformnow uses an internal temperature compute (tempby default), which is only called at times when the bias due tofix deformhas been removed. This also enables new possibilities like turning off the thermostat in some/all directions withcompute temp/partial, and I think will work nicely with our plans for a molecular thermostat. The other benefit here is that we can add akick yes/noflag for applying the streaming velocity profile at the start of the first run, rather than having the user do so manually (which has some pitfalls to the uninitiated due to shear using the lower box corner as the origin, but elongation using the center of the box, calculated as if tilt factors were zero). - For now, I've made the toggle for applying box deformation in
post_integrateto be settingnevery = 0infix deform, since using it only really makes sense (at least to me) if deformation is being applied every time step. I realise that's not particularly intuitive as a user interface though, so I'm happy to switch to a flag plus some error checking or some other solution if that's preferred. I've also made sure this is compatible with respa. - For reversible integration, it's quite a bit simpler (and a little more computationally efficient) to store velocity in the peculiar (thermal) frame rather than the lab frame. To this end, I've added a
peculiar = yes/noflag tofix nvt/sllodto toggle on use of peculiar frame velocities. Using peculiar velocities would be problematic for some calculations like angular momentum which assume velocities are in the lab frame, but when those aren't needed I think it's worth having the option for the performance it gives (and it was withpeculiar yesthat I first got reversibility working properly - lab frame was more sensitive). - I've only implemented handling of mixed flows for the specific combination of x/y/z TRATE + xy/xz/yz ERATE for now, since my focus was on making SLLOD work correctly and it is those styles which give a constant $\nabla\mathbf{u}$. Some thought should probably be given to other deformation styles, but I'm not so familiar with their use cases and didn't want to inadvertently break something. I also haven't looked closely yet at the new
fix deform/pressurestyle to determine how/whether these changes should interact with it (perhaps @jtclemm could weigh in on this?). For handling long-running xy + yz deformation with a box flip, I've required xz to also be set to ERATE (even if it's ERATE 0) so that we can be given explicit control over it. The solution to the ODE for dealing with this is a bit annoying with (many) different cases depending on the structure of the flow tensor; for now I've used some simple equality checks to choose the correct solution since the rates being checked are user inputs and therefore unlikely to be "very close but not exactly the same", but a proper floating pointnearly_equalfunction would be a more robust solution. I've flagged places where I think those should be used, and if that's the way forward then my suggestion is to movenearly_equalfrom the UEF package (possibly with minor modifications for robustness) intomath_extraor similar.
Possible future expansion:
Beyond the scope of this PR (I think), and perhaps this should be a separate issue, but it would be nice to be able to support things like fix npt/sllod and long-running mixed flows. I think this is difficult to do with the current fix */sllod + fix deform setup since both fixes would want to be changing the unit cell, so to that end, I think it would be useful if sllod were to be separated from fix deform and instead control the unit cell itself and integrate it directly instead of trying to solve analytically for its shape (similar to how the UEF package does it). This certainly isn't the only solution though, and I'm suggesting this mainly as a starting point for conversation.
Post Submission Checklist
- [ ] The feature or features in this pull request is complete
- [ ] Licensing information is complete
- [ ] Corresponding author information is complete
- [ ] The source code follows the LAMMPS formatting guidelines
- [ ] Suitable new documentation files and/or updates to the existing docs are included
- [ ] The added/updated documentation is integrated and tested with the documentation build system
- [ ] The feature has been verified to work with the conventional build system
- [ ] The feature has been verified to work with the CMake based build system
- [ ] Suitable tests have been added to the unittest tree.
- [ ] A package specific README file has been included or updated
- [ ] One or more example input decks are included
Further Information, Files, and Links
~~After submitting, I see that this is also showing a bunch of changes from the diff between 29Aug2024_update4 and develop, so here is a compare link to see just my changes: https://github.com/lammps/lammps/compare/stable_29Aug2024_update4...ssande7:lammps:sllod-integration Please let me know if you want me to try and rebase before going further.~~
After submitting, I see that this is also showing a bunch of changes from the diff between 29Aug2024_update4 and develop, so here is a compare link to see just my changes:
I am sorry, but we cannot merge this pull request because it is not based on the develop branch but on the stable release branch. Those branches have diverged since we backport changes and make update releases with bugfixes.
Thus you will have to either rebase your changes to the develop branch or rather create a new feature branch from develop, create the diff for your changes only and then apply it to the feature branch, commit, and resubmit as a new pull request (assuming there were no conflicting changes to develop since).
@akohlmey no worries, I suspected that might be the case - I'll sort that out shortly.
Should be fixed now.
Thanks. Please let us know when you want people to take a closer look. I think @athomps is probably the first person that should check your changes.
I think it's ready now for someone to take a look. I've convinced myself of correctness, so it would be good to get feedback on any code structure/user interface/etc changes before I duplicate into accelerator styles and update the documentation.
@ssande7 This sounds like a significant improvement in the consistency of SLLOD in LAMMPS, particularly recovering the correct work production rate and also time reversibility. For people who want to be able to reproduce previous results, can you provide an option for that (including both psslod yes/no cases)? Will this change affect regression testing for examples/nemd/in.nemd and examples/VISCOSITY/in.nemd.2d? I would be interested in seeing the pre-arXiv manuscript.
Thanks @athomps , I think that should be possible. It wouldn't account for the bug fixes in domain, angmom, etc relating to velocity remapping, so it wouldn't give exactly the same results, but probably very close in large systems. I'll take a look and see what I can get working.
For fix deform, I think a flag to turn on/off handling of coupled directions makes sense. Fix nvt/sllod could also just use a flag, but it might get a little messy, so perhaps it would be simpler to copy the old version as a new fix and call it nvt/sllod/legacy or something like that. Would that fit ok with the LAMMPS way of doing things?
@ssande7 We prefer not to keep legacy versions hanging around. If there is a clean way to provide old functionality, then good, if not, let's not do it.
@athomps that makes sense - I'll see what I can do with just flags. I think at least the flag for fix deform makes sense to preserve the ability to mix x/y/z TRATE with xy/xz/yz ERATE without them coupling.
@athomps I've added some flags which I think should achieve close enough to the previous behaviour - I'll leave it up to you whether it's worth the added complication (which isn't that much) and added risk of people using it incorrectly.
For fix nvt/sllod, I've added an integrator flag which defaults to reversible but can be set to legacy. Switching between the two is fairly easy since they use different hooks during the time step, so error/warning handling was the most complicated part (and not that bad).
For fix deform, I've called the flag couple, and made it default to no since it's not compatible with most deformation styles and warns about those cases. I also added some warnings to fix nvt/sllod in cases where couple yes is needed.
Probably some bike shedding to do on the flag names if we go ahead with them, but let me know what you think. I've checked that couple no gives the same box shape as the old implementation (100 steps under 2-way shear + biaxial extension). I also compared integrator legacy + couple no with the old version, and for that the agreement wasn't exact (as expected) because of some of the other bug fixes, but it looked very close. Note, one change which isn't relevant to this short run but would make some difference in longer runs is around line 1055 of fix deform, where it now accounts for the xy tilt factor affecting the xz tilt when yz is remapped (perhaps this should also be gated behind the couple flag, but I'm not sure).
I've also edited the example input scripts in my earlier comment to use couple yes.
That's a cool demonstration video! I had two immediate questions to help me understand the changes:
- Are the changes to remapping in the domain class equivalent to the changes in this PR #4577?
- Does the coupling/direction keyword accomplish the same thing as the erate/rescale option in fix deform/pressure? If so, is that a more concise way to add broader support beyond trate? I still need to read calc_xz_correction() more closely to understand its effect. Relevant codeblock here: https://github.com/lammps/lammps/blob/abba7fcf5523b7854d0644bec8ced1d10200a4a7/src/EXTRA-FIX/fix_deform_pressure.cpp#L486-L496
Thanks @jtclemm
- Are the changes to remapping in the domain class equivalent to the changes in this PR https://github.com/lammps/lammps/pull/4577?
Yes, it looks like they are (and I think I like that version better), although it doesn't look like that PR passes through velocity in remap_all, which would be needed here, and it also doesn't appear to check deform_groupbit (which from memory is why I didn't just change remap). My changes also have the addition of an unmap function that handles velocity, which is needed for angmom, omega, etc.
- Does the coupling/direction keyword accomplish the same thing as the erate/rescale option in fix deform/pressure? If so, is that a more concise way to add broader support beyond trate? I still need to read calc_xz_correction() more closely to understand its effect.
In part yes, it does accomplish the same thing, but it also handles the effects of elongation in the shear direction (not just perpendicular), as well as the effects of xy shear on the xz component when there's yz shear or tilt (calc_xz_correction() is responsible for that part, taking into account how elongation will affect it all). All of those considerations are important for $\nabla \mathbf{u}$ to stay constant throughout the simulation. It essentially comes from solving $\dot{\mathbf{e}} = \mathbf{e}\cdot\nabla\mathbf{u}$ $\forall \mathbf{e}\in$ { $\mathbf{a}, \mathbf{b}, \mathbf{c}$ } (where $\mathbf{a}$, $\mathbf{b}$ and $\mathbf{c}$ are the lattice vectors), taking the limiting case of a triangular $\nabla\mathbf{u}$ (which is all we can handle with a reduced triclinic unit cell). I have the details of that math written up and could send you the draft paper privately if that would help (or we will hopefully get it onto arXiv in the next few weeks).
Sidenote - I just realised that keyword conflicts with the deform/pressure couple keyword, so I've renamed it to couple/erate (for now).
One other thought - if this is all too much complexity regarding how it may interact with other use cases, an alternative might be to add a style to fix deform for constant $\nabla \mathbf{u}$, which would guarantee that we have control over the full unit cell and just allow setting of each component of $\nabla\mathbf{u}$. In that case, rather than solving analytically for the unit cell, we could just integrate it directly, which would avoid the complexity of calc_xz_correction() and should be easier to extend into an npt/sllod.
but it also handles the effects of elongation in the shear direction (not just perpendicular)
If it is more careful than erate/rescale, there might be an argument to just have it replace that style.
In that case, rather than solving analytically for the unit cell, we could just integrate it directly,
Integrating the current box dimensions is what erate/rescale does in order to be compatible with all other styles, not just trate. If this is what you mean by integrating directly, could the improved sophistication of your coupling then similarly also work with any other style?
If so, maybe the erate/rescale style could just be moved out of fix deform/pressure and into fix deform and perform your improved corrections. Legacy behavior would be preserved in the regular erate style. Just a thought.
If it is more careful than erate/rescale, there might be an argument to just have it replace that style.
It's more careful in the specific case of xx/yy/zz = TRATE + xy/xz/yz = ERATE, but I don't think it works for other cases where the equation for the box evolution would be different. There are probably some cases we could handle, but it's not a general solution.
Integrating the current box dimensions is what erate/rescale does in order to be compatible with all other styles, not just trate. If this is what you mean by integrating directly, could the improved sophistication of your coupling then similarly also work with any other style?
What I meant there is that we could have something like
fix 1 all deform const/gradu xy 0.1 yz 1.0 xx 0.5 yy -0.5
which would set the style of all 6 components to CONST_GRADU and integrate all of them according to $\dot{\mathbf{e}} = \mathbf{e}\cdot\nabla\mathbf{u}$ (handling coupled directions in a reversible way). Part of the idea here is that by using the const/gradu option, the user is giving us explicit control over the entire unit cell rather than just some components of it, and we have a guarantee that all dimensions are deforming as we expect for SLLOD (which is formulated around $\nabla \mathbf{u}$ being applied consistently to the entire unit cell). We could perhaps achieve the same thing with erate/rescale, but it would at least need to check that we have control over the necessary set of xy/xz/yz for the given flow as is currently done in this PR, and I'd have to think about how h_rate is set to make sure that's ok in the general case.
There are probably some cases we could handle, but it's not a general solution.
Could you recast the corrections as functions of the current box geometry so that the solution doesn't depend on the exact type of loading?
Could you recast the corrections as functions of the current box geometry so that the solution doesn't depend on the exact type of loading?
I'll take a look at the math, but my guess is that this would probably work if we're integrating the corrections directly, so long as we're careful about h_rate. It might get tricky to do reversibly since we'd want to be using some sort of Trotter splitting scheme, but we can hopefully take advantage of h not being updated yet to achieve that.