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?
Sorry, I donât agree that the SLLOD equations of motion are not needed here. Lees-Edwards and Lagrangian-rhomboid (deformed) boundary conditions are equivalent and SLLOD IS needed in both cases. See Figure 4 of https://doi.org/10.1080/08927020601026629
To avoid having to go back to the original DOLLS/SLLOD papers by Hoover, Evans, et al., there is much more discussion of this in Todd and Daivisâ book on NEMD: Nonequilibrium Molecular Dynamics
A proper implementation of homogenous shear-driven NEMD would probably be of benefit in GROMACS. Iâm happy to help contribute if this is of interest to the community.
Note that such an implementation has been available in LAMMPS (which also uses fix deform rather than LE BCs) for quite a few years, but it has only recently been corrected and adapted to molecular systems by Stephen Sanderson at UQ: https://doi.org/10.1063/5.0315430
Noting âu ¡ âu = 0 under shear flow reveals that the Sllod equations are equivalent in that case to superimposing the expected velocity profile, qi ¡ âu, and then evolving under Newtonâs laws while simultaneously evolving the periodic unit cell in a manner commensurate with the flow.
This seems to consistent with what I write above and not using SLLOD. Or am I missing something?
This depends how you deform the box and what type of flow you impose. In theory, youâre correct that SLLOD isnât needed for planar shear flow (since itâs just a change in reference frame plus the initial kick in that case), but you certainly need it for elongational flows. In practice, even with planar shear flow you will start running into numerical integration problems, because the streaming velocity of a given atom changes as it moves, and not accounting for that effectively results in the thermal velocity jumping to a new value between half-steps.
The current implementation (similar to the âoldâ version of LAMMPS) may be âgood enoughâ for some cases (i.e. atomic systems at low shear rates under planar Couette flow). However, to compute the viscosity for large molecular systems, high shear rates, or elongational flows I would suggest changes would be needed to this part of the code.
For planar sheer flows the only issue I can see with the current implementation in GROMACS is that the velocity in leap frog is half a time step off from the coordinates. Thus the computed streaming velocity has a location error of 0.5*dt *v. This error is extremely small for any shear rate that is somewhat close to realistic shear rates. We might want to correct for this.
I have not thought of elongational flows with the deform option. Here deform is mainly intended for slow deformation of e.g. membranes to observe deformation and rupture. In such cases one does not want to modify the equations of motion. I will add a note in the manual that the deform option is not intended for simulating elongational flows.
There is also the issue of atomic vs. molecular thermostatting during the NEMD simulations. There will be another paper soon from Stephen looking at this, but for now some more info can be found in the NEMD book linked above. This is mostly an issue for large molecules and high shear rates.
Some other points spotted by Stephen during a quick look at the code. The box update also appears to come after the second velocity half step, which, unless Iâm mistaken, will cause some difficulties with forces not quite being right (and could make it more difficult to get lab-frame velocities). Also, the elongational flow isnât compatible with a constant flow tensor since you have e.g. h_xx(t) = h_xx(0) + t*e_xx rather than h_xx(t) = h_xx(0) * exp(e_xx * t). This might not be a problem for now if you are going to recommend against using for elongational flows.
I will read the publications you mentioned and I am happy to discuss this further. Since I am also using LAMMPS at the moment, I could run a simple benchmark (e.g. a water box or even a LJ box) to see if there is any difference, please let me know what should I look at.
About the elongational flow, if I recall correctly, the choice of using an âengineeringâ strain rate instead of the actual physical rate is simply to avoid instabilities (i.e. the box deforming exponentially fast). LAMMPS uses the âengineeringâ definition too, if I am not mistaken. From the documentation of fix deform:
The final , delta , scale , vel , and erate styles all change the specified dimension of the box via âconstant displacementâ which is effectively a âconstant engineering strain rateâ. This means the box dimension changes linearly with time from its initial to final value.
Maybe this could be pointed out in GROMACS documentation, rather than saying âdo not use the deform code for elongational flowsâ. Also becasue there are several cases where one would want to apply elongational strains in âliquidsâ (visco-elasto-plastic liquids).
I recall attending a talk a few years ago about how to produce âstableâ and âphysicalâ (as opposed to âengineeringâ) strains in MD, I will look for the related publication.
I recall attending a talk a few years ago about how to produce âstableâ and âphysicalâ (as opposed to âengineeringâ) strains in MD, I will look for the related publication.
Turns out it was the KR boundary conditions mentioned in the review by Todd and Daivis.
Hi Michele - I can help you set up some suitable benchmarks to compare the current implementations in LAMMPS and GROMACS. Probably best to handle this via email and then report back once weâve had a look? My email is je274@bath.ac.uk Cheers, James