How to extract the unconvected coordinates and the peculiar velocities in a shear flow system-LAMMPS

Currently, I experiment with a molecular water system. I perform nonequilibrium molecular dynamics runs in LAMMPS, using the SLLOD equations of motion in conjunction with the Lees-Edwards periodic boundary conditions.
I perform NEMD simulationS for several shear rate values under the NVT ensemble to reach a steady-state. After
the system reaches a nonequilibrium steady state, I want to compute the diffusion coefficient by two different
techniques: either the Green-Kubo relations (which relate integrals of the momentum autocorrelation functions over the nonequilibrium steady-state to the diffusion coefficients directly); or the Einstein relations (which relate the unconvected mean
square displacements of the molecules at time t to products of the diffusion coefficients and t); To do so, the unconvected coordinates and the SLLOD peculiar momenta must be used.

I have understood that using the command
dump traj_dump all custom 1 file.lammpstrj id mol element type xu yu zu x y z vx vy vz fx fy fz
the convected coordinates, and the actual momenta are printed out, i.e., the terms corresponding to the external force are included.

I wonder if there is any way of directly extracting the unconvected coordinates and the peculiar velocities of the atoms every N timesteps. Are there any specific commands that I can use?

Below is the input script that I use, setting the nonequilibrium run simulation.
I really appreciate any help you can provide.

units real
atom_style full
log
dimension 3
boundary p p p

read_data restart_file
pair_style lj/cut/coul/long 10 10
kspace_style ewald 1.0e-4
pair_modify tail yes
pair_coeff 1 1 0.155297 3.16562
pair_coeff 2 2 0.0 0.000
pair_modify mix arithmetic
bond_style harmonic
bond_coeff 1 412.00935 1.0

angle_style harmonic
angle_coeff 1 45.739 109.47

special_bonds lj 0.0 0.0 0.0 coul 0.0 0.0 0.0

reset_timestep 0

neighbor 3.0 bin
neigh_modify every 1 delay 0 check yes

variable shear_rate index ${shear_rate}

fix integrator all nvt/sllod temp 300.0 300.0 100.0 # check thermostat
fix 2 all deform 1 xy vel ${shear_rate} remap v units box

dump traj_dump all custom 100 traj_${simname}.lammpstrj id mol element type xu yu zu x y z vx vy vz fx fy fz
dump_modify traj_dump sort id

run_style verlet
timestep 1
run 100000

LAMMPS does not have Lees-Edward boundary conditions. What fix deform does, however, results in equivalent trajectories.

There is no consistent diffusion happening because of the deformation.

There is no external force with applying fix deform.

The data you are looking for should be accessible either via post processing or by defining a sufficiently set up atom style variable for the velocity in shearing direction. You know the shearing rate which corresponds to the velocity at the box boundary; that would be the velocity to be subtracted from atoms at the same box boundary to get the relative velocity within the shearing removed. for atoms on the opposite box boundary this added term would be 0 and for atoms in between it can be obtained from a linear interpolation.

Equivalent to the velocity, you can also compute the displacement due to shearing and subtract it from the unwrapped position data. with a correspondingly processed trajectory, you could then do the kind of calculation of self-diffusion that you want. I cannot comment on the scientific value of such a property, though.

Thank you. Your answer was very helpful. I have tried to compute the streaming velocities corresponding to atoms via linear interpolation. It seems to work. However, I need a clarification regarding the streaming velocity-shear rate formula. Flow is in the x-direction and is linear in y with streaming velocity vx=gamma*y, where gamma is the shear rate. In my system, I perform a non-equilibrium run for vx=0.002 A/fs. After shearing for 20ns the dimensions of the system listed in the restart file are:

1.90329e+01 5.10919e+01 xlo xhi
1.90329e+01 5.10919e+01 ylo yhi
1.90329e+01 5.10919e+01 zlo zhi
-9.632e+00 0.0e+00 0.0e+00 xy xz yz

My question is what is the value of the shear rate. Is it given by gamma=0.002/51.09 1/fs??

Thank you in advance.

As a continuation of your answer regarding the computation of the displacement due to shearing, I have the following question.
The displacement due to shearing is given by
Delta qx(t)=int_{0}^{t}vx(t’)dt’, right?
Thus, I will need to have at my disposal the configuration path between the two instant times, 0-t. But if I understand correctly this is impossible, since it is highly demanding. Is there any other way to to compute the displacement due to shearing? Maybe directly through the LAMMPS?

The box dimensions do not represent in anyway the shearing. Fix deform will deform the box until there is a 45 degree angle and then it will “snap back” to a -45 degree angle and keep deforming. The xy tilt factor is representing the state of this deformation. If you visualize your trajectory, you will see what I mean. For more details, please study the LAMMPS manual, particularly the sections about fix deform and triclinic boxes. The shear rate is an input parameter to fix deform (which can be given in different ways).

Please note that I don’t have time to fill in for your adviser. I am happy to help with technical issues, and point out locations in the documentation, but if you need explanations beyond that, you need to talk to somebody else that is invested in your research.

If you are uncertain about details of the documentation or the method, I highly recommend to apply “the scientific method”: form a hypothesis, set up a (simple, fast) test simulation that would confirm that hypothesis or contradict it, check the results and have your confirmation or start another cycle checking out additional details or possible explanations why you didn’t get the expected result.