[gmx-users] stepsize too small ... but potential energy negative!

Luca Mollica luca.mollica at ibs.fr
Mon May 24 13:57:19 CEST 2010


Another way to try to understand what is going on wrong is to cut away 
the residue 1 (atom 1)
and see what happens: in this case you will understand if the "guilty" 
is the residue 1 (as LINCS seems, at least, to suggest) or anything else 
in the protein.

Cheers

Luca




> Dear Luca, dear all,
> thank you for your hints. I made some trials with my systems and these are
> my answers to your questions:
>
> - my system is a protein (with or w/o ligand) in solvent (water SPC).
> Following your suggestions, I tried to perform an EM on the protein w/o
> ligand after the editconf step (i.e. I created the topology with pdb2gmx,
> created a cubic box with editconf, then I used grompp+mdrun to perform EM).
> I received the same error: the minimization starts, performs 37 steps, then
> stops with the "stepsize too small...." communication. The energy in this
> case is -3.06e+4 (one order of magnitude higher than the solvated system)
> but still negative, and the maximum force is still "inf on atom 1". So
> minimizing in vacuo does not help to solve the problem. This atom 1 is the
> Nter of the protein, it is belonging to a Ser and I did not charge it
> explicitly with pdb2gmx (i.e. I did not use the flag -ter) but it is bound
> to 3 H in the .gro file.
>
> - If I look at the protein with VMD or other visualizing tools, it seems to
> me that no major problems are present on the structure. In particular, atom
> 1 is not "broken" or something similar, I can't see no "aberrations". I
> don't know how to check if something goes wrong apart from visualization, do
> you know some tools that could automatically check the file? To the best of
> my knowledge, gmxcheck performs checks only on trajectories, or can I use it
> also on single structures?
>
> - I also tried to continue anyway after the minimization step with the PR MD
> in NVT, but it did not start because of a lot of LINCS error. It suggests me
> that some distortions are present in my structure, but if I cannot minimize
> it, how to relieve this problem?
>
> I copy&paste below the .log file for this new minimization in vacuo. The
> parameters I used are the same for the minimization in solvent.
> Any help will be appreciated.
> Thank you very much and best regards
>
> Input Parameters:
>     integrator           = steep
>     nsteps               = 50000
>     init_step            = 0
>     ns_type              = Grid
>     nstlist              = 5
>     ndelta               = 2
>     nstcomm              = 1
>     comm_mode            = Linear
>     nstlog               = 100
>     nstxout              = 100
>     nstvout              = 100
>     nstfout              = 0
>     nstenergy            = 100
>     nstxtcout            = 0
>     init_t               = 0
>     delta_t              = 0.001
>     xtcprec              = 1000
>     nkx                  = 70
>     nky                  = 70
>     nkz                  = 70
>     pme_order            = 4
>     ewald_rtol           = 1e-05
>     ewald_geometry       = 0
>     epsilon_surface      = 0
>     optimize_fft         = TRUE
>     ePBC                 = xyz
>     bPeriodicMols        = FALSE
>     bContinuation        = FALSE
>     bShakeSOR            = FALSE
>     etc                  = No
>     epc                  = No
>     epctype              = Isotropic
>     tau_p                = 1
>     ref_p (3x3):
>        ref_p[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>        ref_p[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>        ref_p[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>     compress (3x3):
>        compress[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>        compress[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>      compress[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>     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                = 1
>     rtpi                 = 0.05
>     coulombtype          = PME
>     rcoulomb_switch      = 0
>     rcoulomb             = 1
>     vdwtype              = Cut-off
>     rvdw_switch          = 0
>     rvdw                 = 1.2
>     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                 = no
>     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.1
>     em_tol               = 1000
>     niter                = 20
>     fc_stepsize          = 0
>     nstcgsteep           = 1000
>     nbfgscorr            = 10
>     ConstAlg             = Lincs
> shake_tol            = 0.0001
>     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:        6039
>     ref_t:           0
>     tau_t:           0
> anneal:          No
> ann_npoints:           0
>     acc:            0           0           0
>     nfreeze:           N           N           N
>     energygrp_flags[  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.320163 nm for Ewald
> Cut-off's:   NS: 1   Coulomb: 1   LJ: 1.2
> System total charge: -3.000
> Generated table with 4400 data points for Ewald.
> Tabscale = 2000 points/nm
> Generated table with 4400 data points for LJ6.
> Tabscale = 2000 points/nm
> Generated table with 4400 data points for LJ12.
> Tabscale = 2000 points/nm
> Generated table with 4400 data points for 1-4 COUL.
> Tabscale = 2000 points/nm
> Generated table with 4400 data points for 1-4 LJ6.
> Tabscale = 2000 points/nm
> Generated table with 4400 data points for 1-4 LJ12.
> Tabscale = 2000 points/nm
> Configuring nonbonded kernels...
> Testing x86_64 SSE2 support... present.
>
>
> Removing pbc first time
> Initiating Steepest Descents
> Max number of connections per atom is 27
> Total number of connections is 27956
> Max number of graph edges per atom is 4
> Total number of graph edges is 4086
> Started Steepest Descents on node 0 Mon May 24 12:03:01 2010
> Steepest Descents:
>     Tolerance (Fmax)   =  1.00000e+03
>     Number of steps    =        50000
> Grid: 7 x 7 x 7 cells
>             Step           Time         Lambda
>                0        0.00000        0.00000
>
>     Energies (kJ/mol)
>          G96Bond       G96Angle    Proper Dih.  Improper Dih.          LJ-14
>      2.40716e+03    1.23569e+03    6.85831e+02    8.70819e+01   2.11992e-314
>       Coulomb-14        LJ (SR)        LJ (LR)   Coulomb (SR)   Coul. recip.
>     2.11992e-314   -2.19860e+03   -1.52663e+02   -6.70307e+03   -2.59368e+04
>        Potential Pressure (bar)
>     -3.05754e+04            nan
>             Step           Time         Lambda
>                1        1.00000        0.00000
>
>             Step           Time         Lambda
>                2        2.00000        0.00000
>
>             Step           Time         Lambda
>                3        3.00000        0.00000
>
>             Step           Time         Lambda
>                4        4.00000        0.00000
>
>             Step           Time         Lambda
>                5        5.00000        0.00000
>
>             Step           Time         Lambda
>                6        6.00000        0.00000
>
>             Step           Time         Lambda
>                7        7.00000        0.00000
>
>             Step           Time         Lambda
>                8        8.00000        0.00000
>
>             Step           Time         Lambda
>                9        9.00000        0.00000
>
>             Step           Time         Lambda
>               10       10.00000        0.00000
>
> ..........(same thing until)
>
>            Step           Time         Lambda
>               37       37.00000        0.00000
>
>
> Stepsize too small, or no change in energy.
> Converged to machine precision,
> but not to the requested precision Fmax<  1000
>
> Steepest Descents converged to machine precision in 38 steps,
> but did not reach the requested Fmax<  1000.
> Potential Energy  = -3.05753603909389e+04
> Maximum force     =                   inf on atom 1
> Norm of force     =                   inf
>
> Anna
>
>
> Message: 2
> Date: Fri, 21 May 2010 14:44:04 +0200
> From: Luca Mollica<luca.mollica at ibs.fr>
> Subject: Re: [gmx-users] stepsize too small ... but potential energy
> 	negative!
> To: Discussion list for GROMACS users<gmx-users at gromacs.org>
> Message-ID:<4BF68014.1060706 at ibs.fr>
> Content-Type: text/plain; charset="iso-8859-1"
>
> Dear Anna,
> are you talking about a system "in vacuo" or in solvent ?
> If you have placed the protein in water w/o minimizing it before
> solvation, probabily an "in vacuo" minimization could be useful for your
> system before moving into the solvated case. Moreover, are you sure that
> the protein does not have "broken" residue or something like that ?
> Sometimes, the completion of topology creation step goes fine but
> something wrong (on the gromacs side, apart from visualization) with
> sidechains or portion of the proteins that are generated by other
> softwares/servers.
> I do not think it's a ligand problem, BTW.
>
> Cheers
>
> Luca
>
> Message: 3
> Date: Fri, 21 May 2010 14:46:28 +0200
> From: Luca Mollica<luca.mollica at ibs.fr>
> Subject: Re: [gmx-users] stepsize too small ... but potential energy
> 	negative!
> To: gmx-users at gromacs.org
> Message-ID:<4BF680A4.3060702 at ibs.fr>
> Content-Type: text/plain; charset="iso-8859-1"
>
> Moreover,
> the problem is not the potential energy, but the force that must converge.
> The "inf" there about force and about atom 1 tells me that the computed
> force has problems, indeed.
> Which is atom 1 ? Is the protein Nter neutral or charged ? Or is atom 1
> an atom from the ligand ?
>
> Cheers
>
> L
>
>
>
>
>    




More information about the gromacs.org_gmx-users mailing list