[gmx-users] Need some help with membrane simulation

Ali Khan akk5r at virginia.edu
Thu Jun 5 23:17:54 CEST 2014 

Is there a difference between Potential-switch and Potential-switch-verlet?


On Thu, Jun 5, 2014 at 5:12 PM, Justin Lemkul <jalemkul at vt.edu> wrote:

>
>
> On 6/5/14, 4:57 PM, Ali Khan wrote:
>
>> Justin,
>>
>> 4.6.5 does not have Potential-switch-Verlet. When I added the other
>> changes, I still had the same error.
>>
>>
> Likely something more fundamental is wrong with the system.  What was the
> outcome of EM?  Exact numbers produced by mdrun, please.
>
>
>  Are these recommendations ones that I should generally use when
>> implementing charmm in gromacs?
>>
>
> Yes.  I'm still experimenting with all the new Gromacs features, but as
> best I can tell, in terms of reproducing a "normal" CHARMM setup:
>
> cutoff-scheme = Verlet
> vdwtype = cutoff
> vdw-modifier = Potential-switch
> coulombtype = PME
> rlist = 1.2
> rvdw = 1.2
> rvdw-switch = 1.0
> rcoulomb = 1.2
>
> -Justin
>
>
>
>>
>> On Thu, Jun 5, 2014 at 4:29 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>
>>
>>>
>>> On 6/5/14, 2:47 PM, Ali Khan wrote:
>>>
>>>  I am using the charmm36 force field.
>>>>
>>>> *My energy minimization .mdp file:*
>>>>
>>>> ; minim.mdp - used as input into grompp to generate em.tpr
>>>> ; Parameters describing what to do, when to stop and what to save
>>>> integrator = steep ; Algorithm (steep = steepest descent minimization)
>>>> emtol = 1000.0   ; Stop minimization when the maximum force < 1000.0
>>>> kJ/mol/nm
>>>> emstep          = 0.01          ; Energy step size
>>>> nsteps = 50000000   ; Maximum number of (minimization) steps to perform
>>>>
>>>> ; Parameters describing how to find the neighbors of each atom and how
>>>> to
>>>> calculate the interactions
>>>> nstlist = 1 ; Frequency to update the neighbor list and long range
>>>> forces
>>>> ns_type = grid ; Method to determine neighbor list (simple, grid)
>>>> rlist = 1.2 ; Cut-off for making neighbor list (short range forces)
>>>> coulombtype = PME ; Treatment of long range electrostatic interactions
>>>> rcoulomb = 1.2 ; Short-range electrostatic cut-off
>>>> rvdw = 1.2 ; Short-range Van der Waals cut-off
>>>> pbc = xyz ; Periodic Boundary Conditions (yes/no)
>>>> vdw-type        = cutoff
>>>> vdw-modifier = Potential-shift-Verlet
>>>> cutoff-scheme = Verlet
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> *Temperature equilibration .mdp file:*
>>>>
>>>> title = NVT equilibration
>>>> define = -DPOSRES ; position restrain the protein
>>>> ; Run parameters
>>>> integrator = md ; leap-frog integrator
>>>> nsteps = 500000 ; 2 * 500000 = 1000 ps
>>>> dt = 0.002 ; 2 fs
>>>> ; Output control
>>>> nstxout = 0
>>>> nstvout = 0
>>>> nstenergy = 500
>>>> nstlog = 0
>>>> ; Bond parameters
>>>> continuation = no ; first dynamics run
>>>> constraint_algorithm = lincs ; holonomic constraints
>>>> constraints = all-bonds ; all bonds (even heavy atom-H bonds)
>>>> constrained
>>>> lincs_iter = 1 ; accuracy of LINCS
>>>> lincs_order = 4 ; also related to accuracy
>>>> ; Neighborsearching
>>>> ns_type = grid ; search neighboring grid cells
>>>> nstlist = 5 ; 10 fs
>>>> rlist = 1.2 ; short-range neighborlist cutoff (in nm)
>>>> rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
>>>> rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
>>>> rcoulomb-switch          = 0
>>>> rvdw-switch              = 0.8
>>>>
>>>>
>>> Use rvdw-switch = 1.0.
>>>
>>>
>>>   ; long-range cut-off for switched potentials
>>>
>>>> rlistlong                = 1.4
>>>>
>>>>
>>> Probably never had the desired effect in terms of using CHARMM force
>>> fields and isn't necessary.  The Verlet scheme with rlist = 1.2 should
>>> provide adequate buffer space.
>>>
>>>
>>>   ; Electrostatics
>>>
>>>> coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
>>>> pme_order = 4 ; cubic interpolation
>>>> fourierspacing = 0.16 ; grid spacing for FFT
>>>> ; Relative dielectric constant for the medium and the reaction field
>>>> epsilon_r                = 1
>>>> epsilon_rf               = 1
>>>> ; Temperature coupling is on
>>>> tcoupl = V-rescale ; modified Berendsen thermostat
>>>> tc-grps = POPC Water_and_ions ; two coupling groups - more accurate
>>>> tau_t = 0.1 0.1 ; time constant, in ps
>>>> ref_t = 310 310 ; reference temperature, one for each group, in K
>>>> ; Pressure coupling is off
>>>> pcoupl = no ; no pressure coupling in NVT
>>>> ; Periodic boundary conditions
>>>> pbc = xyz ; 3-D PBC
>>>> ; Dispersion correction
>>>> DispCorr = No ; account for cut-off vdW scheme
>>>> ; Method for doing Van der Waals
>>>> vdw-type                 = cutoff
>>>> vdw-modifier = Potential-shift-Verlet
>>>>
>>>>
>>> Use Potential-switch-Verlet.  I don't remember if this was added to the
>>> code in time for 4.6.5, but it works well in any of the 5.0 releases.
>>>
>>> -Justin
>>>
>>>
>>>   cutoff-scheme = Verlet
>>>
>>>> ; Velocity generation
>>>> gen_vel = yes ; assign velocities from Maxwell distribution
>>>> gen_temp = 310 ; temperature for Maxwell distribution
>>>> gen_seed = -1 ; generate a random seed
>>>>
>>>>
>>>>
>>>> On Thu, Jun 5, 2014 at 2:36 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>>
>>>>
>>>>
>>>>> On 6/5/14, 2:11 PM, Ali Khan wrote:
>>>>>
>>>>>   Hi all,
>>>>>
>>>>>>
>>>>>> I don't know if anyone has experienced this problem before, but I
>>>>>> figure I
>>>>>> ask. I have several times run a simulation on a 135A x 135A POPC
>>>>>> membrane
>>>>>> acquired fro the charmm-gui website. I wanted now to simulate one that
>>>>>> was
>>>>>> slightly larger (150A x 150A). I made the larger membrane on
>>>>>> charmm-gui,
>>>>>> and I did everything the same way as I did before (even used the same
>>>>>> .mdp
>>>>>> files). After energy minimization, I get an error during my
>>>>>> temperature
>>>>>> equilibration as follows:
>>>>>>
>>>>>> *Back Off! I just backed up step0c_n1.pdb to ./#step0c_n1.pdb.4#*
>>>>>> *Wrote pdb files with previous and current coordinates*
>>>>>>
>>>>>> *WARNING: Listed nonbonded interaction between particles 31617 and
>>>>>> 31619*
>>>>>> *at distance 3f which is larger than the table limit 3f 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.*
>>>>>>
>>>>>>
>>>>>>
>>>>>> *-------------------------------------------------------*
>>>>>> *Program mdrun, VERSION 4.6.5*
>>>>>> *Source code file: /home/ali/Downloads/gromacs-4.6.5/src/mdlib/pme.c,
>>>>>> line:
>>>>>> 851*
>>>>>>
>>>>>> *Fatal error:*
>>>>>> *5 particles communicated to PME node 1 are more than 2/3 times the
>>>>>> cut-off
>>>>>> out of the domain decomposition cell of their charge group in
>>>>>> dimension
>>>>>> x.*
>>>>>>
>>>>>> *This usually means that your system is not well equilibrated.*
>>>>>> *For more information and tips for troubleshooting, please check the
>>>>>> GROMACS*
>>>>>> *website at http://www.gromacs.org/Documentation/Errors
>>>>>> <http://www.gromacs.org/Documentation/Errors>*
>>>>>> *-------------------------------------------------------*
>>>>>>
>>>>>>
>>>>>> I referred to the gromacs website, and it stated that this error means
>>>>>> that
>>>>>> my system is blowing up. I have tried several things to fix it
>>>>>> (including:
>>>>>> remaking the system on charmm-gui, running energy minimization with
>>>>>> small
>>>>>> energy steps [emstep = 0.01, 0.001, 0.0001], energy minimizing to Fmax
>>>>>> of
>>>>>> 100 kJ kJ/mol/nm instead of 1000, and I have done all these
>>>>>> combinations
>>>>>> in
>>>>>> both single and double precision) and none of these methods worked. I
>>>>>> also
>>>>>> should add that sometimes I am able to perform the temperature
>>>>>> equilibration, but then I would get the same error during the pressure
>>>>>> equilibration. I doubt it is charmm-gui, because they are usually very
>>>>>> reliable in generating membranes.
>>>>>>
>>>>>>
>>>>>>   On the contrary, we've had numerous reports of crashing membranes,
>>>>>> due
>>>>>>
>>>>> to
>>>>> some clashes in the starting structures.  Manual modifications to the
>>>>> coordinates of the offending atoms (fractions of an Angstrom, really)
>>>>> usually solves it.  The output of EM tells you where the maximum force
>>>>> is.
>>>>>    How big is this force?  What atom is it acting upon?  That's where
>>>>> you
>>>>> start your investigation.
>>>>>
>>>>>
>>>>>    Does anyone have any suggestions? I can post any of my input files
>>>>> if
>>>>>
>>>>>  requested.
>>>>>>
>>>>>>
>>>>>>   Always post an .mdp file when a run is crashing.  Membranes can be
>>>>>>
>>>>> particularly challenging to get right; small errors in settings can
>>>>> have
>>>>> catastrophic consequences.
>>>>>
>>>>> -Justin
>>>>>
>>>>> --
>>>>> ==================================================
>>>>>
>>>>> Justin A. Lemkul, Ph.D.
>>>>> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>>>>>
>>>>> Department of Pharmaceutical Sciences
>>>>> School of Pharmacy
>>>>> Health Sciences Facility II, Room 601
>>>>> University of Maryland, Baltimore
>>>>> 20 Penn St.
>>>>> Baltimore, MD 21201
>>>>>
>>>>> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
>>>>> http://mackerell.umaryland.edu/~jalemkul
>>>>>
>>>>> ==================================================
>>>>>
>>>>> --
>>>>> Gromacs Users mailing list
>>>>>
>>>>> * Please search the archive at http://www.gromacs.org/
>>>>> Support/Mailing_Lists/GMX-Users_List before posting!
>>>>>
>>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>>>
>>>>> * For (un)subscribe requests visit
>>>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>>>> send a mail to gmx-users-request at gromacs.org.
>>>>>
>>>>>
>>>>>  --
>>> ==================================================
>>>
>>> Justin A. Lemkul, Ph.D.
>>> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>>>
>>> Department of Pharmaceutical Sciences
>>> School of Pharmacy
>>> Health Sciences Facility II, Room 601
>>> University of Maryland, Baltimore
>>> 20 Penn St.
>>> Baltimore, MD 21201
>>>
>>> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
>>> http://mackerell.umaryland.edu/~jalemkul
>>>
>>> ==================================================
>>> --
>>> Gromacs Users mailing list
>>>
>>> * Please search the archive at http://www.gromacs.org/
>>> Support/Mailing_Lists/GMX-Users_List before posting!
>>>
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>
>>> * For (un)subscribe requests visit
>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>> send a mail to gmx-users-request at gromacs.org.
>>>
>>>
> --
> ==================================================
>
> Justin A. Lemkul, Ph.D.
> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>
> Department of Pharmaceutical Sciences
> School of Pharmacy
> Health Sciences Facility II, Room 601
> University of Maryland, Baltimore
> 20 Penn St.
> Baltimore, MD 21201
>
> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
> http://mackerell.umaryland.edu/~jalemkul
>
> ==================================================
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>

More information about the gromacs.org_gmx-users mailing list