[gmx-users] Constant-Force Pulling of Ubiquitin

[Weiwei Tao]

Hi GMX Users,

I am using Gromacs (Version 4.5.5) to do constant-force pulling of
ubiquitin and it's a implicit model. My mdp file for pulling is shown as
following.

integrator          =  md
dt                  =  0.001    ; ps !
nsteps              =  500000 ; total 500 ps.

nstxout             =  100
nstvout             =  100
nstfout             =  100
nstlist             =  10
nstlog              =  100
nstcalcenergy =100

rlist               =  5
rvdw              =  5
rcoulomb = 5
coulombtype = cut-off
vdwtype = cut-off
table-extension =  5
bd_fric             =  0
ld_seed             =  -1
pbc                 =  no
ns_type             =  simple
constraints         = all-bonds
lincs_order         = 4
lincs_iter          = 1
lincs-warnangle     = 30

Tcoupl              = v-rescale
tau_t               = 1.0
tc-grps             =  Protein
ref_t               =  300

Pcoupl              =  no

gen_vel             =  yes
gen_temp            =  300
gen_seed            =  173529

comm_mode = Angular
nstcomm      =100


; IMPLICIT SOLVENT ALGORITHM
implicit_solvent         = gbsa

; GENERALIZED BORN ELECTROSTATICS
; Algorithm for calculating Born radii
gb_algorithm             = Still
; Frequency of calculating the Born radii inside rlist
nstgbradii               = 1
; Cutoff for Born radii calculation; the contribution from atoms
; between rlist and rgbradii is updated every nstlist steps
rgbradii                 = 5
; Dielectric coefficient of the implicit solvent
gb_epsilon_solvent       = 80
; Salt concentration in M for Generalized Born models
gb_saltconc              = 0
; Scaling factors used in the OBC GB model. Default values are OBC(II)
gb_obc_alpha             = 1
gb_obc_beta              = 0.8
gb_obc_gamma             = 4.85
gb_dielectric_offset     = 0.009
sa_algorithm             = Ace-approximation
; Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA
; The value -1 will set default value for Still/HCT/OBC GB-models.
sa_surface_tension       = 2.05016

; Pull code
pull            = constant_force
;Center of mass pulling using a linear potential and therefore a constant
force.
pull_geometry   = direction
pull_start      = yes       ; define initial COM distance > 0
pull_ngroups    = 1
pull_group1     = Chain_B
pull_group0     = Chain_A
pull_k1         = -500      ; kJ mol^-1 nm^-2
pull_vec1 = 0.0 0.0 1.0

However, after pulling simulation, it turns out the potential of this
system becomes lower rather than higher (from -10000 to -20000). It's very
wired since potential should become larger after pulling.
Here is the notification after g_energy:

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Potential                  -20734.5        230    1268.47   -1385.71
 (kJ/mol)

You may want to use the -driftcorr flag in order to correct for spurious
drift in the graphs. Note that this is not
a substitute for proper equilibration and sampling! You should select the
temperature in order to obtain fluctuation properties.

I wonder whether there is any problem with my mdp file.
Thank you so much!!

Best,
Vivian