# Calculating heat of vaporization of solvent with MD simulation.

I was asked to calculate the solvation free energy of some solutes in Dimethylacetamide (DMA) with MD simulation. I disliked to develop own force field for DMAC. It’s mainly because I didn’t think my method to develop the force field was accurate. Since I used AMBER all the time, I tried to use GAFF for DMA at the beginning to testify whether it’s good for DMAC. I planed to use density and heat of vaporization to testify it.

Dr. Justin A. Lemkul at University of Maryland helped me a lot for this (http://mackerell.umaryland.edu/~jalemkul/).

So I did the following steps to prepare the PRMTOP and INPCRD files for the simulation.

1. Build a DMAC structure with Molefacture in VMD and optimize the structure with SQM.

2. Use Gaussian to optimize the geometry by

MP2/6-31G* opt.

3. Use Gaussian to calculate ESP by

HF/6-31G* SCF=Tight Geom=AllCheck Guess=Read Pop=MK IOp(6/33=2,6/41=10,6/42=17)

4. Use Antechamber to get PREP and FRCMOD file

antechamber -i esp.log -fi gout -o dma.prep -fo prepi -at gaff

antechamber -i dma.prep -f prepi -o dma.frcmod -p \$AMBERHOME/dat/leap/parm/gaff.dat

5. Edit dma.prep, change MOL to DMA. Edit NEWPDB.PDB, change MOL to DMA.

6. Use TLEAP to build a liquid and vacuum system.

#the script to build a liquid and vacuum system

source leaprc.gaff

check a

charge a

saveamberparm a dma_vac.prmtop dma_vac.inpcrd  #files for vacuum system

solvatebox a DMA 15.0

check a

charge a

saveamberparm a dma_liquid.prmtop dma_liquid.inpcrd    #files for liquid system

quit

The following was to do short MD simulation with NPT. I did 10ns simulation for liquid system and 100ns for vacuum system at 300K. And then restart a new system for both 1ns. And we would use the 1ns simulation results for heat of vaporization calculation.

The end of AMBER output were below.

liquid averaged result
117 molecule

NSTEP = 1000000 TIME(PS) = 11020.000 TEMP(K) = 300.22 PRESS = 19.0
Etot = 3017.4941 EKtot = 1570.5383 EPtot = 1446.9559
BOND = 579.1624 ANGLE = 766.5262 DIHED = 397.1170
1-4 NB = 298.1792 1-4 EEL = 7051.5797 VDWAALS = -942.9803
EELEC = -6702.6283 EHBOND = 0.0000 RESTRAINT = 0.0000
EKCMT = 104.7960 VIRIAL = 97.3657 VOLUME = 18100.2857
Density = 0.9352

#############################
vacuum averaged result
1 molecular

NSTEP = 1000000 TIME(PS) = 101000.000 TEMP(K) = 298.44 PRESS = 0.0
Etot = 37.5408 EKtot = 13.3439 EPtot = 24.1970
BOND = 5.0303 ANGLE = 6.4772 DIHED = 3.3979
1-4 NB = 2.5018 1-4 EEL = 60.2780 VDWAALS = 0.7379
EELEC = -54.2261 EHBOND = 0.0000 RESTRAINT = 0.0000

The following is the equation for the calculation.

dHvap=Hgas-Hliq=Ugas-Uliq+p(Vgas-Vliq)

Assuming Vgas >> Vliq

dHvap=Ugas – Uliq + pVgas = Ugas – Uliq + RT

U is the sum of kinetic and potential energies. So the heat of vaporization of DMA in this case is

Ugas = (13.3439 + 24.1970) / 1 = 37.5409 kcal/mol

Uliq = ( 1570.5383 + 1446.9559 ) / 117 = 25.7905 kcal/mol

RT = 8.314 * 300 / 1000 / 4.184 = 0.5961 kcal/mol

dHvap= 37.5409 – 25.7905 + 0.5961 = 12.3465 kcal/mol

The experimental data is 12.00 kcal/mol at 298K (James S. Chickos, “Enthalpies of Vaporization of Organic and Organometallic Compounds, 1880–2002,” Journal of Physical and Chemical Reference Data 32, no. 2 (2003): 519, doi:10.1063/1.1529214.)

The density is 0.9352g/ml and the experimental data is 0.937g/ml.

So the density and heat of vaporization both fit well with experimental data.

GAFF is good for DMA!!

Continue the tricky calculation of solvation free energies with TI : (

## 2 thoughts on “Calculating heat of vaporization of solvent with MD simulation.”

1. Hi Xia,

How you are simulating single molecule of DMA. Could you explain your system in details and md input parameter.

Thanks
Rahman

Like

• AMBER has a tutorial for DNA chains in vacuum, in the implicit solvent and in a TIP3P water box. I followed the one in vacuum for the simulation of single DMA.

Like