[gmx-users] setting up a simulation of an ionic liquid

[Justin Lemkul]

On 5/13/13 11:50 AM, Laura Leay wrote:
> All,
>
> I've seen a few threads about simulations in ionic liquds but have not come across anything that tells me what settings I should use in my mdp file. The system is nitric acid which has fully dissociated into NO3- and HO3+. The simulation will run fine with just the ions at low density under NVT. However, when I solvate the box with SPC water (using the Amber force field) the simulation energy minimises with a really maximum force, of the order of 10^5. If I try to run an NPT simulation it immediately crashes and the md.log fle reveals that the electrostatic potential was unreasonable high, resulting in NaN.
>

You should investigate which atom bears the maximum force; that will point to 
the source of your problem.

> Below is the mdp file I've been using. I've been using PME electrostatics. If anyone can suggest some changes to make I would appreciate it.
>

The force field dictates most of the settings, most notably the cutoff.  The 
value of rvdw seems wrong, at least.

-Justin

> Laura
>
> PS, Justin, you were right earlier, there was a problem with my toplogy causing renaming of residues, I just hadn't spotted it. Thanks for your reply.
>
> -------------
>
>
> ; VARIOUS PREPROCESSING OPTIONS
> title                    = Yo
> cpp                      = /usr/bin/cpp
> include                  =
> ;define                   = -DFLEXIBLE
>
> ; RUN CONTROL PARAMETERS
> integrator               = steep ;md for simulation, steep for Emin
> ; Start time and timestep in ps
> tinit                    = 0
> dt                       = 0.0001
> nsteps                   = 50000000 ;remove 3
> ; For exact run continuation or redoing part of a run
> init_step                = 0
> ; mode for center of mass motion removal
> comm-mode                = Linear
> ; number of steps for center of mass motion removal
> nstcomm                  = 1
> ; group(s) for center of mass motion removal
> comm-grps                =
>
> ; LANGEVIN DYNAMICS OPTIONS
> ; Temperature, friction coefficient (amu/ps) and random seed
> ;bd-temp                  = 300
> bd-fric                  = 0
> ld-seed                  = 1993
>
> ; ENERGY MINIMIZATION OPTIONS
> ; Force tolerance and initial step-size
> emtol                    = 100
> emstep                   = 0.01
> ; Max number of iterations in relax_shells
> niter                    = 20
> ; Step size (1/ps^2) for minimization of flexible constraints
> fcstep                   = 0
> ; Frequency of steepest descents steps when doing CG
> nstcgsteep               = 1000
> nbfgscorr                = 10
>
> ; OUTPUT CONTROL OPTIONS
> ; Output frequency for coords (x), velocities (v) and forces (f)
> nstxout                  = 0
> nstvout                  = 0
> nstfout                  = 0
> ; Checkpointing helps you continue after crashes
> nstcheckpoint            = 1000
> ; Output frequency for energies to log file and energy file
> nstlog                   = 50
> nstenergy                = 50
> ; Output frequency and precision for xtc file
> nstxtcout                = 50
> xtc-precision            = 1000
> ; This selects the subset of atoms for the xtc file. You can
> ; select multiple groups. By default all atoms will be written.
> xtc-grps                 =
> ; Selection of energy groups
> energygrps               =
>
> ; NEIGHBORSEARCHING PARAMETERS
> ; nblist update frequency
> nstlist                  = 10
> ; ns algorithm (simple or grid)
> ns_type                  = grid
> ; Periodic boundary conditions: xyz (default), no (vacuum)
> ; or full (infinite systems only)
> pbc                      = xyz
> ; nblist cut-off
> rlist                    = 0.9
> domain-decomposition     = no
>
> ; OPTIONS FOR ELECTROSTATICS AND VDW
> ; Method for doing electrostatics
> coulombtype              = PME
> rcoulomb-switch          = 0
> rcoulomb                 = 0.9
> ; Dielectric constant (DC) for cut-off or DC of reaction field
> epsilon-r                = 1
> ; Method for doing Van der Waals
> vdw-type                 = Cut-off
> ; cut-off lengths
> rvdw-switch              = 0
> rvdw                     = 1.5
> ; Apply long range dispersion corrections for Energy and Pressure
> DispCorr                 = EnerPres
> ; Extension of the potential lookup tables beyond the cut-off
> table-extension          = 1
> ; Spacing for the PME/PPPM FFT grid
> fourierspacing           = 0.12
> ; FFT grid size, when a value is 0 fourierspacing will be used
> fourier_nx               = 0
> fourier_ny               = 0
> fourier_nz               = 0
> ; EWALD/PME/PPPM parameters
> pme_order                = 4
> ewald_rtol               = 1e-05
> ewald_geometry           = 3d
> epsilon_surface          = 0
> optimize_fft             = no
>
> ; GENERALIZED BORN ELECTROSTATICS
> ; Algorithm for calculating Born radii
> gb_algorithm             = Still
> ; Frequency of calculating the Born radii inside rlist
> nstgbradii               = 1
> ; Cutoff for Born radii calculation; the contribution from atoms
> ; between rlist and rgbradii is updated every nstlist steps
> rgbradii                 = 2
> ; Salt concentration in M for Generalized Born models
> gb_saltconc              = 0
>
> ; IMPLICIT SOLVENT (for use with Generalized Born electrostatics)
> implicit_solvent         = No
>
> ; OPTIONS FOR WEAK COUPLING ALGORITHMS
> ; Temperature coupling
> Tcoupl                   = berendsen
> ; Groups to couple separately
> tc-grps                  = System
> ; Time constant (ps) and reference temperature (K)
> tau_t                    = 0.1
> ref_t                    = 300
> ; Pressure coupling
> Pcoupl                   = berendsen
> Pcoupltype               = semiisotropic
> ; Time constant (ps), compressibility (1/bar) and reference P (bar)
> tau_p                    = 1.0
> compressibility          = 0 4.5e-5
> ref_p                    = 1.0 1.0
> ; Random seed for Andersen thermostat
> andersen_seed            = 815131
>
> ; SIMULATED ANNEALING
> ; Type of annealing for each temperature group (no/single/periodic)
> annealing                = no
> ; Number of time points to use for specifying annealing in each group
> annealing_npoints        =
> ; List of times at the annealing points for each group
> annealing_time           =
> ; Temp. at each annealing point, for each group.
> annealing_temp           =
>
> ; GENERATE VELOCITIES FOR STARTUP RUN
> gen_vel                  = yes
> gen_temp                 = 300
> gen_seed                 = 1993
>
> ; OPTIONS FOR BONDS
> ;constraints              = all-bonds
> ; Type of constraint algorithm
> constraint-algorithm     = Lincs
> ; Do not constrain the start configuration
> unconstrained-start      = no
> ; Use successive overrelaxation to reduce the number of shake iterations
> Shake-SOR                = no
> ; Relative tolerance of shake
> shake-tol                = 1e-04
> ; Highest order in the expansion of the constraint coupling matrix
> lincs-order              = 4
> ; Number of iterations in the final step of LINCS. 1 is fine for
> ; normal simulations, but use 2 to conserve energy in NVE runs.
> ; For energy minimization with constraints it should be 4 to 8.
> lincs-iter               = 1
> ; Lincs will write a warning to the stderr if in one step a bond
> ; rotates over more degrees than
> lincs-warnangle          = 30
> ; Convert harmonic bonds to morse potentials
> morse                    = no
>
> ; ENERGY GROUP EXCLUSIONS
> ; Pairs of energy groups for which all non-bonded interactions are excluded
> energygrp_excl           =
>
>

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

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

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