[gmx-users] Tabulated electrostatics not giving expected interactions

Kimberly Bowal klb83 at cam.ac.uk
Thu Oct 5 13:53:04 CEST 2017


Hello,


I have set up a simple simulation with 1000 water molecules and one K+
ion.  The ion is fixed in place at the centre of the box and I expect to
see clustering of the water molecules around the ion over time.  Instead I
only see a few weak interactions between one molecule and the ion and am
unsure what is wrong.

I am using a user-defined potential for both VdW and Coulomb interactions
(LJ and point charge) and have included the first lines below (the table
continues until r = 4).


I see in previous posts (
https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2015-October/101544.html)
that Gromacs will automatically neutralise a non-neutral system.  I
understand why this is done for Ewald methods but does this happen even
with a tabulated Coulomb potential?


Any suggestions as to what am I missing or doing incorrectly here would be
greatly appreciated!

I have attached my key input files as well.



Thank you,

Kimberly



# Parameters for Gromacs table

# Table : A/r^12 - C/r^6 + Dqq/r

# Units in nm kJ/mol

# LJ form

#   A   1.144406E-06   C   7.925970E-05

#

# r   f( r) = 1/r   -f'( r) = 1/r^2   g( r) = -1/r^6    g'( r) = -6/r^7
h( r) = 1/r^12   - h'( r) = 12/r^13

0.0000000000E+00   0.0000000000E+00   0.0000000000E+00   0.0000000000E+00
0.0000000000E+00   0.0000000000E+00   0.0000000000E+00

5.0000000000E-04   2.0000000000E+03   4.0000000000E+06
-6.4000000000E+19   -7.6800000000E+23   4.0960000000E+39   9.8304000000E+43

1.0000000000E-03   1.0000000000E+03   1.0000000000E+06
-1.0000000000E+18   -6.0000000000E+21   1.0000000000E+36   1.2000000000E+40

1.5000000000E-03   6.6666666667E+02   4.4444444444E+05
-8.7791495199E+16   -3.5116598080E+20   7.7073466293E+33   6.1658773034E+37

2.0000000000E-03   5.0000000000E+02   2.5000000000E+05
-1.5625000000E+16   -4.6875000000E+19   2.4414062500E+32   1.4648437500E+36

2.5000000000E-03   4.0000000000E+02   1.6000000000E+05
-4.0960000000E+15   -9.8304000000E+18   1.6777216000E+31   8.0530636800E+34
-------------- next part --------------
; VARIOUS PREPROCESSING OPTIONS
include                  = 
define                   = 

; RUN CONTROL PARAMETERS
integrator               = md-vv 
; Start time and timestep in ps
tinit                    = 0
dt                       = 0.001
nsteps                   = 1000000
; For exact run continuation or redoing part of a run
init_step                = 0
; Part index is updated automatically on checkpointing (keeps files separate) ##added by klb83 (in TT's)
simulation_part          = 1
; 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                =  ; default is whole syst


; LANGEVIN DYNAMICS OPTIONS
; Temperature, friction coefficient (amu/ps) and random seed
bd-fric                  = 0
ld-seed                  = 1093

; ENERGY MINIMIZATION and SHELL MD OPTIONS
; Force tolerance and initial step-size
emtol                    = 1.0 
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

; TEST PARTICLE INSERTION OPTIONS
rtpi                     = 0.05

; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout                  = 1000 
nstvout                  = 1000 
nstfout                  = 1000
; Output frequency for energies to log file and energy file
nstlog                   = 100
nstcalcenergy			 = -1
nstenergy                = 100
; Output frequency and precision for xtc file
nstxout-compressed       = 1000

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
cutoff-scheme       = group
nstlist             = 1
ns_type                  =  grid
; Periodic boundary conditions:
pbc                      = xyz
periodic_molecules       = no
; nblist cut-off        
rlist                    = 3.0
; long-range cut-off for switched potentials
rlistlong                = -1

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype              = User ; cut-off ; PME ; Ewald
; coulomb-modifier 	 	 = Potential-shift-Verlet ; Potential-shift ; None
rcoulomb_switch          = 0.0
rcoulomb                 = 3.0
; Relative dielectric constant for the medium and the reaction field
epsilon_r                = 1
epsilon_rf               = 1
; Method for doing Van der Waals
vdw-type                 = User ; cut-off ; PME
; vdw-modifier 			 = Potential-shift-Verlet ; Potential-shift ; None
; cut-off lengths       
rvdw_switch              = 0.0
rvdw                     = 3.0
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                 = no ; EnerPres
; Extension of the potential lookup tables beyond the cut-off
table-extension          = 1

; EWALD/PME/PPPM parameters
; 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
pme_order                = 4 
ewald_rtol               = 1e-05
ewald_geometry           = 3d
epsilon_surface          = 0

; Tables
energygrps = O H K
; Separate tables between energy group pairs
energygrp_table = 

; 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                 = 3.0
; Dielectric coefficient of the implicit solvent
gb_epsilon_solvent       = 80
; Salt concentration in M for Generalized Born models
gb_saltconc              = 0
gb_obc_alpha             = 1
gb_obc_beta              = 0.8
gb_obc_gamma             = 4.85
gb_dielectric_offset     = 0.009
sa_algorithm             = Ace-approximation
; Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA
; The value -1 will set default value for Still/HCT/OBC GB-models.
sa_surface_tension       = -1

; IMPLICIT SOLVENT (for use with Generalized Born electrostatics)
implicit_solvent         = no

; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling  
Tcoupl                   = nose-hoover
nsttcouple               = 1
nh-chain-length          = 10
; Groups to couple separately
tc-grps                  = System
; Time constant (ps) and reference temperature (K)
tau_t                    = 0.5 
ref_t                    = 500
; Pressure coupling     
Pcoupl                   = no 
Pcoupltype               = isotropic
nstpcouple               = -1
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau_p                    = 0.1 
compressibility          = 4.5e-5
ref_p                    = 1.0
; Scaling of reference coordinates, No, All or COM
refcoord_scaling         = no

; SIMULATED ANNEALING  
; Type of annealing for each temperature group
annealing                = no ; single ; periodic
; 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                 = 500
gen_seed                 = 1086

; OPTIONS FOR BONDS    
constraints              = all-bonds
; Type of constraint algorithm
constraint-algorithm     = Lincs
; Do not constrain the start configuration
continuation      = 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           = 

; WALLS                
; Number of walls, type, atom types, densities and box-z scale factor for Ewald
nwall                    = 0
wall_type                = 9-3
wall_r_linpot            = -1
wall_atomtype            = 
wall_density             = 
wall_ewald_zfac          = 3

; NMR refinement stuff 
; Distance restraints type: No, Simple or Ensemble
disre                    = no
; Force weighting of pairs in one distance restraint: Conservative or Equal
disre-weighting          = Conservative
; Use sqrt of the time averaged times the instantaneous violation
disre-mixed              = no
disre-fc                 = 1000
disre-tau                = 0
; Output frequency for pair distances to energy file
nstdisreout              = 100
; Orientation restraints: No or Yes
orire                    = no
; Orientation restraints force constant and tau for time averaging
orire-fc                 = 0
orire-tau                = 0
orire-fitgrp             = 
; Output frequency for trace(SD) to energy file
nstorireout              = 100

; Free energy control stuff
free-energy              = no
init-lambda              = 0
delta-lambda             = 0 
sc-alpha                 = 0
sc-power                 = 0
sc-sigma                 = 0.3
nstdhdl                  = 10
separate-dhdl-file       = yes
dhdl-derivatives         = yes
dh_hist_size             = 0
dh_hist_spacing          = 0.1
couple-moltype           = 
couple-lambda0           = vdw-q
couple-lambda1           = vdw-q
couple-intramol          = no

; Non-equilibrium MD stuff
acc-grps                 = 
accelerate               = 
freezegrps               = 
freezedim                = 
cos-acceleration         = 0

; Electric fields      
; Format is number of terms (int) and for all terms an amplitude (real)
; and a phase angle (real)
E-x                      = 
E-xt                     = 
E-y                      = 
E-yt                     = 
E-z                      = 
E-zt                     = 

; User defined thingies
user1-grps               = 
user2-grps               = 
userint1                 = 0
userint2                 = 0
userint3                 = 0
userint4                 = 0
userreal1                = 0
userreal2                = 0
userreal3                = 0
userreal4                = 0


More information about the gromacs.org_gmx-users mailing list