AMBER Archive (2008)

Subject: AMBER: Is boiling point a good way to validate solvent paramters?

From: Mike Wykes (wykesm_at_averell.umh.ac.be)
Date: Wed Feb 13 2008 - 08:08:25 CST

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_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu