[gmx-users] A question about simulating a box of ethane in GROMACS

[Phil (Yang) Song]

Hi, all

I am learning GROMACS 4.5.1 recently and encountered a problem that puzzled
me for weeks. I am wondering if someone can help me with this.

Here is the problem: I was trying to perform a MD simulation for 216 ethane
molecules in a box. I have first generated a box of ethane with random
position and orientation. Then, I wanted to run a minimization for the box
of ethane to get rid of close contact of the randomly ethane molecules.
However, I have found that some of the ethane molecules laid on top of each
other after minimization and infinity forces were resulted.

Since my mdp file was successfully used to generate other type of molecules
such as acetylene and ethylene and they should be correct. Also, I have
tried to use conjugate gradient instead of steepest decent for the
minimization, but this effort came out to be in vain.

Later, I thought the topology file should be the source of error. I have
carefully checked the topology file for couple of times but did not find any
obvious error. I have also tried to change the charges on each atom in the
ethane to 0 and then run the minimization. Again, the ethane was attracted
into each other. Hence, I think the vdw could be the problem. However, the
parameters of the vdw interaction are all from the opls-aa files that come
with gromacs and it should be fine.

Eventually, I am puzzled...

I am iterating the content of the content of mdp input file and itp file in
the email and also attaching the input and output of my result as a tarball.
Hope someone can help me to solve this problem.

Thank you in advance!

Best,
Phil Song

================================================================================
MDP file

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; preprocessing options

; title of the simulation
title        =

; include path
include      =

; cpp flag
define       =

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; energy minimization
integrator    = steep

; Start time and timestep in ps
tinit         = 0
dt            = 0.001
nsteps        = 500

; run continuation or reperforming
init_step     = 0

emtol         = 1000.0
emstep        = 0.01

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; periodic boundary conditions
pbc           = xyz

periodic_molecules = no

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; electrostatics and vdw options
; electrostatics method
coulombtype    = PME
rcoulomb       = 1.3

; fourier setup for PME
fourierspacing = 0.12
fourier_nx     = 0
fourier_ny     = 0
fourier_nz     = 0
pme_order      = 4
ewald_rtol     = 1e-5
optimize_fft   = yes

ewald_geometry = 3d

; vdw cutoff
rvdw-switch    = 0
rvdw           = 1.31

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; neighbor list options
; neighbor list frequency
nstlist       = 10

; neighbor search algorithm
ns_type       = grid

; neighbor list cut-off
rlist         = 1.30

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; thermostat and barostat
tcoupl            = v-rescale
tc_grps         = System
ref_t             = 50.5
tau_t             = 0.5

pcoupl            = berendsen
pcoupltype        = isotropic
nstpcouple        = 10
ref_p             = 1.01325
tau_p             = 1.0

; Using compressibility of ethanol at 273.1 (approximately)
; from J. Phys. Chem., 1963, 67 (10), pp 2160.2164

compressibility = 1.0e-5

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; velocty generation
gen_vel       = yes
gen_temp      = 50.5
; gen_seed      = 173529

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; output control options
; output frequency for coords, velocities and forces
nstxout       = 5000
nstvout       = 5000
nstfout       = 5000

; output frequency for energy
nstenergy     = 5000
; energygrps    =

; output frequency for log
nstlog        = 5000

; output options for trajectory
nstxtcout     = 5000
xtc_precision = 100000
; xtc_grps      =


================================================================================
ITP file

; Ethane:
; Assign of the atom index
;
;      H2   H6
;      |    |
; H3 - C1 - C5 - H7
;      |    |
;      H4   H8
;

[ moleculetype ]
;     name        nrexcl
    Ethane             3

[ atoms ]
;   nr      type   resnr   residu    atom    cgnr      charge          mass
     1  opls_135       1      ETH     CE1      1        -0.18
     2  opls_140       1      ETH     HE2      1         0.06
     3  opls_140       1      ETH     HE3      1         0.06
     4  opls_140       1      ETH     HE4      1         0.06
     5  opls_135       1      ETH     CE5      2        -0.18
     6  opls_140       1      ETH     HE6      2         0.06
     7  opls_140       1      ETH     HE7      2         0.06
     8  opls_140       1      ETH     HE8      2         0.06

[ bonds ]
;  ai  aj   funct                 c0            c1
    1   5       1
    2   1       1
    3   1       1
    4   1       1
    6   5       1
    7   5       1
    8   5       1

[ pairs ]
;  ai  aj   funct
    2   6
    2   7
    2   8
    3   6
    3   7
    3   8
    4   6
    4   7
    4   8

[ angles ]
;  ai    aj    ak  funct          c0            c1
;  H3-C-C
    2     1     5      1
    3     1     5      1
    4     1     5      1
;  H-C-H
    2     1     3      1
    2     1     4      1
    3     1     4      1
;  C-C-H3
    6     5     1      1
    7     5     1      1
    8     5     1      1
;  H-C-H
    6     5     7      1
    6     5     8      1
    7     5     8      1

[ dihedrals ]
;  ai    aj    ak    al  funct       c0        c1        c2        c3
c4
    2     1     5     6      3
    3     1     5     6      3
    4     1     5     6      3
    2     1     5     7      3
    3     1     5     7      3
    4     1     5     7      3
    2     1     5     8      3
    3     1     5     8      3
    4     1     5     8      3



-- 
Phil (Yang) Song
Room 509, 590 Comm. Ave.
Chem. Dept., Boston Univ.
Boston, MA, 02215, USA
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110906/82207392/attachment.html>