[gmx-users] Negative pulling forces, periodic molecules

Aykut Erbas aerbas at ph.tum.de
Fri May 28 15:57:43 CEST 2010


Hi,

Thanks for your answer,
my pull_pbcatoms look like
pull_pbcatom0 belongs to the surface
pull_pbcatom2 belongs to the peptide

I will be happy if you could be more clear about
In principle pull_pbcatom0 can be any atom, it should not change a 
predefined pulling rate

thanks
A.

Here is the log file

Input Parameters:
   integrator           = md
   nsteps               = 1000000
   init_step            = 0
   ns_type              = Grid
   nstlist              = 40
   ndelta               = 2
   nstcomm              = 1
   comm_mode            = Linear
   nstlog               = 5000
   nstxout              = 5000
   nstvout              = 5000
   nstfout              = 5000
   nstenergy            = 5000
   nstxtcout            = 2000
   init_t               = 0
   delta_t              = 0.002
   xtcprec              = 1000
   nkx                  = 70
   nky                  = 39
   nkz                  = 49
   pme_order            = 4
   ewald_rtol           = 1e-05
   ewald_geometry       = 0
   epsilon_surface      = 0
   optimize_fft         = FALSE
   ePBC                 = xyz
   bPeriodicMols        = TRUE
   bContinuation        = FALSE
   bShakeSOR            = FALSE
   etc                  = Berendsen
   epc                  = Berendsen
   epctype              = Semiisotropic
   tau_p                = 1
   ref_p (3x3):
      ref_p[    0]={ 1.00000e+00,  0.00000e+00,  0.00000e+00}
      ref_p[    1]={ 0.00000e+00,  1.00000e+00,  0.00000e+00}
      ref_p[    2]={ 0.00000e+00,  0.00000e+00,  1.00000e+00}
   compress (3x3):
      compress[    0]={ 2.50000e-08,  0.00000e+00,  0.00000e+00}
      compress[    1]={ 0.00000e+00,  2.50000e-08,  0.00000e+00}
      compress[    2]={ 0.00000e+00,  0.00000e+00,  4.50000e-05}
   refcoord_scaling     = No
   posres_com (3):
      posres_com[0]= 0.00000e+00
      posres_com[1]= 0.00000e+00
      posres_com[2]= 0.00000e+00
   posres_comB (3):
      posres_comB[0]= 0.00000e+00
      posres_comB[1]= 0.00000e+00
      posres_comB[2]= 0.00000e+00
   andersen_seed        = 815131
   rlist                = 0.8
   rtpi                 = 0.05
   coulombtype          = PME
   rcoulomb_switch      = 0
   rcoulomb             = 0.8
   vdwtype              = Cut-off
   rvdw_switch          = 0
   rvdw                 = 0.8
   epsilon_r            = 1
   epsilon_rf           = 1
   tabext               = 1
   implicit_solvent     = No
   gb_algorithm         = Still
   gb_epsilon_solvent   = 80
   nstgbradii           = 1
   rgbradii             = 2
   gb_saltconc          = 0
   gb_obc_alpha         = 1
   gb_obc_beta          = 0.8
   gb_obc_gamma         = 4.85
   sa_surface_tension   = 2.092
   DispCorr             = No
   free_energy          = no
   init_lambda          = 0
   sc_alpha             = 0
   sc_power             = 0
   sc_sigma             = 0.3
   delta_lambda         = 0
   nwall                = 0
   wall_type            = 9-3
   wall_atomtype[0]     = -1
   wall_atomtype[1]     = -1
   wall_density[0]      = 0
   wall_density[1]      = 0
   wall_ewald_zfac      = 3
   pull                 = umbrella
   pull_geometry        = position
   pull_dim (3):
      pull_dim[0]=1
      pull_dim[1]=1
      pull_dim[2]=0
   pull_r1              = 1
   pull_r0              = 1.5
   pull_constr_tol      = 1e-06
   pull_nstxout         = 100
   pull_nstfout         = 100
   pull_ngrp            = 1
   pull_group 0:
     atom (12558):
        atom[0,...,12557] = {0,...,12557}
     weight: not available
     pbcatom              = 6278
     vec (3):
        vec[0]= 0.00000e+00
        vec[1]= 0.00000e+00
        vec[2]= 0.00000e+00
     init (3):
        init[0]= 0.00000e+00
        init[1]= 0.00000e+00
        init[2]= 0.00000e+00
     rate                 = 0
     k                    = 0
     kB                   = 0
   pull_group 1:
     atom (8):
        atom[0,...,7] = {12665,...,12672}
     weight: not available
     pbcatom              = 12668
     vec (3):
        vec[0]= 1.00000e+00
        vec[1]= 0.00000e+00
        vec[2]= 0.00000e+00
     init (3):
        init[0]=-5.39178e-01
        init[1]= 1.09576e+00
        init[2]= 0.00000e+00
     rate                 = 0.01
     k                    = 100
     kB                   = 100
   disre                = No
   disre_weighting      = Conservative
   disre_mixed          = FALSE
   dr_fc                = 1000
   dr_tau               = 0
   nstdisreout          = 100
   orires_fc            = 0
   orires_tau           = 0
   nstorireout          = 100
   dihre-fc             = 1000
   em_stepsize          = 0.01
   em_tol               = 10
   niter                = 20
   fc_stepsize          = 0
   nstcgsteep           = 1000
   nbfgscorr            = 10
   ConstAlg             = Lincs
   shake_tol            = 1e-04
   lincs_order          = 4
   lincs_warnangle      = 30
   lincs_iter           = 1
   bd_fric              = 0
   ld_seed              = 1993
   cos_accel            = 0
   deform (3x3):
      deform[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      deform[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      deform[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
   userint1             = 0
   userint2             = 0
   userint3             = 0
   userint4             = 0
   userreal1            = 0
   userreal2            = 0
   userreal3            = 0
   userreal4            = 0
grpopts:
   nrdf:     36476.3     321.985     27568.7
   ref_t:         300         300         300
   tau_t:         0.1         0.1         0.1
anneal:          No          No          No
ann_npoints:           0           0           0
   acc:               0           0           0
   nfreeze:           N           N           N
   energygrp_flags[  0]: 0 0 0
   energygrp_flags[  1]: 0 0 0
   energygrp_flags[  2]: 0 0 0
   efield-x:
      n = 0
   efield-xt:
      n = 0
   efield-y:
      n = 0
   efield-yt:
      n = 0
   efield-z:
      n = 0
   efield-zt:
      n = 0
   bQMMM                = FALSE
   QMconstraints        = 0
   QMMMscheme           = 0
   scalefactor          = 1
qm_opts:
   ngQM                 = 0
Table routines are used for coulomb: TRUE
Table routines are used for vdw:     FALSE
Will do PME sum in reciprocal space.

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
U. Essman, L. Perela, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen
A smooth particle mesh Ewald method
J. Chem. Phys. 103 (1995) pp. 8577-8592
-------- -------- --- Thank You --- -------- --------

Using a Gaussian width (1/beta) of 0.25613 nm for Ewald
Cut-off's:   NS: 0.8   Coulomb: 0.8   LJ: 0.8
System total charge: -0.000
Generated table with 900 data points for Ewald.
Tabscale = 500 points/nm
Generated table with 900 data points for LJ6.
Tabscale = 500 points/nm
Generated table with 900 data points for LJ12.
Tabscale = 500 points/nm
Generated table with 900 data points for 1-4 COUL.
Tabscale = 500 points/nm
Generated table with 900 data points for 1-4 LJ6.
Tabscale = 500 points/nm
Generated table with 900 data points for 1-4 LJ12.
Tabscale = 500 points/nm

Enabling SPC water optimization for 4595 molecules.

Configuring nonbonded kernels...
Testing x86_64 SSE support... present.



Will apply umbrella COM pulling in geometry 'position'
between a reference group and 1 group
Pull group 0: 12558 atoms, mass 137674.546
Pull group 1:     8 atoms, mass   113.117

Initializing LINear Constraint Solver

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
B. Hess and H. Bekker and H. J. C. Berendsen and J. G. E. M. Fraaije
LINCS: A Linear Constraint Solver for molecular simulations
J. Comp. Chem. 18 (1997) pp. 1463-1472
-------- -------- --- Thank You --- -------- --------

The number of constraints is 1219

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
S. Miyamoto and P. A. Kollman
SETTLE: An Analytical Version of the SHAKE and RATTLE Algorithms for Rigid
Water Models
J. Comp. Chem. 13 (1992) pp. 952-962
-------- -------- --- Thank You --- -------- --------

Center of mass motion removal mode is Linear
We have the following groups for center of mass motion removal:
  0:  DIAM

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
H. J. C. Berendsen, J. P. M. Postma, A. DiNola and J. R. Haak
Molecular dynamics with coupling to an external bath
J. Chem. Phys. 81 (1984) pp. 3684-3690
-------- -------- --- Thank You --- -------- --------

There are: 26458 Atoms

Constraining the starting coordinates (step 0)

Constraining the coordinates at t0-dt (step 0)
RMS relative constraint deviation after constraining: 1.09e-06
Initial temperature: 300.134 K

Started mdrun on node 0 Thu May 27 12:03:34 2010

           Step           Time         Lambda
              0        0.00000        0.00000

Grid: 17 x 9 x 11 cells
   Energies (kJ/mol)
           Bond        G96Bond          Angle       G96Angle    Proper Dih.
    4.90240e+04    1.53403e+02    9.28348e+04    2.23220e+02    1.41487e+02
 Ryckaert-Bell.  Improper Dih.          LJ-14     Coulomb-14        LJ (SR)
    4.83454e+05    5.14344e+01   -1.62957e+01    1.71089e+03   -2.68920e+04
   Coulomb (SR)   Coul. recip.   COM Pull En.      Potential    Kinetic En.
   -1.88099e+05   -3.41786e+04    2.58240e-12    3.78407e+05    8.04666e+04
   Total Energy    Temperature Pressure (bar)  Cons. rmsd ()
    4.58874e+05    3.00708e+02    2.08080e+03    1.09884e-06




chris.neale at utoronto.ca wrote:
> Read about the .mdp options pull_pbcatom0 and pull_pbcatom1
>
> -- original message --
>
> Hi everyone,
>
> I have a pulling code running on single machine with Gromacs 4.0.7
> The systems composes of a surface, a protein  and water
> I pull the protein from the terminal group on the surface laterally.
>
> With Gromacs 3.3.3 (on single processors ) the setup works perfectly and
> generates correct results for pulling positions compared to g_traj...
> To be able use periodic molecules (for the surface), I switched to
> gromacs 4.0.7 and I use the "position" option to reproduce the previous
> v3 results.
> The problem  is that the pullf.xvg gives negative forces...
> Normally, the pulling force should reach a pseudo-constant level  but in
> my case it just keep going toward minus infinity (unless I cease the run)
> Also when I checked position of the pulled group I see that the tangent
> of the displacement-time curve (which is th pulling rate) is also
> increasing afte.
> This is bit irrational since the pulling rate is set constant at the
> first place...
> For instance, with a pulling rate of 10nm/ns in 1ns you should have
> displacement of roughly 10nm.
> Untill 0.4 ns everything looks ok but then (when the peptide moves to
> the next periodic box) velocity increases and I have a displacement of
> more 150nm within 1ns...
> Also, I must say that in the trajectory the increasing velocity of the
> pulled can be seen clearly with naked-eye.
>
> Here is my mdp file;
>
> integrator               = md
> nsteps                   = 1000000
> dt                       = 0.002
> nstlist                  = 40
> nstxout                  = 5000
> nstvout                  = 5000
> nstfout                  = 5000
> nstxtcout                = 2000
> nstlog                   = 5000
> nstenergy                = 5000
> constraints              = hbonds
> ns_type                  = grid
> coulombtype              = pme
> pme_order                = 4
> fourierspacing           = 0.12
> rlist                    = 0.8
> rvdw                     = 0.8
> rcoulomb                 = 0.8
> energygrps               = DIAM Protein SO
> tcoupl                   = berendsen
> tc_grps                  = DIAM Protein SOL
> tau_t                    = 0.1 0.1 0.1
> ref_t                    = 300 300 300
> gen_vel                  = no
> Pcoupl                   = berendsen
> Pcoupltype               = semiisotropic
> compressibility          = 2.5E-8 4.5E-5
> tau_p                    = 1.0 1.0
> ref_p                    = 1.0 1.0
> comm_mode                = angular
> comm_grps                = DIAM
> periodic_molecules       = yes
> pbc                      = xyz
> ;PULLING
> pull                     = umbrella
> pull_start               = yes
> pull_geometry            = position
> pull_nstxout             = 100
> pull_nstfout             = 100
> pull_ngroups             = 1
> pull_group0              = DIAM
> pull_group1              = AA1
> pull_vec1                = 1.0 0.0 0.0
> pull_init1               = 0.0 0.0 0.0
> pull_rate1               = 0.01
> pull_k1                  = 100
>
>
>




More information about the gromacs.org_gmx-users mailing list