[gmx-users] forcefield problems - advice needed
Lubos Vrbka
lists at vrbka.net
Sun Feb 4 15:54:30 CET 2007
hi guys,
for quite a long time i am trying to 'create' a forcefield for aliphatic
alcohols (right now, for ethanol) to use with a water potential we
employ in our simulations.
there is a problem - the water potential uses combination rule 2, so i
thought i would just try to use FFGMX and OPLS and change it a bit. it
didn't seem to work. i also tried the FFAMBER, that uses exactly this
combination rule. unfortunately, it fails as well.
i'm getting all sorts of errors - segfaults, lincs failures, ... in most
cases, the whole simulation 'freezes' (kinetic energy of order of 1e-5,
for example) and nothing happens. i am also quite frequently getting
close contacts between atoms (within one molecule, or between two
different molecules).
it seems i am missing something in the forcefield, but i am not able to
identify what is it... i would be very grateful if anyone of you could
help me and shed more light in this.
if anyone has working alcohol potential (united or all atom type) with
comb. rule 2, i would be more than happy, if you could share it with me...
also a 'side question' - when using FFAMBER, there are dihedrals defined
with all atom names specified (e.g., CT CT CT CT), and also dihedrals
with wildcards (X CT CT X). for a CT-CT-CT-CT dihedrals, both of these
would be applicable - are they both used in the simulation, or does the
latter (X CT CT X) override the former definition (CT CT CT CT)?
once again, i would be very glad if anyone of you could help me. i went
through the manual many times, but i am becoming more and more confused...
i attach the files for a simulation that is failing for me. it's 2+2
pentanol molecules in a prismatic simulation box (i have it setup like
this because we are normally doing slab simulations). all necessary
files are included.
thanks. with best regards,
lubos
-------------- next part --------------
;
; forcefield by Lubos Vrbka for (contaminated) ice simulations
; Lubos Vrbka, 2006
; lubos (@) vrbka (.) net
;
; water parameters - nada/van der eerden
; implementation by marcelo carignano (cari (@) purdue (.) edu)
; nacl parameters - j. chem. phys. 100, 3757
; implementation by marcelo carignano (cari (@) purdue (.) edu)
; used with spc/e - consistency?
; alcohol parameters - AMBER parm99 parameters
; implementation by lubos vrbka (lubos (@) vrbka (.) net)
;
; changed gen-pairs to yes (consistently with amber FF) - should not affect
; NE6 water potential
#define _FF_ICE
[ defaults ]
; nbfunc comb-rule gen-pairs FudgeLJ FudgeQQ
1 2 yes 0.5 0.8333
; *************************************************************************
[ atomtypes ]
;name type mass charge ptype c6 c12
; original AMBER atom types: CT, H1, HC, OH, HO
; CT sp3 carbon (charge corresponding to CH3!)
; H1 hydrogen bound to C with el. withdrawing substituent
; HC hydrogen bound to C
; HO hydroxyl hydrogen
; OH hydroxyl oxygen
; name btype mass charge ptype sigma epsilon
CT CT 12.010 0.1116 A 3.39967e-01 4.57730e-01
HC HC 1.0080 0.0372 A 2.64953e-01 6.56888e-02
H1 H1 1.0080 0.0372 A 2.47135e-01 6.56888e-02
;HO HO 1.0080 0.4215 A 2.47135e-01 6.56888e-02
HO HO 1.0080 0.4215 A 0.00000e+00 0.00000e+00
OH OH 16.0000 -0.6497 A 3.06647e-01 8.80314e-01
; NE6 water model
EP EP 0.00000 -0.866 D 0.00000E+00 0.00000E+00
OW OW 15.99940 0.000 A 3.11500E-01 0.714850050
HW HW 1.00800 0.477 A 6.73000E-02 0.115419008
LP LP 0.00000 -0.477 D 0.00000E+00 0.00000E+00
; sections [ ????types ] contain generic alcohol contaminant params
[ bondtypes ]
; i j func b0 kb
HO OH 1 0.09600 462750.4
CT OH 1 0.14100 267776.0
CT CT 1 0.15260 259408.0
CT H1 1 0.10900 284512.0
CT HC 1 0.10900 284512.0
[ angletypes ]
; i j k func th0 cth
CT OH HO 1 108.500 460.240
CT CT OH 1 109.500 418.400
CT CT CT 1 109.500 334.720
CT CT HC 1 109.500 418.400
CT CT H1 1 109.500 418.400
HC CT HC 1 109.500 292.880
H1 CT H1 1 109.500 292.880
H1 CT OH 1 109.500 418.400
[ dihedraltypes ]
; i j k l func coefficients
;CT CT OH HO 3 1.71544 0.96232 0.00000 -2.67776 0.00000 0.00000
;CT CT CT OH
;CT CT CT CT 3 3.68192 3.09616 -2.09200 -3.01248 0.00000 0.00000
;CT CT CT HC 3 0.66944 2.00832 0.00000 -2.67776 0.00000 0.00000
;CT CT CT H1 3 0.66944 2.00832 0.00000 -2.67776 0.00000 0.00000
;H1 CT OH HO
;HC CT CT H1
;HC CT CT HC 3 0.62760 1.88280 0.00000 -2.51040 0.00000 0.0000
;HC CT CT OH 3 1.04600 -1.04600 0.00000 0.00000 0.00000 0.00000
X CT CT X 3 0.65084 1.95253 0.00000 -2.60338 0.00000 0.00000
X CT OH X 3 0.69733 2.09200 0.00000 -2.78933 0.00000 0.00000
; *************************************************************************
[ moleculetype ]
; molname nrexcl
5OL 3
[ atoms ]
; id at type res nr resname at name cg nr charge
1 HO 1 5OL HO 1 0.420
2 OH 1 5OL OH 1 -0.648
3 CT 1 5OL C01 1 0.154
4 H1 1 5OL H11 1 0.037
5 H1 1 5OL H12 1 0.037
6 CT 1 5OL C02 1 -0.074
7 HC 1 5OL H21 1 0.037
8 HC 1 5OL H22 1 0.037
9 CT 1 5OL C03 1 -0.074
10 HC 1 5OL H31 1 0.037
11 HC 1 5OL H32 1 0.037
12 CT 1 5OL C04 1 -0.074
13 HC 1 5OL H41 1 0.037
14 HC 1 5OL H42 1 0.037
15 CT 1 5OL C05 1 -0.111
16 HC 1 5OL H51 1 0.037
17 HC 1 5OL H52 1 0.037
18 HC 1 5OL H53 1 0.037
[ bonds ]
; values defined after the defaults section
; i j funct b0 kb
1 2 1
2 3 1
3 4 1
3 5 1
3 6 1
6 7 1
6 8 1
6 9 1
9 10 1
9 11 1
9 12 1
12 13 1
12 14 1
12 15 1
15 16 1
15 17 1
15 18 1
[ angles ]
; values defined after the defaults section
; i j k funct th0 cth
1 2 3 1
2 3 4 1
2 3 5 1
2 3 6 1
3 6 7 1
3 6 8 1
3 6 9 1
6 9 10 1
6 9 11 1
6 9 12 1
9 12 13 1
9 12 14 1
9 12 15 1
12 15 16 1
12 15 17 1
12 15 18 1
4 3 5 1
4 3 6 1
7 6 8 1
7 6 9 1
10 9 11 1
10 9 12 1
13 12 14 1
13 12 15 1
16 15 17 1
16 15 18 1
[ dihedrals ]
; values defined after the defaults section
; i j k l func coefficients
1 2 3 4 3
1 2 3 5 3
1 2 3 6 3
2 3 6 7 3
2 3 6 8 3
2 3 6 9 3
3 6 9 10 3
3 6 9 11 3
3 6 9 12 3
6 9 12 13 3
6 9 12 14 3
6 9 12 15 3
9 12 15 16 3
9 12 15 17 3
9 12 15 18 3
4 3 6 7 3
4 3 6 8 3
4 3 6 9 3
5 3 6 7 3
5 3 6 8 3
5 3 6 9 3
7 6 9 10 3
7 6 9 11 3
7 6 9 12 3
8 6 9 10 3
8 6 9 11 3
8 6 9 12 3
10 9 12 13 3
10 9 12 14 3
10 9 12 15 3
11 9 12 13 3
11 9 12 14 3
11 9 12 15 3
13 12 15 16 3
13 12 15 17 3
13 12 15 18 3
14 12 15 16 3
14 12 15 17 3
14 12 15 18 3
;[ pairs ]
; i j funct c6 c12
;1 4 1
;1 5 1
;1 6 1
;1 7 1
;2 5 1
;2 6 1
;2 7 1
;3 6 1
;3 7 1
;4 7 1
; *************************************************************************
[ moleculetype ]
; molname nrexcl
NE6 1
[ atoms ]
; We use a strange order of atoms to make things go faster in GROMACS (?)
; id at type res nr resname at name cg nr charge
1 HW 1 NE6 HW1 1 0.477
2 HW 1 NE6 HW2 1 0.477
3 LP 1 NE6 LP1 1 -0.044
4 LP 1 NE6 LP2 1 -0.044
5 EP 1 NE6 EP 1 -0.866
6 OW 1 NE6 OW 1 0.000
[ constraints ]
; i j funct distance
1 6 1 0.0980 ; HW1 - OW
2 6 1 0.0980 ; HW1 - OW
1 2 1 0.15857 ; HW1 - HW2
[ dummies3 ]
; For the EP:
; Dummy from funct a b
5 6 1 2 1 0.199642537 0.199642537
; For the LPs:
; Dummy from funct a b c
3 6 1 2 4 -0.437172 -0.437172 8.022961206
4 6 1 2 4 -0.437172 -0.437172 -8.022961206
[ exclusions ]
1 2 3 4 5 6
2 1 3 4 5 6
3 1 2 4 5 6
4 1 2 3 5 6
5 1 2 3 4 6
6 1 2 3 4 5
-------------- next part --------------
;
; sample top file for NE6 water with contaminant simulation
;
; include parameters and topology
#include "_ffice.itp"
[ system ]
; Name
NE6 water with contaminant
[ molecules ]
; Compound #mols
;NE6 192
;NA 0
;CL 0
5OL 4
;2OL 10
-------------- next part --------------
NE6 water with contaminant
72
15OL HO 1 0.300 0.299 6.499
15OL OH 2 0.300 0.242 6.576
15OL C01 3 0.300 0.322 6.692
15OL H11 4 0.211 0.385 6.692
15OL H12 5 0.389 0.385 6.692
15OL C02 6 0.300 0.234 6.816
15OL H21 7 0.389 0.171 6.816
15OL H22 8 0.211 0.171 6.816
15OL C03 9 0.300 0.321 6.941
15OL H31 10 0.211 0.384 6.941
15OL H32 11 0.389 0.384 6.942
15OL C04 12 0.300 0.232 7.065
15OL H41 13 0.389 0.169 7.065
15OL H42 14 0.211 0.170 7.065
15OL C05 15 0.300 0.319 7.190
15OL H51 16 0.389 0.382 7.191
15OL H52 17 0.300 0.255 7.278
15OL H53 18 0.211 0.382 7.191
25OL HO 19 0.900 1.199 6.499
25OL OH 20 0.900 1.142 6.576
25OL C01 21 0.900 1.222 6.692
25OL H11 22 0.811 1.285 6.692
25OL H12 23 0.989 1.285 6.692
25OL C02 24 0.900 1.134 6.816
25OL H21 25 0.989 1.071 6.816
25OL H22 26 0.811 1.071 6.816
25OL C03 27 0.900 1.221 6.941
25OL H31 28 0.811 1.284 6.941
25OL H32 29 0.989 1.284 6.942
25OL C04 30 0.900 1.132 7.065
25OL H41 31 0.989 1.069 7.065
25OL H42 32 0.811 1.070 7.065
25OL C05 33 0.900 1.219 7.190
25OL H51 34 0.989 1.282 7.191
25OL H52 35 0.900 1.155 7.278
25OL H53 36 0.811 1.282 7.191
35OL HO 37 0.300 0.299 2.801
35OL OH 38 0.300 0.242 2.724
35OL C01 39 0.300 0.322 2.608
35OL H11 40 0.389 0.385 2.608
35OL H12 41 0.211 0.385 2.608
35OL C02 42 0.300 0.234 2.484
35OL H21 43 0.211 0.171 2.484
35OL H22 44 0.389 0.171 2.484
35OL C03 45 0.300 0.321 2.359
35OL H31 46 0.389 0.384 2.359
35OL H32 47 0.211 0.384 2.358
35OL C04 48 0.300 0.232 2.235
35OL H41 49 0.211 0.169 2.235
35OL H42 50 0.389 0.170 2.235
35OL C05 51 0.300 0.319 2.110
35OL H51 52 0.211 0.382 2.109
35OL H52 53 0.300 0.255 2.022
35OL H53 54 0.389 0.382 2.109
45OL HO 55 0.900 1.199 2.801
45OL OH 56 0.900 1.142 2.724
45OL C01 57 0.900 1.222 2.608
45OL H11 58 0.989 1.285 2.608
45OL H12 59 0.811 1.285 2.608
45OL C02 60 0.900 1.134 2.484
45OL H21 61 0.811 1.071 2.484
45OL H22 62 0.989 1.071 2.484
45OL C03 63 0.900 1.221 2.359
45OL H31 64 0.989 1.284 2.359
45OL H32 65 0.811 1.284 2.358
45OL C04 66 0.900 1.132 2.235
45OL H41 67 0.811 1.069 2.235
45OL H42 68 0.989 1.070 2.235
45OL C05 69 0.900 1.219 2.110
45OL H51 70 0.811 1.282 2.109
45OL H52 71 0.900 1.155 2.022
45OL H53 72 0.989 1.282 2.109
1.34724 1.55566 10.00000
-------------- next part --------------
;
; File 'mdout.mdp' was generated
; By user: spoel (291)
; On host: chagall
; At date: Mon Dec 15 13:13:06 2003
;
; VARIOUS PREPROCESSING OPTIONS
title = test for pentanol
cpp = /lib/cpp
include =
define =
; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.001
nsteps = 200000
; 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 = 50
; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout = 100
nstvout = 0
nstfout = 0
; Checkpointing helps you continue after crashes
nstcheckpoint = 0
; Output frequency for energies to log file and energy file
nstlog = 10
nstenergy = 0
; Output frequency and precision for xtc file
nstxtcout = 0
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 = 5
; 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.55
domain-decomposition = no
; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = pme
rcoulomb-switch = 0
rcoulomb = 0.55
; 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 = 0.55
; 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 = 3dc
epsilon_surface = 0
optimize_fft = no
; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
Tcoupl = nose-hoover
; Groups to couple separately
tc-grps = System
; Time constant (ps) and reference temperature (K)
tau_t = 0.1
ref_t = 10
; Pressure coupling
Pcoupl = no
; SIMULATED ANNEALING
; Type of annealing for each temperature group (no/single/periodic)
annealing = single
; Number of time points to use for specifying annealing in each group
annealing_npoints = 3
; List of times at the annealing points for each group
annealing_time = 0 10 20
; Temp. at each annealing point, for each group.
annealing_temp = 10 100 250
; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel = yes
gen_temp = 10
gen_seed = 1993
; OPTIONS FOR BONDS
constraints = none
; 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
More information about the gmx-users
mailing list