Reversible integration of SLLOD

[ssande7]

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/rescale ${xy} y trate ${yy} x trate ${xx} remap none
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/rescale $(-v_xy) y trate $(-v_yy) x trate $(-v_xx) remap v
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:

1. Account for the change in streaming velocity caused by updating particle positions. Without this, the thermal velocity is implicitly changed.
2. Add an option to update the unit cell with fix deform during post_integrate rather than end_of_step to 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).
3. Apply the position/velocity updates due to SLLOD in a reversible way (important for mixed flows).
4. 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.
5. Fix a bug where velocities were not re-mapped after a box flip or when unwrapping coordinates for calculations like angular momentum.
6. Set h_rate from the start for TRATE deform style so that the correct flow rate tensor is used by fix nvt/sllod in the first time step.

See the implementation details section below for additional notes on some of these changes (numbers match up).

Still to-do:

1. Updates to documentation - I've left this for now until the user interface is finalised, but see below for an example input script.
2. 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.
3. Fix merge conflicts with develop - I based this on the 29Aug2024_update4 release 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/rescale ${gamma_xy} &
  xz erate/rescale ${gamma_xz} &
  yz erate/rescale ${gamma_yz} &
  remap ${remap}

# 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

1. The cleanest way I found to handle velocity bias changing with position updates was to modify compute temp/deform to 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 of temp/deform by sllod became a hard requirement, so to preserve the ability to use other biases (and reduce duplicate code), temp/deform now uses an internal temperature compute (temp by default), which is only called at times when the bias due to fix deform has been removed. This also enables new possibilities like turning off the thermostat in some/all directions with compute temp/partial, and I think will work nicely with our plans for a molecular thermostat. The other benefit here is that we can add a kick yes/no flag 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).
2. For now, I've made the toggle for applying box deformation in post_integrate to be setting nevery = 0 in fix 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.
3. 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/no flag to fix nvt/sllod to 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 with peculiar yes that I first got reversibility working properly - lab frame was more sensitive).
4. 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/pressure style 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 point nearly_equal function 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 move nearly_equal from the UEF package (possibly with minor modifications for robustness) into math_extra or 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:
stable_29Aug2024_update4...ssande7:lammps:sllod-integration
Please let me know if you want me to try and rebase before going further.

[akohlmey]

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).

[ssande7]

@akohlmey no worries, I suspected that might be the case - I'll sort that out shortly.

[ssande7]

Should be fixed now.

[akohlmey]

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.

[ssande7]

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.

[athomps]

@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.

[ssande7]

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?

[athomps]

@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.

[ssande7]

@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.

[ssande7]

@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.

[Copilot]

Pull request overview

Copilot reviewed 48 out of 48 changed files in this pull request and generated 3 comments.

[Copilot]

The arguments to MPI_Allreduce are swapped. The source and destination buffers should be &tspin (send) and &tspin_all (receive), not the other way around. The current code tries to reduce from uninitialized tspin_all and place the result in tspin, which is incorrect.

Suggested change

	MPI_Allreduce(&tspin_all,&tspin,1,MPI_DOUBLE,MPI_SUM,world);
	MPI_Allreduce(&tspin,&tspin_all,1,MPI_DOUBLE,MPI_SUM,world);

[Copilot]

The arguments to MPI_Allreduce are swapped. The source and destination buffers should be tspin (send) and tspin_all (receive), not the other way around. The current code tries to reduce from uninitialized tspin_all and place the result in tspin, which is incorrect.

Suggested change

	MPI_Allreduce(tspin_all,tspin,3,MPI_DOUBLE,MPI_SUM,world);
	MPI_Allreduce(tspin,tspin_all,3,MPI_DOUBLE,MPI_SUM,world);

[Copilot]

Missing space between directive and version. The reStructuredText directive should be .. versionadded:: TBD (with double colon and space before TBD), not .. versionadded::TBD (no space).

Suggested change

	.. versionadded::TBD
	.. versionadded:: TBD

[ssande7]

@athomps all tests are passing now, including some extra ones I've added for the new integration method and peculiar frame option. I've also updated the docs to hopefully make everything clear, and updated the relevant nemd examples and HowTo guide. Note, using N = 0 as the flag for fix deform to update the box during post_integrate was initially just a placeholder, so I'm happy to change that to some other interface if it's not intuitive enough. I've referenced the arxiv paper in the fix nvt/sllod doc page for now, but I'm hoping it will be accepted by the journal soon so I can swap that to the published version.

For moving erate/rescale, the simplest solution I found was to just move cumulative_shift across to fix deform and leave the rest of set_extra alone, except for cumulative_remap which was no longer needed. I also noticed in the process though that fix deform/pressure wasn't calling force->kspace->setup() when the box was updated, which I think is a bug so I've added that call in. @jtclemm is that correct, or was it omitted for a reason? I've also made some small changes to fix deform/pressure to (I think) make sure it works correctly with post_integrate updates + run_style respa, but if there's anything I've missed that would fundamentally break in that scenario then we might need to disallow that case.

For now, when the tilt factors are mapped back to the closest unit cell, I've made it so that the correction to xz due to xy tilt only applies if we have erate/rescale for xz to keep status quo behavior in other cases, but if similar handling were extended to other cases then (combined with the updated flip code) it might allow the restriction on flips that would induce a continuous xz change to be dropped in cases where fix deform has control over the xz tilt. That wouldn't play nicely with styles like xz final/scale/delta though if users are expecting to end up with a specific tilt value, so I'm not sure what the preference is here.

I haven't had a chance to test the Kokkos version on GPU yet, but I will try to do that soon, and it's working fine with the OpenMP backend. It would be worth someone double checking what I've done though in case anything is suboptimal. On that note - is the sync at fix_nvt_sllod.cpp:247 needed? I included it copying from line 385 to be safe, but I think it could be skipped unless I'm missing something. It also might be worth having a fused kernel version of compute temp/deform/kk in the case that the internal temperature compute is compute temp/kk, since that will probably be the case in most simulations, but I'm not sure how much difference that would make or whether it would be worth the added complexity. Same for fix nvt/sllod/kk, and possibly also some places in the non-Kokkos versions - we know the streaming velocity from the position so we could duplicate and inline code instead of having multiple loops over nlocal, but I'm not sure at this stage whether it's worthwhile.

For the /eff variants, I think the integration scheme is correct/equivalent to the previous version, and I didn't notice any problems in the comparison tests I ran, but it might be worth one of the original authors double checking since I'm not too familiar with that package. I don't think fix nvt/sllod/eff gets much use though since it's apparently had a typo since 2021 (creating the temperature compute). To handle general internal temperature computes in compute temp/deform/eff, I've just checked for the /eff suffix and used that to decide whether the spin part needs to be calculated as a correction or whether it's already included. I think that's fine with the current version, but it would be broken by accelerator variants so maybe something else is preferred. I've also assumed that the check for mass being non-null in compute temp/deform/eff compute_scalar was needed for some reason so I kept something equivalent.

[athomps]

@ssande7 This is an impressive body of work. Thank you for doing it.

[ssande7]

Thanks @athomps , it's been good to see this all come together finally.

[jtclemm]

I also noticed in the process though that fix deform/pressure wasn't calling force->kspace->setup() when the box was updated, which I think is a bug so I've added that call in. @jtclemm is that correct, or was it omitted for a reason?

That was an oversight, so good catch!

I've also made some small changes to fix deform/pressure to (I think) make sure it works correctly with post_integrate updates + run_style respa, but if there's anything I've missed that would fundamentally break in that scenario then we might need to disallow that case.

I just finished reading through all the code and it looks good to me. I'm going to run a few tests early next week to check a few combinations of deform commands, but if I don't remember to reply that means all the tests passed.

Overall, really impressive to get all this in and working! Certainly not a trivial improvement.

[ssande7]

Thanks @jtclemm !
Happy to help out if anything ends up failing in those tests.