Calculating viscosity


I need to compute the viscosity of TNB (tri naphtalene benzene) fluid at different temperatures ranging from 100 to 400 C. Initially, the periodic simulation box containing 50 molecules was subject to a 10 ns NPT run for equilibration at the desired temperature. Then the following commands were activated to impose the deformation along x direction:

fix 1 all deform 100000 x scale 2.0 remap v
fix 2 all nvt/sllod temp 473 473 100

The deformation is applied over a total run-time of 10 ns. In addition to Pxy (and Pxz), the x-component of atoms’ velocity is outputted for the purpose of estimating its average gradient along a perpendicular direction (either y or z).
The problem now is the huge fluctuations observed for the pressure component, which doesn’t seem to be approaching a constant value. Could you please let me know what is causing this problem?
Regarding the alternative approach, would coupling “fix nve” and "fix viscosity” work? Is the scalar output of “fix viscosity” divided by (the time interval)*(cross sectional area) yield the flux in the prependecular direction?
I’d truly appreciate your comments and inputs.
Thanks in advance.

Best regards,

The examples/VISCOSITY folder of the LAMMPS distribution contains example inputs for different ways to determine viscosity. Please see the README file and the referenced sections in the manual and corresponding text book literature for more details.


Dear Axel,

Thanks a lot for making reference to the example directory, found it very helpful.

Best regards,


Further to my previous question, I have come across some other issues that I would like to hear your comments on.
Using the command lines below, I’m trying to estimate the viscosity through the Muller-Plathe approach:

read_restart Restart.10000000

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
timestep 1.0

fix 1 all nve
fix 2 all setforce NULL 0.0 NULL
velocity all set NULL 0.0 NULL

fix 3 all viscosity 500000 x z 6
compute vel all property/atom vx

Averaging of vx within each slab is done in MATLAB. The system was previously equilibrated through an NPT run, and the attempt here is to monitor the transfer of the x-component of the momentum along the z direction. Therefore, the particles’ velocity and force are set to zero in the y direction. However this causes the temperature to drop extensively. This issue can be resolved by scaling the velocity to a much higher value than the target temperature, but I’m not sure if that is the right solution.
In addition, regardless of the chosen frequency for the perturbation, the velocity gradient is in the order of 10^-6 (1/fs), being much bigger than what is expected and yielding a small value for the shear viscoity. Could you please let me know what I’m doing wrong here?
Thank you in advance for the valuable inputs.

Best regards,