GROMACS version: 2021.3
GROMACS modification: No
Here post your question
I am simulating a multilayer graphene sheet. After a few ns of production run, the individual graphene layers move apart and resulting in loss of structural integrity of the multilayer/multisheet structure. I would like the multilayer sheet to be mobile to enable membrane permeation so position restraints can not be used to keep the structure intact.
I was wondering if distance restraints can be applied to keep the different layers within a certain distance (0.33-0.34 nm) so as to keep the structure intact? If so, if [constraints] type 1 or type 2 can be used as listed below?
[ constraints ]
1 1001 1 0.34
atoms 1 and 1001 are in adjacent sheets assuming each sheet has 1000 atoms, type 1, b0 = 0.34
If the above is not correct, can someone please suggest me a better approach.
A constraint is not the best idea, it will remove the degree of freedom completely and will keep the two atoms at a constant distance. That’s not very physical.
In Gromacs’ pull code, it’s possible to apply e.g. a semi-harmonic (flat-bottom) restraining potential to the centers of masses of two groups only in a chosen direction. It could look like the following if you want to, say, semi-harmonically restraint the Z-distance of groups 1 and 2 (defined separately):
pull-coord1-type = flat-bottom ; semi-harmonic potential
pull-coord1-geometry = distance
pull-coord1-groups = 1 2 ; define separately
pull-coord1-init = 0.34 ; above this value the groups will feel the force
pull-coord1-dim = N N Y ; only in Z dimension
pull-coord1-k = 1000 ; force constant
You can then define a separate term for each pair of consecutive layers.
Now, I understand that you want to model permeation through a membrane, so perhaps you don’t want your molecule to be aligned along the Z-axis and still have lateral mobility of the sheets.
At this point, you have a choice whether to (a) reparametrize/reinforce the interactions between the sheets, (b) find a variant where you can define the reference vector as the normal to the graphene plane (doable perhaps with virtual atoms but not trivial), (c) restrain the number of contacts (coordination) between the layers using Plumed, (d) settle for restraining the layers without lateral mobility, so just using distance between COMs in all dimensions (they will still be able to rotate with respect to each other though), or (e) using milder pulling conditions if you’re pulling the graphene through the bilayer (consider if maybe the dissolution of the layers could be a physical phenomenon too).