[gmx-users] work from pulling simulation

[chris.neale at utoronto.ca]

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.