[lammps-users] Calculations of PMF using 'fix smd'

Dear all,

I have some questions about the calculations of PMF using ‘fix smd’. I used the following command in the input file after equilibration to calculate PMF.

fix pmf na smd cvel 20.0 0.0001 tether 0 0 10 0

variable dt equal 0.5

variable step equal 5000000

variable tt equal dt*step

variable distance equal 10

variable dr equal f_pmf[6]

fix out all print 1000 “{tt} {dr}” file umbrella_${distance}.pmf screen no

thermo_style custom step f_pmf[1] f_pmf[2] f_pmf[3] f_pmf[4] f_pmf[5] f_pmf[6] f_pmf[7]

run 5000000

Based on the manual, 'This fix computes a vector list of 7 quantities, the x-, y-, and z-component of the pulling force, the total force in direction of the pull, the equilibrium distance of the spring, the distance between the two reference points, and finally the accumulated PMF.

The values for the total force in direction of the pull seem to be incorrect since it decreases during all steps and not in agreement with the predicted PMF plots. Do you have any ideas where I made a mistake?

5000000 0 0 0 0 62.260006 62.259956 0

5001000 0.96437519 0.26431287 -0.010918119 -1 62.310006 62.259956 0.025025

5002000 1.9287504 0.52862575 -0.021836238 -2 62.360006 62.259956 0.10005

5003000 2.8931256 0.79293862 -0.032754357 -3 62.410006 62.259956 0.225075

5004000 3.8575007 1.0572515 -0.043672476 -4 62.460006 62.259956 0.4001

5005000 4.8218759 1.3215644 -0.054590595 -5 62.510006 62.259956 0.625125

5006000 5.7862511 1.5858772 -0.065508714 -6 62.560006 62.259956 0.90015

5007000 6.7506263 1.8501901 -0.076426833 -7 62.610006 62.259956 1.225175

5008000 7.7150015 2.114503 -0.087344951 -8 62.660006 62.259956 1.6002

5009000 8.6793767 2.3788159 -0.09826307 -9 62.710006 62.259956 2.025225

5010000 9.6437519 2.6431287 -0.10918119 -10 62.760006 62.259956 2.50025

5011000 10.608127 2.9074416 -0.12009931 -11 62.810006 62.259956 3.025275

5012000 11.572502 3.1717545 -0.13101743 -12 62.860006 62.259956 3.6003

5013000 12.536877 3.4360674 -0.14193555 -13 62.910006 62.259956 4.225325

5014000 13.501253 3.7003802 -0.15285367 -14 62.960006 62.259956 4.90035

5015000 14.465628 3.9646931 -0.16377178 -15 63.010006 62.259956 5.625375

5016000 15.430003 4.229006 -0.1746899 -16 63.060006 62.259956 6.4004

5017000 16.394378 4.4933188 -0.18560802 -17

The spring force is negative because it wants to contract, i.e. pull the moving group against the direction of pulling.
That seems logical because the target distance is increasing, but the actual distance is not.
This is the expected behavior at the beginning of a simulation when the group being pulled is at rest.
However at some point, it should pick up speed and then the force should oscillate.

There is no information about the details of your setup, so it is difficult to make any recommendations,
but I suggest you read up on the publications describing steered MD and particularly look up
how to determine a suitable force constant for the spring which is a crucial adjustable parameter for accurate results.

Thank you for your prompt response.

Here is more information about my system:
My model system consists of an electrolyte solution (1 M) confined within two planar graphene surfaces parallel to each other located 30 ̊A apart along the z-axis.
The graphene sheets form the boundaries of the simulation box along the z-direction. Force constant and velocity are considered 20 Kcal/mol-Angstrom and 0.0001 Angstroms/femtosecond, respectively. The PMF output is after 5 ns equilibration.

Kind regards,
Nasim

Thanks, but this doesn’t really help much because you don’t explain what you are pulling where. Also, I can read the input and thus I know what force constant and pulling velocity you have chosen, however it would be more important to know why you chose those values.

My assessment of what is happening and my recommendation to read up on the method and how to choose effective settings is unchanged.