Many thanks for the suggestions and references David. Upon receipt of
your email I did a quick calculation of the enthalpy of vaporization
using the average potential energy (EPTOT in the md output file) in the
equation Delta-H(vap) = -<U>/N + RT and got values completely different
to experimental values. I then read the reference you sent me ( J. Chem.
Phys.V 120, P 9665-9678) and realized that it was the average
intermolecular potential energy that I was supposed to be using, which
admittedly made a lot more sense. So I repeated the calculation using
the sum of EELEC and VDWAALS as the intermolecular potential energy, the
results are summarised below, all energies in Kcal/mol. The averages are
taken from the last 500ps of a 1ns MD and I checked that all the
quantities had equilibrated.
<EEL> <NB> <PEinter> N mol <PE>/N RT Hvap
Hvap(expt) <density> density (expt)
benzene 285.767 -1007.47 -2014.94 157 -12.834
0.5961 13.43 8.10229446 0.855647 0.87381
acetonitrile -1001.11 -1011.52 -2012.63 291
-6.91625 0.5961 7.5124 8.01 0.721866 0.777
It seems that my acetonitrile parameters are not too bad, but that there
is clearly a problem with the intermolecular forces in my benzene model.
Searching for organic solvent in the amber archive, I found a link to
the Amber Parameter Database (
http://pharmacy.man.ac.uk/amber/ ) with a
box of all-atom acetonitrile provided by C. Jaime, X. Grabuleda, & P. A.
Kollman and a reference to their paper (J. Comp. Chem. 2000,21,901) in
which they calculate its density, heat of vaporization, and isothermal
compressibility and obtain good agreement with available experimental
data. Unfortunately nobody has uploaded a validated benzene model, so I
guess I'm going to have to tweak mine until I get better results.
On that note, how do I go about tweaking the benzene parameters? If I
want to get an enthalpy of vaporization value close to the experimental
value, my average intermolecular potential energy per molecule is going
have to change from -12.8 Kcal/mol to around -8.7 Kcal/mol. I can either
increase the positive and hence repulsive electrostatic interactions by
scaling the charges or play around with the VDW constants to weaken the
negative attractive VDWs interactions or do a combination of both.
I guess I'd have to try scaling the charges, then check the agreement of
the liquid properties with experimental results, then change the VDWs
params and check again, and then try combinations of changes in charges
& VDWs and check the properties again. This seems like a huge job
though.
Does anyone, perhaps with experience of tweaking solvent parameters have
any better ideas? Any help would be much appreciated..
Many thanks
Mike
-----Original Message-----
From: owner-amber.scripps.edu [mailto:owner-amber.scripps.edu] On Behalf
Of David A. Case
Sent: Wednesday, February 13, 2008 9:34 PM
To: amber.scripps.edu
Subject: Re: AMBER: Is boiling point a good way to validate solvent
paramters?
On Wed, Feb 13, 2008, Mike Wykes wrote:
>
> 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.
I would actually be a little concerned that the densities are so low --
this might be a signal that the vdW parameters need to be tweaked.
What one usually checks for energies is the enthalpy of vaporization,
which can be computed easily from the average total potential energy per
molecule of a liquid calculation at ordinary temperatures. Getting this
right helps to ensure that the attractive forces holding the liquid
together are about right.
>
> 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.
This is probably not a good way to get the boiling point, since it would
take a very long time for the phase equilibration to take place. A
better way would be to compute the free energy of disappearing a
particle from the liquid. This, combined with the (easily-computed)
free energy of re-creating it in the gas phase, would give you a free
energy of vaporization. Then doing a temperature scan would be more
likely to work. But, as I indicated above, go for the enthalphy of
vaporization first, since it is almost as good, and is
*much* easier to compute.
See this paper for the full theory here:
%A H.W. Horn
%A W.C. Swope
%A J.W. Pitera
%A J.D. Madura
%A T.J. Dick
%A G.L. Hura
%A T. Head-Gordon
%T Development of an improved four-site water model for biomolecular
simulations: TIP4P-Ew
%J J. Chem. Phys.
%V 120
%P 9665-9678
%D 2004
but note that most of the complicated corrections are small for rigid
solvent molecules. So, you can start with just
Delta-H(vap) = -<U>/N + RT
Where <U> is the average potential energy, N is the number of solvent
molecules. (The "polarization" correction can be important for polar
molecules like water, and maybe acetonitrile, less so for benzene. The
equation above will get you close, but the full theory is needed if you
want a very careful result.)
If you want to see how to get the vapor pressure and boiling point,
check out the second paper in this series:
%A H.W. Horn
%A W.C. Swope
%A J.W. Pitera
%T Characterization of the TIP4P-Ew water model: Vapor pressure and
boiling point %J J. Chem. Phys.
%V 123
%P 194504
%D 2005
...good luck....dac
--
==================================================================
David A. Case | e-mail: case.scripps.edu
Dept. of Molecular Biology, TPC15 | fax: +1-858-784-8896
The Scripps Research Institute | phone: +1-858-784-9768
10550 N. Torrey Pines Rd. | skype: dacase
La Jolla CA 92037 USA | http://www.scripps.edu/case
==================================================================
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
-----------------------------------------------------------------------
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:36 PST
