GROMACS version: 2018.8
GROMACS modification: No
I have a molecule, and I need to apply pulling forces with opposite directions to a head atom and a tail atom to stretch the molecule to an almost straight configuration, similar to that used in Figure 1 of the paper J. Am. Chem. Soc. 2019, 141, 29, 11603–11613.
I tried with these pull settings, but it does not work as I expected.
pull = yes
pull-ngroups = 2
pull-group1-name = Left
pull-group2-name = Right
pull-ncoords = 2
pull-coord1-type = constant-force
pull-coord1-geometry = direction-periodic
pull_coord1_vec = -1 0 0
pull-coord1-dim = Y Y Y
pull_coord1_k = 10000
pull-coord1-groups = 1 2
pull-coord2-type = constant-force
pull-coord2-geometry = direction-periodic
pull_coord2_vec = 1 0 0
pull-coord2-dim = Y Y Y
pull_coord2_k = 10000
pull-coord2-groups = 1 2
The Left group is the head atom, the Right group is the tail atom.
Could you please suggest to me the correct pull settings? So I can apply pulling forces with opposite directions at the two ends of a molecule.
If the goal is to apply constant force to pull the 2 groups away from each other along a single direction, you should only need 1 pull coordinate (coord1 refers to the coordinate along which the groups are pulled, not the groups themselves). Also, if you’re using direction-periodic, you should only include pull_coord1_vec (setting pull_coord1_vec to -1 0 0 means group 1 is pulled in the negative x-direction and group 2 is pulled in the positive x-direction), no need to include pull_coord1-dim (this is used for pull_coord1_geometry being set to distance, setting it to Y Y Y means it’s pulling in all 3 dimensions). For an example, see the mdp setup in this tutorial, and for a list of all possible mdp parameters, you can check the mdp options listed here.
Hello,
Thank you very much for the detailed reply! I tried with your suggestions and made a new .mdp file, and it works.
I have a molecule in a simulation box, and under NVT simulation, it will shrink or collapse. My goal is to pull the two ends of the molecule and to make the molecule straight again. In addition, I would like to know how much work is done by the pull force or the (potential) energy change of the molecule between the two states. However, I found it is hard to understand the output for me. Could you please help me further with this issue?
The COM Pull En. from gmx energy is about -335.808960 kJ/mol at the last timestep. And from pull_pullf.xvg and pull_pullx.xvg, I found the force is recorded as - 100 kJ/mol/nm and the moving distance of the pulling force is -3.35809 nm. I think the work done by the pulling force is calculated by F*r, and it seems it is reflected by COM Pull En.. However, I do not understand why COM Pull En. is a negative value. The potential energy of the initial state is about 1050 kJ/mol, and the final value after pulling is about 700 kJ/mol, so roughly the new potential energy equals initial potential energy plus COM Pull En., i.e., 700 = 1050 + (-335).
Could you please tell me how can I calculate the work done by the pulling force or the potential energy difference between the state before pulling and after pulling? I expected the potential energy to increase after pulling, but it decreases according to the output from gmx energy.