I would like to try out the deform mdp command to calculate the viscosity of an isotropic liquid (just water as a model system for now) but I am at a loss when it comes to some details. When I start with a rectangular box, and put in my mdp file, e.g.:
deform 0 0 0 0 0 0.005
I understand this means that i exert shear along the c(y) component of the box (which, i belive, is the y-component of the box-vector that originally points in the z-direction). As discussed in a previous post by MichelePellegrino and mattiafelicepalermo I then get the viscosity by dividing the off diagonal stress tensor by the corresponding shear rate. three questions:
with the deform command above, which element of the stress tensor (PXX, PXY, …) and which box dimension (together with the deform rate (here 0.005) to get the shear rate) would i then use, to get the viscosity?
with deform, do i want to use a non-cubic initial box to improve convergence? if so, which box dimension would I elongate here?
In the documentation it is mentioned that deform can be combined with anisotropic or semi-isotropic pressure scaling. For an isotropic liquid such as bulk water this would make no sense, and I just use normal isotropic pressure scaling, correct?
It depends, it’s the pressure component corresponding to the deformation rate component. In your case PYZ or PZY.
It should not matter, since viscosity is the response to the deformation rate and not the deformation itself.
I think that’s mostly due to the system not being in thermodynamic equilibrium, and not due to the consititutive relation for the liquid. I prefer to run NVT when computing viscosity.
thank you for your swift reply! - one more clarification: the shear rate that i need to calculate the viscosity from PYZ would in this case be the deform parameter 0.005 nm ps-1 divided by a particular box dimension - is that boxXX? If not, which one would that be?
are you sure about NVT? … you can do that with cos-acceleration, but with deform your viscosity is directly calculated from the pressure (or more precisely one component of the stress tensor), and when you do NVT you will basically never end up exactly with your desired pressure (here 1 atm), and often enough you get negative pressures (and components) which would make the viscosity negative … i seem to be missing something here…
sorry, please ignore my last msg … I assume the result depends on the direction of the shear (can be positive or negative), …
I guess the documentation for this feature is somewhat limited - I wonder, the deform algorithm that is implemented in gromacs - are there any papers, or other documentation outside the gromacs world that discuss the details there?
Unfortunately… I don’t think so. Most papers refer to DOLLS or SLLOD algorithms, which are employed in cases where th triclinic box is not subject to a shear deformation (if I am not mistaken). In Gromacs, the box is actually deformed, so no need for SLLOD/DOLLS (but you can look them up, if you want). This is the merge request with the latest updates of the deform code, if you which to have a look: https://gitlab.com/gromacs/gromacs/-/merge_requests/3141
Michele’s right. It, unfortunately, took me a decade to realize that DOLLS and SLLOD are constructs to work around a non-deforming box and that that is why I had never understood them. I simply couldn’t imagine that you would want to simulate shear flow without deforming the unit cell. With actual box deformation no modifications of the equations of motion are needed, one only needs to make sure that one removes the flow profile from the kinetic energy calculation and thermostat velocity scaling.
Thanks a lot for your insight! I was wondering if you could elaborate a bit on how to remove the flow profile from the kinetic energy calculation and thermostat velocity scaling in practice. How can this be achieved through the settings in the “mdp” file?