Dear all
I am investigating the conformation of a molecule in both explicit
benzene (BNZ) and acetronitrile (ACN) solvent boxes. Experimental work
has shown the change in conformation to be dependent on the polarity of
the solvent. Before doing md with my molecule in the solvent, I wanted
to check the solvent was accurately described by the force field (GAFF)
and that the charges I used are reliable.
I obtained the charges for each solvent from the RESP database (
http://q4md-forcefieldtools.org/REDDB/index.php ) project number W-46
and then used Antechamber to read in the charges and assign the GAFF
parameters.
After initial minimization of the box that leap generated, and NVT
temperature equilibration from 0-300K I ran a 2 ns NPT MD to equilibrate
the density (with SHAKE turned on, dt = 0.002). I obtained a reasonable
match to the experimental values: For ACN a density of 0.71 g/cm^3 at
300K (experimental value is 0.78 at 293K) while for BNZ I got 0.84,
close to the experimental value of 0.88.
Then I tried to check the boiling point. I restarted the 300K MD run in
NPT, but adding a restraint on the temperature so that it varied from
300K to 400k over 2ns. I kept SHAKE on, but reduced the time step to
0.001.
Plotting the density as a function of time (and hence temperature) I
expected to see a linear decrease up to the boiling point, then a
sharper decrease in the density as the liquid boiled. Ideally this would
be at 355K for ACN 353 for BNZ . But when plotting the density vs time,
the density decreased more or less linearly with increasing temperature
and did not plummet at any point between 300 and 400K. Hence either the
boiling point is above 400K, or something is going wrong in my
simulation. Please find my output file for the BNZ 300K-400K MD run
below. Please tell me if I am doing something wrong.
Is boiling point a good indicator of whether the solvent is accurately
described, or would anyone else recommend something else?
Does anyone have solvent parameters that they have checked give a good
match to experimental values of density, boiling point etc, that I could
use?
Many thanks
Mike
-------------------------------------------------------
Amber 9 SANDER 2006
-------------------------------------------------------
| Run on 02/12/2008 at 21:38:52
[-O]verwriting output
File Assignments:
| MDIN: boiling.in
| MDOUT: boiling.out
|INPCRD: equil.rst
| PARM: bnz.prmtop
|RESTRT: boiling.rst
| REFC: refc
| MDVEL: mdvel
| MDEN: mden
| MDCRD: boiling.mdcrd
|MDINFO: mdinfo
|INPDIP: inpdip
|RSTDIP: rstdip
Here is the input file:
equilibration 2ns boiling
&cntrl
imin = 0, irest = 1, ntx = 5,
ntb = 2, pres0 = 1.0, ntp = 1,
taup = 2.0,
cut = 10, ntr = 0,
ntc = 2, ntf = 2,
temp0 = 300.0,
ntt = 3, gamma_ln = 1000.0,
nstlim = 2000000, dt = 0.001,
ntpr = 1000, ntwx = 10000, ntwr = 1000
nmropt=1
/
&wt type='TEMP0', istep1=0,istep2=2000000,
value1=300.0, value2=400.0
/
&wt type='END'
/
------------------------------------------------------------------------
--------
1. RESOURCE USE:
------------------------------------------------------------------------
--------
| Flags: MPI
getting new box info from bottom of inpcrd
| INFO: Old style inpcrd file read
| peek_ewald_inpcrd: Box info found
|Largest sphere to fit in unit cell has radius = 12.109
| New format PARM file being parsed.
| Version = 1.000 Date = 02/11/08 Time = 19:50:53
NATOM = 1884 NTYPES = 2 NBONH = 942 MBONA = 942
NTHETH = 1884 MTHETA = 942 NPHIH = 3768 MPHIA = 942
NHPARM = 0 NPARM = 0 NNB = 8164 NRES = 157
NBONA = 942 NTHETA = 942 NPHIA = 942 NUMBND = 2
NUMANG = 2 NPTRA = 2 NATYP = 2 NPHB = 0
IFBOX = 1 NMXRS = 12 IFCAP = 0 NEXTRA = 0
NCOPY = 0
| Memory Use Allocated
| Real 179286
| Hollerith 11463
| Integer 156938
| Max Pairs 542592
| nblistReal 22608
| nblist Int 87783
| Total 4697 kbytes
| Duplicated 0 dihedrals
| Duplicated 0 dihedrals
BOX TYPE: RECTILINEAR
------------------------------------------------------------------------
--------
2. CONTROL DATA FOR THE RUN
------------------------------------------------------------------------
--------
BNZ
General flags:
imin = 0, nmropt = 1
Nature and format of input:
ntx = 5, irest = 1, ntrx = 1
Nature and format of output:
ntxo = 1, ntpr = 1000, ntrx = 1, ntwr =
1000
iwrap = 0, ntwx = 10000, ntwv = 0, ntwe =
0
ioutfm = 0, ntwprt = 0, idecomp = 0, rbornstat=
0
Potential function:
ntf = 2, ntb = 2, igb = 0, nsnb =
25
ipol = 0, gbsa = 0, iesp = 0
dielc = 1.00000, cut = 10.00000, intdiel = 1.00000
scnb = 2.00000, scee = 1.20000
Frozen or restrained atoms:
ibelly = 0, ntr = 0
Molecular dynamics:
nstlim = 2000000, nscm = 1000, nrespa = 1
t = 0.00000, dt = 0.00100, vlimit = 20.00000
Langevin dynamics temperature regulation:
ig = 71277
temp0 = 300.00000, tempi = 0.00000, gamma_ln=1000.00000
Pressure regulation:
ntp = 1
pres0 = 1.00000, comp = 44.60000, taup = 2.00000
SHAKE:
ntc = 2, jfastw = 0
tol = 0.00001
NMR refinement options:
iscale = 0, noeskp = 1, ipnlty = 1, mxsub =
1
scalm = 100.00000, pencut = 0.10000, tausw = 0.10000
Ewald parameters:
verbose = 0, ew_type = 0, nbflag = 1, use_pme =
1
vdwmeth = 1, eedmeth = 1, netfrc = 1
Box X = 32.151 Box Y = 24.217 Box Z = 30.380
Alpha = 90.000 Beta = 90.000 Gamma = 90.000
NFFT1 = 32 NFFT2 = 24 NFFT3 = 30
Cutoff= 10.000 Tol =0.100E-04
Ewald Coefficient = 0.27511
Interpolation order = 4
------------------------------------------------------------------------
--------
3. ATOMIC COORDINATES AND VELOCITIES
------------------------------------------------------------------------
--------
BNZ
begin time read from input coords = 2200.000 ps
Begin reading energy term weight changes/NMR restraints
WEIGHT CHANGES:
TEMP0 02000000 300.000000 400.000000 0 0
RESTRAINTS:
** No restraint defined **
Done reading weight changes/NMR restraints
Number of triangulated 3-point waters found: 0
| Atom division among processors:
| 0 948 1884
Sum of charges from parm topology file = 0.00000000
Forcing neutrality...
| Running AMBER/MPI version on 2 nodes
------------------------------------------------------------------------
--------
4. RESULTS
------------------------------------------------------------------------
--------
---------------------------------------------------
APPROXIMATING switch and d/dx switch using CUBIC SPLINE INTERPOLATION
using 5000.0 points per unit in tabled values
TESTING RELATIVE ERROR over r ranging from 0.0 to cutoff
| CHECK switch(x): max rel err = 0.2738E-14 at 2.422500
| CHECK d/dx switch(x): max rel err = 0.8314E-11 at 2.736960
---------------------------------------------------
| Local SIZE OF NONBOND LIST = 279574
| TOTAL SIZE OF NONBOND LIST = 536086
NSTEP = 1000 TIME(PS) = 2201.000 TEMP(K) = 299.56 PRESS =
97.0
Etot = 2362.1625 EKtot = 1401.8845 EPtot =
960.2780
BOND = 295.9920 ANGLE = 412.1808 DIHED =
385.8918
1-4 NB = 616.8442 1-4 EEL = -25.9557 VDWAALS =
-1006.8309
EELEC = 282.1558 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 141.6988 VIRIAL = 92.1274 VOLUME =
23672.3981
Density =
0.8602
Ewald error estimate: 0.1781E-02
------------------------------------------------------------------------
------
NMR restraints: Bond = 0.000 Angle = 0.000 Torsion =
0.000
========================================================================
=======
NSTEP = 2000 TIME(PS) = 2202.000 TEMP(K) = 309.54 PRESS =
-173.8
Etot = 2430.7749 EKtot = 1448.5996 EPtot =
982.1753
BOND = 298.9117 ANGLE = 387.4329 DIHED =
433.6618
1-4 NB = 620.2818 1-4 EEL = -26.4547 VDWAALS =
-1016.9076
EELEC = 285.2494 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 136.8269 VIRIAL = 225.6200 VOLUME =
23664.8139
Density =
0.8605
Ewald error estimate: 0.1839E-02
------------------------------------------------------------------------
------
Etc etc...
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Sun Feb 17 2008 - 06:07:06 PST
