[gmx-users] Inconsistent results in different clusters and cores
Justin Lemkul
jalemkul at vt.edu
Fri May 17 13:23:07 CEST 2013
On 5/17/13 6:20 AM, tarak karmakar wrote:
> I have read the 'implementation of charmm in gromacs' by bjelkmar,
> JCTC. There they have used following cut-offs
> coulombtype=PME
> rcoulomb=1.2
> vdwtype=switch
> rvdw=1.2
> rvdw-switch=1.0
>
> I am not sure about rlist.
>
It's really unnecessary to spam the list with a dozen requests that ask the same
question.
By definition, with PME, rlist = rcoulomb. For CHARMM, use rlistlong = 1.4; the
message that is printed below is an error that has been fixed in newer versions.
-Justin
> On Fri, May 17, 2013 at 1:46 PM, tarak karmakar <tarak20489 at gmail.com> wrote:
>> Now if I increase the rlist it is showing that
>>
>> With coulombtype = PME, rcoulomb must be equal to rlist
>> If you want optimal energy conservation or exact integration
>> use PME-Switch
>>
>> I don't know exactly what will be the best set up with charmm force field.
>>
>> Any suggestion please!!
>>
>> Thanks,
>>
>> Tarak
>>
>>
>>
>>
>> On Fri, May 17, 2013 at 1:34 PM, tarak karmakar <tarak20489 at gmail.com>
>> wrote:
>>>
>>> What about Dispersion Correction ?
>>> But if I use this set of informations
>>>
>>> ; 7.3.3 Run Control
>>> integrator = md ; md integrator
>>> tinit = 0 ; [ps] starting time for
>>> run
>>> dt = 0.001 ; [ps] time step for
>>> integration
>>> nsteps = 5000000 ; maximum number of steps
>>> to integrate, 0.001 * 50,00,000 =5 ns
>>> nstcomm = 1 ; [steps] frequency of
>>> mass motion removal
>>> comm_grps = system ; group(s) for center of
>>> mass motion removal
>>>
>>> comm_mode = linear
>>>
>>>
>>> ; 7.3.8 Output Control
>>> nstxout = 5000 ; [steps] freq to write coordinates
>>> to trajectory
>>> nstvout = 5000 ; [steps] freq to write velocities to
>>> trajectory
>>> nstfout = 5000 ; [steps] freq to write forces to
>>> trajectory
>>> nstlog = 1000 ; [steps] freq to write energies
>>> to log file
>>> nstenergy = 1000 ; [steps] freq to write energies
>>> to energy file
>>> nstxtcout = 1000 ; [steps] freq to write
>>> coordinates to xtc trajectory
>>> xtc_precision = 1000 ; [real] precision to write xtc
>>> trajectory
>>> xtc_grps = System ; group(s) to write to xtc
>>> trajectory
>>> energygrps = protein ligand
>>>
>>> ; 7.3.9 Neighbor Searching
>>> nstlist = 10 ; [steps] freq to update neighbor
>>> list
>>> ns_type = grid ; method of updating neighbor list
>>> pbc = xyz ; periodic boundary conditions in
>>> all directions
>>> rlist = 1.2 ; [nm] cut-off distance for the
>>> short-range neighbor list
>>> rlistlong = 1.4
>>>
>>> ; 7.3.10 Electrostatics
>>> coulombtype = PME ; Particle-Mesh Ewald
>>> electrostatics
>>> rcoulomb = 1.2 ; [nm] distance for Coulomb
>>> cut-off
>>> fourierspacing = 0.16 ; [nm] grid spacing for FFT grid
>>> when using PME
>>> pme_order = 4 ; interpolation order for PME, 4 =
>>> cubic
>>> ewald_rtol = 1e-5 ; relative strength of
>>> Ewald-shifted potential at rcoulomb
>>>
>>> ; 7.3.11 VdW
>>> vdwtype = switch ; twin-range cut-off with rlist
>>> where rvdw >= rlist
>>> rvdw = 1.2 ; [nm] distance for LJ cut-off
>>> rvdw-switch = 1.0
>>>
>>> DispCorr = Ener ; apply long range dispersion
>>> corrections for energy
>>>
>>>
>>> ; 7.3.14 Temperature Coupling
>>> tcoupl = nose-hoover ; temperature coupling
>>> tc_grps = system ; groups to couple seperately to
>>> temperature bath
>>> tau_t = 1.0 ; [ps] time constant for coupling
>>> ref_t = 300 ; [K] reference temperature for
>>> coupling
>>>
>>> ; 7.3.15 Pressure Coupling
>>> pcoupl = parrinello-rahman ; pressure coupling where
>>> box vectors are variable
>>> pcoupltype = isotropic ; pressure coupling in
>>> x-y-z directions
>>> tau_p = 1.0 ; [ps] time constant for
>>> coupling
>>> compressibility = 4.5e-5 ; [bar^-1] compressibility
>>> ref_p = 1.0 ; [bar] reference pressure
>>> for coupling
>>>
>>> gen_vel = yes ; velocity generation
>>>
>>> gen_temp = 300
>>> gen_seed = 8877691
>>>
>>> ; 7.3.18 Bonds
>>> constraints = h-bonds ; covalent h-bonds constraints
>>> constraint_algorithm = LINCS ; LINear Constraint Solver
>>> continuation = yes ; apply constraints to the start
>>> configuration
>>> lincs_order = 4 ; highest order in the expansion
>>> of the contraint coupling matrix
>>> lincs_iter = 1 ; number of iterations to correct
>>> for rotational lengthening
>>> lincs_warnangle = 30 ; [degrees] maximum angle that a
>>> bond can rotate before LINCS will complain
>>>
>>>
>>> It is showing the following warning
>>>
>>> "For energy conservation with switch/shift potentials, rlist should be 0.1
>>> to 0.3 nm larger than rvdw."
>>>
>>>
>>>
>>>
>>>
>>>
>>> On Sun, May 12, 2013 at 11:57 PM, tarak karmakar <tarak20489 at gmail.com>
>>> wrote:
>>>>
>>>> Oh !
>>>> Thanks a lot Justin. I'll rerun all my jobs with this corrected mdp.
>>>> Restrains things I didn't follow properly, anyway I'll read about this.
>>>>
>>>>
>>>> On Sun, May 12, 2013 at 11:27 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>>>
>>>>>
>>>>>
>>>>> On 5/12/13 1:53 PM, tarak karmakar wrote:
>>>>>>
>>>>>> Thanks,
>>>>>>
>>>>>> I have used CGENFF force field parameters for the ligand generated from
>>>>>> PARMCHEM with 0 penalties. For protein I have used CHARMM36 force
>>>>>> field.
>>>>>> my npt.mdp file is as follows,
>>>>>>
>>>>>> ; 7.3.3 Run Control
>>>>>> integrator = md
>>>>>
>>>>>
>>>>> Bug 1021 was only relevant with md-vv, so it is not your problem here.
>>>>>
>>>>>
>>>>>> tinit = 0
>>>>>> dt = 0.001
>>>>>> nsteps = 5000000
>>>>>> nstcomm = 1
>>>>>> comm_grps = system
>>>>>> comm_mode = linear
>>>>>>
>>>>>>
>>>>>> ; 7.3.8 Output Control
>>>>>> nstxout = 5000
>>>>>> nstvout = 5000
>>>>>> nstfout = 5000
>>>>>> nstlog = 1000
>>>>>> nstenergy = 1000
>>>>>> nstxtcout = 1000
>>>>>> xtc_precision = 1000
>>>>>> xtc_grps = System
>>>>>> energygrps = lIG Protein Water
>>>>>>
>>>>>> ; 7.3.9 Neighbor Searching
>>>>>> nstlist = 10
>>>>>> ns_type = grid
>>>>>> pbc = xyz
>>>>>> rlist = 1.2
>>>>>>
>>>>>> ; 7.3.10 Electrostatics
>>>>>> coulombtype = PME
>>>>>> rcoulomb = 1.2
>>>>>>
>>>>>> ; 7.3.11 VdW
>>>>>> vdwtype = cut-off
>>>>>> rvdw = 1.2
>>>>>> DispCorr = EnerPres
>>>>>>
>>>>>
>>>>> Your short-range settings are incorrect for strict use of CHARMM. You
>>>>> should set:
>>>>>
>>>>> vdwtype = switch
>>>>> rvdw-switch = 1.0
>>>>> rlistlong = 1.4
>>>>>
>>>>> Your other settings for rlist, rcoulomb, and rvdw are fine.
>>>>>
>>>>>
>>>>>> ; 7.3.13 Ewald
>>>>>> fourierspacing = 0.12
>>>>>> pme_order = 4
>>>>>> ewald_rtol = 1e-5
>>>>>>
>>>>>> ; 7.3.14 Temperature Coupling
>>>>>> tcoupl = nose-hoover
>>>>>> tc_grps = system
>>>>>> tau_t = 1.0
>>>>>> ref_t = 300
>>>>>>
>>>>>> ; 7.3.15 Pressure Coupling
>>>>>> pcoupl = parrinello-rahman
>>>>>> pcoupltype = isotropic
>>>>>> tau_p = 1.0
>>>>>> compressibility = 4.5e-5
>>>>>> ref_p = 1.0
>>>>>>
>>>>>> gen_vel = yes
>>>>>
>>>>>
>>>>> In the absence of any restraints, initial velocity generation can
>>>>> produce incorrect dynamics. This is why we use restraints. Thus far, your
>>>>> observations simply seem consistent with random effects of improper
>>>>> nonbonded parameters and/or equilibration.
>>>>>
>>>>> -Justin
>>>>>
>>>>>> gen_temp = 300
>>>>>> gen_seed = 8877691
>>>>>>
>>>>>> ; 7.3.18 Bonds
>>>>>> constraints = h-bonds
>>>>>> constraint_algorithm = LINCS
>>>>>> continuation = yes
>>>>>> lincs_order = 4
>>>>>> lincs_warnangle = 30
>>>>>>
>>>>>>
>>>>>> On Sun, May 12, 2013 at 11:11 PM, Justin Lemkul <jalemkul at vt.edu>
>>>>>> wrote:
>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> On 5/12/13 1:34 PM, tarak karmakar wrote:
>>>>>>>
>>>>>>>> Thanks Justin for the Quick and Helpful reply.
>>>>>>>>
>>>>>>>> Yes. If I am right, the chaotic behavior of the simulations is
>>>>>>>> inherent
>>>>>>>> and can be assessed statistically by generating several independent
>>>>>>>> trajectories and analyzing their similar outcomes. But with the same
>>>>>>>> '.mdp'
>>>>>>>> file I am getting TOO much different results, and that's where I
>>>>>>>> worry.
>>>>>>>> I'll surely try with the recent version of gromacs. But, for now, can
>>>>>>>> you
>>>>>>>> give me a little more informations about problems (bugs) with the
>>>>>>>> 4.5.5
>>>>>>>> version, related to my context?
>>>>>>>>
>>>>>>>>
>>>>>>> The proposed relationship to bug 1012 is unclear to me. The issue
>>>>>>> there
>>>>>>> was an incompatibility between an integrator and thermostat, with
>>>>>>> obvious
>>>>>>> differences in thermodynamic output. Assessing your system in this
>>>>>>> context
>>>>>>> is not helpful. If you want to assess whether different core counts
>>>>>>> or
>>>>>>> hardware produce problems, then you need a very simple test case (like
>>>>>>> a
>>>>>>> box of water) that shows significant differences in averaged
>>>>>>> observables.
>>>>>>> As I said before, maybe your ligand parameters are insufficiently
>>>>>>> accurate
>>>>>>> (how did you generate them?) or .mdp settings are incorrect. Without
>>>>>>> such
>>>>>>> information, there is little point in trying to debug anything.
>>>>>>>
>>>>>>>
>>>>>>> -Justin
>>>>>>>
>>>>>>> --
>>>>>>> ==============================**==========
>>>>>>>
>>>>>>>
>>>>>>> Justin A. Lemkul, Ph.D.
>>>>>>> Research Scientist
>>>>>>> Department of Biochemistry
>>>>>>> Virginia Tech
>>>>>>> Blacksburg, VA
>>>>>>> jalemkul[at]vt.edu | (540) 231-9080
>>>>>>>
>>>>>>> http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justin<http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin>
>>>>>>>
>>>>>>> ==============================**==========
>>>>>>>
>>>>>>> --
>>>>>>> gmx-users mailing list gmx-users at gromacs.org
>>>>>>>
>>>>>>> http://lists.gromacs.org/**mailman/listinfo/gmx-users<http://lists.gromacs.org/mailman/listinfo/gmx-users>
>>>>>>> * Please search the archive at http://www.gromacs.org/**
>>>>>>>
>>>>>>> Support/Mailing_Lists/Search<http://www.gromacs.org/Support/Mailing_Lists/Search>before
>>>>>>> posting!
>>>>>>>
>>>>>>> * Please don't post (un)subscribe requests to the list. Use the www
>>>>>>> interface or send it to gmx-users-request at gromacs.org.
>>>>>>> * Can't post? Read
>>>>>>> http://www.gromacs.org/**Support/Mailing_Lists<http://www.gromacs.org/Support/Mailing_Lists>
>>>>>>>
>>>>>
>>>>> --
>>>>> ========================================
>>>>>
>>>>> Justin A. Lemkul, Ph.D.
>>>>> Research Scientist
>>>>> Department of Biochemistry
>>>>> Virginia Tech
>>>>> Blacksburg, VA
>>>>> jalemkul[at]vt.edu | (540) 231-9080
>>>>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>>>>>
>>>>> ========================================
>>>>> --
>>>>> gmx-users mailing list gmx-users at gromacs.org
>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>>>> * Please search the archive at
>>>>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>>>> * Please don't post (un)subscribe requests to the list. Use the www
>>>>> interface or send it to gmx-users-request at gromacs.org.
>>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>>
>>>>
>>>>
>>>>
>>>>
>>>
>>>
>>
>>
--
========================================
Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
More information about the gromacs.org_gmx-users
mailing list