On 7/09/2011 12:40 AM, Phil (Yang) Song wrote:
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.
That's almost certainly inappropriate. Atomic clashes are inevitable,
and EM will not fix severe problems.
Get a single molecule in a small box, and use genconf to replicate it.
Then use EM and equilibrate thoroughly to remove the residual ordering.
Judicious use of NPT will fix your density to whatever you want.
Mark
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
--
gmx-users mailing list gmx-users@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-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
