[gmx-users] work from pulling simulation

[chris.neale at utoronto.ca]

Quoting:
> So for the above sample data, at t = .002, I get F = Spring Force =  
> AfmK * (18.042686- 18.042795) == AfmK * -0.000109
> dZ = Change in Position of spring = (18.042795-18.042684) == 0.000111
> Work by spring on Pulled group = F * dZ = AfmK * -0.000109 * 0.000111 ???

*********

I get something entirely different:

F=AfmK*(z-z_restraint). Therefore, over the interval between 0 fs and 2 fs:
F_umbrella_0.000 = AfmK * (18.042684-18.042684) = 0.0
F_umbrella_0.002 = AfmK * (18.042795-18.042686) = 0.000111 * AfmK = 0.327
F_umbrella_0.000_to_0.002 = (F_umbrella_0.000 + F_umbrella_0.002)/2 = 0.1635
dz = 18.042795 - 18.042684 = 0.000111
F*dz = 0.1635 * 0.000111 = 0.0000181485

Then sum all of these segments.

But I have never done this so you should not just rely on what I say.

Also note that you need Jarzynski and many repeats to get free  
energies from your work:
exp(-BdG)=lim[N->infinity] <exp(-BW_N)>, where
W_N=sum[F_i(z)dz]

When you get it working, please report back to the list.

Chris.

--- original message ---

Hi,
Thanks for the quick response, and I apologize for being too vague.
Here's a second attempt.

In Gromacs 3.3.1, I am running a steered MD simulation with commands,
mdrun -s topol -o traj.trr -e ener -c conf -g log \
       -pi pull.ppa -pn pull.ndx \
       -pd pull.pdo -po pullout.ppa

Here is my pull.ppa input file:
;
; GENERAL
verbose                  = yes
; Skip steps               = 1
; Runtype: afm, constraint, umbrella
runtype                  = afm
; Number of pull groups
ngroups                  = 1
; Groups to be pulled
group_1                  = hx23loop_A
; The group for the reaction force.
; reference_group          = Protein
; Weights for all atoms in each group (default all 1)
weights_1                =
reference_weights        =
; Ref. type: com, com_t0, dynamic, dynamic_t0
reftype                  = com
; Use running average for reflag steps for com calculation
reflag                   = 1
; Select components for the pull vector. default: Y Y Y
pull_dim                 = N N Y
; DYNAMIC REFERENCE GROUP OPTIONS
; Cylinder radius for dynamic reaction force groups (nm)
r                        = 0
; Switch from r to rc in case of dynamic reaction force
rc                       = 0
; Update frequency for dynamic reference groups (steps)
update                   = 1

; AFM OPTIONS
; Pull rates in nm/ps  = 1nm/1ns , if rate = 0.0010
afm_rate1                = 0.001
; Force constants in kJ/(mol*nm^2)
afm_k1                   = 3000
; Directions
afm_dir1                 = 0.0 0.0 1.0
; Initial spring positions
afm_init1                = 3.128236  2.156135  18.042684

The first few lines of the PDO output file are:
# AFM   3.0
# Component selection: 0 0 1
# nSkip 1
# Ref. Group ''
# Nr. of pull groups 1
# Group 1 'hx23loop_A'  afmVec 0.000000 0.000000 1.000000  AfmRate  
0.001000  AfmK 3000.000000
#####
0.000000        0.000000        18.042684       18.042684
0.002000        0.000000        18.042795       18.042686
0.004000        0.000000        18.042913       18.042688

Where the columns indicate the time, reference group, COM of the  
pulled group, and position of the spring.

So, my original question given this type of output, is how can one  
calculate the work done by the spring on the pulled protein group?

So for the above sample data, at t = .002, I get
F = Spring Force = AfmK * (18.042686- 18.042795) == AfmK * -0.000109
dZ = Change in Position of spring = (18.042795-18.042684) == 0.000111
Work by spring on Pulled group = F * dZ = AfmK * -0.000109 * 0.000111 ???

I can do this at every time step, and then sum the incremental work  
values to get the total work done.

Originally, I thought it would be more accurate to get the work by
using an average force.
For example, at t=1,
W = (F(t=1)+F(t=0))/2 * dZ(t=1),
but the results obtained this way seem to be worse -- meaning they
were very different than what I expected in my control systems.

If someone could please comment on whether this seems to be the
correct/incorrect way to calculate things, I would appreciate it.

I hope this was more clear, but please ask for more  
clarification/details if I something if needed.
Thank you,
Danny.


-- original message --

1. What is Z_bead?
2. Could you give us a real .pdo file segment?
3. I can only assume AFM, or something similar, but need to see a pull
code input file. Can you provide one?
4. What version of gromacs?

This is all to say I don't have enough information. But from what I can see:

1. If you are using force averages, your equation for W seems totally
incorrect. You should be averaging forces in a given histogram bin of
z and then integrating over the general reaction coordinate of z. But
again, I have no idea what Z_bead means.

2. If you are doing average force by current position vs umbrella then
be sure that you take the negative (ensemble average of system force
equals negative ensemble average of umbrella force: see Neale et al.
Chem Phys Lett 460, 375-381 (2008). -- and lots of other papers I
imagine).

Chris.

-- original message --

Hello GMXers,

I am performing steered MD simulations using the pull code in
1-Dimension (z) and I'd like to calculate the work done on my system (in
the z-direction) by the pulling spring.

If my output PDO pulling file looks something like this:

Time    Z_pulled_group	Z_bead
0	0		0
1	.5		1
2	1		2

Then my forces and work should be this?

Time    dZ	Spring Force	Work_on_system
0	0	0  * k		0  * 0  * k  ???
1	.5	.5 * k		.5 * .5 * k  ???
2	1	1  * k 		1  * 1  * k  ???

where
dZ = (Z_bead-Z_pulled_group)
F  = k*dz
W  = F*dz

Originally, I thought it would be more accurate to get the work by
using an average force.
For example, at t=1,
W = (F(t=1)+F(t=0))/2 * dZ(t=1),
but the results obtained this way seem to be worse -- meaning they
were very different than what I expected in my control systems.

If someone could please comment on whether this seems to be the
correct/incorrect way to calculate things, I would appreciate it.
Thanks,
Danny.