I have a molecule with a long chain in a simulation box. First, the system is run under NVT condition and the molecule chain will shrink or collapse to a more compact configuration. Then 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 compact and the straight configuration. However, I found it hard to understand the output.
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 configurations before pulling and after pulling? I expected the potential energy to increase after pulling, but it decreases according to the output from gmx energy.
Any help would be greatly appreciated!
I realized there are some more things that could be clarified. I would recommend outputting pullf.xvg data reasonably often - or use the pull-fout-average mdp option to get the average pull force since the last output. The distance can easily be calculated based on the pull-coord1-rate and the number of picoseconds between each pullf.xvg output.
Hello,
Thank you very much for the detailed replies and for sharing the papers! Sorry for the late reply.
I used a constant force for the pulling and the distance moved for the pulling force I think can be obtained from the pullx.xvg file. In this case, calculating the work done by the force can be done directly from F*r.
My confusion is that I found the output of COM Pull En. is a negative value. The absolute value of COM Pull En. is almost the same as the calculation from F*r, but I do not know why COM Pull En. is negative. The problem is that I found after the pulling, the potential energy of the system is smaller since COM Pull En. is negative. But I am expecting the potential energy of the molecule to be larger since the work done by the pulling force should increase the potential energy of the molecule.
I don’t think you should compare the COM Pull En. output with the work of the pulling potential. If they are “almost the same” but one with negative value, I would think it might be more of a coincidence. Does What does "COM-pull-energy" mean in the pulling process help?
Otherwise, could you post your COM Pull En. plot (and the output from gmx energy) along with the pullx.xvg plot and the value of the constant force you are pulling with?
Hello,
Thank you for the reply and for sharing the previous related post!
In my case, I used a constant pulling force, so I think calculating the work done by the pulling force is relatively easy.
I have checked several simulations with an identical initial condition but with different k values for pulling. It seems the COM Pull En. from gmx energy is also calculated by integrating F*r, but the value is negative. Maybe it is not a coincidence from my point of view.
I have attached the plots of COM Pull En. and pullx.xvg, pullf.xvg. Since the pulling force is constant, the figure of pullf.xvg is a straight line. By the way, in pullf.xvg, the sign of the pulling force is negative which also confuses me.
Any further suggestions would be greatly appreciated!
Hello,
Sure, here are the pulling settings I used. I tried with a few different settings and got this one based on insights from a response to a previous post of mine.
See Molecular dynamics parameters (.mdp options) - GROMACS 2024.3 documentation :
“For constant force pulling this is the force constant of the linear potential, and thus the negative (!) of the constant force in kJ mol-1 nm-1 (or kJ mol-1 rad-1 for angles).”
I think you will get the expected behaviour if you change to:
pull_coord1_vec = 1 0 0
pull_coord1_k = -100
But it’s possible you should keep pull_coord1_vec as you’ve set it. I’m not quite sure right now.
It seems the COM Pull En. is always negative, no matter the settings of the vector or k of the pulling force in the .mdp file.
The key problem is that I want to find a value like potential energy or similar ones to quantify the energy change of the molecule before and after the pulling force is applied. It would be greatly appreciated if you could provide any further suggestions.
