Hi everyone,
I would like to perform a MD simulation in the (N, V, σa = 0, T) ensemble (based on the scheme 4 from this paper: https://onlinelibrary.wiley.com/doi/10.1002/adts.201800177).
Therefore, the cell shape should be able to change, while the volume of the system is kept constant.
In another article, they mentioned that the volume and deviatoric stress (σa) were controlled with the Martyna-Tobias-Tuckerman-Klein (MTTK) barostat.
I was looking in the LAMMPS documentation for possible keywords to run the simulation under these conditions, but I could not set up the input properly. Is this type of simulation even possible in LAMMPS? Any advice is really appreciated!
Best regards
Alexa
You would also need to apply a constraint to the deformation to maintain the volume. And while at it, one should also implement a constant area deformation constraint, which can be useful for membranes and bilayers. I am not 100%, but I believe NAMD has those. I recall colleagues talking about it (a long time ago, when we were still at Penn) when preparing simulations for transmembrane proteins.
You will have to confirm that it does what you think it’s supposed to do.
Note that MTTK barostatting is performed natively in LAMMPS (specified by mtk yes which is the default option – mtk no implements the older Hoover barostat).
Note also that the mtk terms are required for the desired behaviour. I was able to quickly jerry-rig a compute pressure implementation that returns deviatoric stress (making the returned pressure vector trace-free). When the deviatoric pressure compute was coupled to an unmodified implementation of FixNH I observed:
with mtk no (Hoover barostat), no box dimensions changed at all
with mtk yes (MTTK barostat) volume was not preserved
I’m happy to reconstruct the modified compute pressure if desired, but I’d rather break modify one fundamental part of LAMMPS at a time.
This is a valid point but it should also be considered that there are two computes in compute pressure which are scalar and vector, both of which can be used depending on the coupling (iso or not) by fixNH. I don’t want a pointless discussion about whose quick and dirty modified compute pressure is correct or not, but modifying both compute_scalar (which should return 0 anyway) and compute_vector to compute the deviatoric stress, I could actually have a stable volume preserving dynamics with cell dimensions variation using aniso 0. 0. 5. mtk no options. The mtk yes, on the other hand, gives bogus results.
I must say that I am a bit confused at the moment and do not have much time to dig into this, but there is something to be careful about in the implementation.
I didn’t change my hacked compute pressure’s scalar output so you are probably right (that hacking the scalar and vector output together with mtk no gives sensible volume-preserving dynamics).
To be explicit about this – at the moment this has just been a bit of a “nerd snipe” for me. I haven’t actually read the deviatoric barostat paper linked by @Alexa. For now, it is up to them to check what I wrote in the PR or propose corrections (I may have more time in a few weeks).