[gmx-users] free energy binding problem

Justin Lemkul jalemkul at vt.edu
Thu Jun 4 13:41:22 CEST 2015



On 6/4/15 4:38 AM, Daniele Veclani wrote:
> Dear Useres
>
> I am studying an antibiotic in water (box with antibiotic and 1000
> molecules of water), I am currently doing a calculation of free
> energy binding. My problem is that, for some values of lambda, I have this
> error:
>
> Step 1447079, time 2894.16 (ps)  LINCS WARNING
> relative constraint deviation after LINCS:
> rms 437.814673, max 1855.867310 (between atoms 3031 and 3033)
> bonds that rotated more than 30 degrees:
>   atom 1 atom 2  angle  previous, current, constraint length
>     3001   3002   53.3    0.1000   0.1000      0.1000
>     3017   3018   40.1    0.2721   0.1090      0.1090
>     3023   3024   90.0    0.1100   5.3888      0.1100
>     3023   3025   90.0    0.1100   3.4338      0.1100
>     3026   3027   90.0    0.1000   0.1903      0.1000
>     3031   3032   90.0    0.1105   5.8869      0.1090
>     3031   3033   90.0    0.3034 204.2554      0.1100
>
> Back Off! I just backed up step1447079b.pdb to ./#step1447079b.pdb.1#
>
> Back Off! I just backed up step1447079c.pdb to ./#step1447079c.pdb.1#
> Wrote pdb files with previous and current coordinates
>
> WARNING: Listed nonbonded interaction between particles 3014 and 3024
> at distance 2.258 which is larger than the table limit 2.200 nm.
>
> This is likely either a 1,4 interaction, or a listed interaction inside
> a smaller molecule you are decoupling during a free energy calculation.
> Since interactions at distances beyond the table cannot be computed,
> they are skipped until they are inside the table limit again. You will
> only see this message once, even if it occurs for several interactions.
>
> IMPORTANT: This should not happen in a stable simulation, so there is
> probably something wrong with your system. Only change the table-extension
> distance in the mdp file if you are really sure that is the reason.
>
>
> My .MDP file is:
>
>   Run control
> integrator               = sd       ; Langevin dynamics
> tinit                    = 0
> dt                       = 0.002
> nsteps                   = 2500000   ; 1 ns
> nstcomm                  = 100
> ; Output control
> nstxout                  = 500
> nstvout                  = 500
> nstfout                  = 0
> nstlog                   = 500
> nstenergy                = 500
> nstxout-compressed       = 0
> ; Neighborsearching and short-range nonbonded interactions
> cutoff-scheme            = verlet
> nstlist                  = 20
> ns_type                  = grid
> pbc                      = xyz
> rlist                    = 1.2
> ; Electrostatics
> coulombtype              = PME
> rcoulomb                 = 1.2
> ; van der Waals
> vdwtype                  = cutoff
> vdw-modifier             = potential-switch
> rvdw-switch              = 1.0
> rvdw                     = 1.2
> ; Apply long range dispersion corrections for Energy and Pressure
> DispCorr                  = EnerPres
> ; Spacing for the PME/PPPM FFT grid
> fourierspacing           = 0.12
> ; EWALD/PME/PPPM parameters
> pme_order                = 6
> ewald_rtol               = 1e-06
> epsilon_surface          = 0
> ; Temperature coupling
> ; tcoupl is implicitly handled by the sd integrator
> tc_grps                  = system
> tau_t                    = 1.0
> ref_t                    = 300
> ; Pressure coupling is on for NPT
> Pcoupl                   = Parrinello-Rahman
> tau_p                    = 1.0
> compressibility          = 4.5e-05
> ref_p                    = 1.0
> ; Free energy control stuff
> free_energy              = yes
> init_lambda_state        = 6
> delta_lambda             = 0
> calc_lambda_neighbors    = 1        ; only immediate neighboring windows
> ; Vectors of lambda specified here
> ; Each combination is an index that is retrieved from init_lambda_state for
> each simulation
> ; init_lambda_state      = 6
> vdw_lambdas              = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
> 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
> coul_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> ; We are not transforming any bonded or restrained interactions
> bonded_lambdas           = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> restraint_lambdas        = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> ; Masses are not changing (particle identities are the same at lambda = 0
> and lambda = 1)
> mass_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> ; Not doing simulated temperting here
> temperature_lambdas      = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> ; Options for the decoupling
> sc-alpha                 = 0.5
> sc-coul                  = no       ; linear interpolation of Coulomb (none
> in this case)
> sc-power                 = 1.0
> sc-sigma                 = 0.3
> couple-moltype           = QUI  ; name of moleculetype to decouple
> couple-lambda0           = vdw      ; only van der Waals interactions
> couple-lambda1           = none     ; turn off everything, in this case
> only vdW
> couple-intramol          = no
> nstdhdl                  = 10
> ; Do not generate velocities
> gen_vel                  = no
> ; options for bonds
> constraints              = h-bonds  ; we only have C-H bonds here
> ; Type of constraint algorithm
> constraint-algorithm     = lincs
> ; Constrain the starting configuration
> ; since we are continuing from NPT
> continuation             = yes
> ; Highest order in the expansion of the constraint coupling matrix
> lincs-order              = 12
>
> Is it a symptom of poor stability of the system (not good energy
> minimization, or not good nvt npt)?
>

Does the system run normally (i.e. without using free energy options)?  Did the 
charge decoupling proceed correctly (assuming you're doing this in separate 
stages because here you're only doing vdW)?  Do other lambda states run, or do 
they all crash?

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================


More information about the gromacs.org_gmx-users mailing list