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
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
a = loadpdb NEWPDB.PDB
saveamberparm a dma_vac.prmtop dma_vac.inpcrd #files for vacuum system
solvatebox a DMA 15.0
saveamberparm a dma_liquid.prmtop dma_liquid.inpcrd #files for liquid system
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
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
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.
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 : (