[gmx-users] Inconsistent results in different clusters and cores

tarak karmakar tarak20489 at gmail.com
Fri May 17 12:20:25 CEST 2013


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.

Regards,
Tarak


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
>>>
>>>
>>>
>>>
>>>
>>
>>
>
>



More information about the gromacs.org_gmx-users mailing list