Here's using ccCA in NWChem on acetic acid/acetate and formic acid/formate.
More about how it works in another post.
Basically, the way I am using it the results are very, very poor with ccCA. All I can think is that I must be doing something wrong.
INPUT files
Acetic acid input:
Title "aceticacid" Start aceticacid echo charge 0 geometry autosym units angstrom C -0.312051 -1.36877 0.00000 H -0.929226 -1.55822 -0.878253 H -0.929226 -1.55822 0.878253 H 0.548700 -2.02934 0.00000 C 0.150590 0.0606620 0.00000 O -0.897092 0.922315 0.00000 H -0.521850 1.81528 0.00000 O 1.29371 0.435169 0.00000 end basis * library "cc-pvtz" end dft mult 1 direct noio XC b3lyp grid fine mulliken end driver default end ccca optimize end task ccca
Acetate input:
Title "acetate_ccca" Start acetate_ccca echo charge -1 geometry autosym units angstrom C -0.0311237 -1.36218 0.00000 H 0.501926 -1.73727 0.878691 H 0.501926 -1.73727 -0.878691 H -1.05131 -1.75101 0.00000 C -0.00500996 0.204086 0.00000 O 1.14247 0.706045 0.00000 O -1.12049 0.771493 0.00000 end basis * library "cc-pvtz" end dft mult 1 direct noio XC b3lyp grid fine mulliken end driver default end ccca optimize end task ccca
Formic acid input:
Title "formicacid" Start formicacid echo charge 0 geometry autosym units angstrom C 0.410955 -0.132154 0.00000 H 1.50430 -0.0475164 0.00000 O -0.134104 1.09718 0.00000 H -1.09846 0.988665 0.00000 O -0.203188 -1.15938 0.00000 end basic * library "cc-pvtz" end dft mult 1 direct noio XC b3lyp grid fine mulliken end driver end ccca optimize end task ccca
Formate input:
Title "formate" Start formate echo charge -1 geometry autosym units angstrom C 0.00000 0.00000 0.329396 H 0.00000 0.00000 1.47310 O -1.13532 0.00000 -0.189103 O 1.13532 0.00000 -0.189103 end basis "ao basis" spherical print * library "cc-pvtz" END dft mult 1 direct noio XC b3lyp grid fine mulliken end driver end ccca optimize end task ccca
OUTPUT
Acetic acid
Temperature = 298.15K frequency scaling parameter = 0.9889 Zero-Point correction to Energy = 38.155 kcal/mol ( 0.060805 au) Thermal correction to Energy = 41.060 kcal/mol ( 0.065434 au) Thermal correction to Enthalpy = 41.653 kcal/mol ( 0.066378 au) Total Entropy = 69.467 cal/mol-K - Translational = 38.179 cal/mol-K (mol. weight = 60.0211) - Rotational = 23.830 cal/mol-K (symmetry # = 1) - Vibrational = 7.458 cal/mol-K Cv (constant volume heat capacity) = 14.439 cal/mol-K - Translational = 2.979 cal/mol-K - Rotational = 2.979 cal/mol-K - Vibrational = 8.480 cal/mol-K ccCA: calculations done, now printing results ccCA-P reference energy = -228.82035086993045 ccCA-S3 reference energy = -228.82800135561300 ccCA-S4 reference energy = -228.82030423449530 ccCA-PS3 reference energy= -228.82417611277174 DK correction = -0.13322049012506909 CCSD(T) correction = -4.8762979862686961E-002 CV correction = -0.20936881324035994 --------------------------- Total ccCA-P energy = -229.21170315315857 Total ccCA-S3 energy = -229.21935363884111 Total ccCA-S4 energy = -229.21165651772341 Total ccCA-PS3 energy = -229.21552839599985 Thermochemistry available: ZPE = 6.0851792771826778E-002 ccCA-P E+ZPE = -229.15085136038675 ccCA-S3 E+ZPE = -229.15850184606930 ccCA-S4 E+ZPE = -229.15080472495160 ccCA-PS3 E+ZPE = -229.15467660322804 Wrote ccCA-P energy to the RTDB Leaving ccCA module... Task times cpu: 5565.6s wall: 5650.7s
Acetate
Temperature = 298.15K frequency scaling parameter = 0.9889 Zero-Point correction to Energy = 29.591 kcal/mol ( 0.047157 au) Thermal correction to Energy = 31.853 kcal/mol ( 0.050762 au) Thermal correction to Enthalpy = 32.446 kcal/mol ( 0.051706 au) Total Entropy = 64.067 cal/mol-K - Translational = 38.129 cal/mol-K (mol. weight = 59.0133) - Rotational = 23.739 cal/mol-K (symmetry # = 1) - Vibrational = 2.199 cal/mol-K Cv (constant volume heat capacity) = 11.235 cal/mol-K - Translational = 2.979 cal/mol-K - Rotational = 2.979 cal/mol-K - Vibrational = 5.276 cal/mol-K ccCA: calculations done, now printing results ccCA-P reference energy = -228.25857936176124 ccCA-S3 reference energy = -228.26625083689740 ccCA-S4 reference energy = -228.25849407721080 ccCA-PS3 reference energy= -228.26241509932930 DK correction = -0.13318127658752132 CCSD(T) correction = -4.4728554700242285E-002 CV correction = -0.20921905251765338 --------------------------- Total ccCA-P energy = -228.64570824556665 Total ccCA-S3 energy = -228.65337972070282 Total ccCA-S4 energy = -228.64562296101622 Total ccCA-PS3 energy = -228.64954398313472 Thermochemistry available: ZPE = 4.7193435242008613E-002 ccCA-P E+ZPE = -228.59851481032464 ccCA-S3 E+ZPE = -228.60618628546081 ccCA-S4 E+ZPE = -228.59842952577421 ccCA-PS3 E+ZPE = -228.60235054789271 Wrote ccCA-P energy to the RTDB Leaving ccCA module... Task times cpu: 3859.1s wall: 3910.2s
Formic acid
Temperature = 298.15K frequency scaling parameter = 0.9889 Zero-Point correction to Energy = 20.909 kcal/mol ( 0.033320 au) Thermal correction to Energy = 22.902 kcal/mol ( 0.036497 au) Thermal correction to Enthalpy = 23.495 kcal/mol ( 0.037441 au) Total Entropy = 59.329 cal/mol-K - Translational = 37.387 cal/mol-K (mol. weight = 46.0055) - Rotational = 21.008 cal/mol-K (symmetry # = 1) - Vibrational = 0.934 cal/mol-K Cv (constant volume heat capacity) = 8.703 cal/mol-K - Translational = 2.979 cal/mol-K - Rotational = 2.979 cal/mol-K - Vibrational = 2.744 cal/mol-K ccCA: calculations done, now printing results ccCA-P reference energy = -189.56748775122853 ccCA-S3 reference energy = -189.57364633780318 ccCA-S4 reference energy = -189.56733835209894 ccCA-PS3 reference energy= -189.57056704451585 DK correction = -0.11856238070660652 CCSD(T) correction = -3.0831132609506540E-002 CV correction = -0.16057296161548607 --------------------------- Total ccCA-P energy = -189.87745422616013 Total ccCA-S3 energy = -189.88361281273478 Total ccCA-S4 energy = -189.87730482703054 Total ccCA-PS3 energy = -189.88053351944745 Thermochemistry available: ZPE = 3.3346398704552728E-002 ccCA-P E+ZPE = -189.84410782745556 ccCA-S3 E+ZPE = -189.85026641403022 ccCA-S4 E+ZPE = -189.84395842832598 ccCA-PS3 E+ZPE = -189.84718712074289 Wrote ccCA-P energy to the RTDB Leaving ccCA module... Task times cpu: 1369.3s wall: 1407.5s
Formate
Temperature = 298.15K frequency scaling parameter = 0.9889 Zero-Point correction to Energy = 12.385 kcal/mol ( 0.019737 au) Thermal correction to Energy = 14.252 kcal/mol ( 0.022713 au) Thermal correction to Enthalpy = 14.845 kcal/mol ( 0.023656 au) Total Entropy = 56.927 cal/mol-K - Translational = 37.321 cal/mol-K (mol. weight = 44.9976) - Rotational = 19.229 cal/mol-K (symmetry # = 2) - Vibrational = 0.377 cal/mol-K Cv (constant volume heat capacity) = 7.310 cal/mol-K - Translational = 2.979 cal/mol-K - Rotational = 2.979 cal/mol-K - Vibrational = 1.351 cal/mol-K ccCA: calculations done, now printing results ccCA-P reference energy = -189.01189831122844 ccCA-S3 reference energy = -189.01808726799254 ccCA-S4 reference energy = -189.01171261173857 ccCA-PS3 reference energy= -189.01499278961049 DK correction = -0.11851294005154500 CCSD(T) correction = -2.6430727035545942E-002 CV correction = -0.16057463040127118 --------------------------- Total ccCA-P energy = -189.31741660871680 Total ccCA-S3 energy = -189.32360556548090 Total ccCA-S4 energy = -189.31723090922694 Total ccCA-PS3 energy = -189.32051108709885 Thermochemistry available: ZPE = 1.9751889903778755E-002 ccCA-P E+ZPE = -189.29766471881302
ccCA-S3 E+ZPE = -189.30385367557713 ccCA-S4 E+ZPE = -189.29747901932316 ccCA-PS3 E+ZPE = -189.30075919719508 Wrote ccCA-P energy to the RTDB Leaving ccCA module...
Solvation energy
Solvation energy may seem easy to calculate, but difficult to calculate accurately using implicit methods, in particular for ions. I used the optimized structures from above, and then did a single-point COSMO (rsolv 0. Not ideal) at RDFT(b3lyp)/cc-pVTZ/
Acetic acid: -8.59 kcal/mol
Acetate: -72.33 kcal/mol
Formic acid: -8.90
Formate: -72.59 kcal/mol
H+: -624.61 kcal/mol (lit. value)
Free energies:
G: ccCA-P+(Hcorr-T*Scorr)
G(acetic acid): -229.21170315315857*627.503+(41.653-298.15*69.467/1000)-8.59
G(acetate): -228.64570824556665*627.503+(32.446-298.15*64.067/1000)-72.33
G(formic acid): -189.87745422616013*627.503+(23.495-298.15*59.329/1000)-8.90
G(formate): -189.31741660871680*627.503+(14.845-298.15*56.927/1000)-72.59
G(H+): -6.28 kcal/mol (lit. value) -264.61 (lit. value)=-270.89 kcal/mol
Results: Direct approach:
Don't forget to account for the standard state. R=1.9858775(34)×10−3 kcal/(K.mol)
DG(acetic/acetate)=G(acetate)+G(H+)-G(acetic)+RT*ln(1/24.46)=
(-228.64570824556665*627.503+(32.446-298.15*64.067/1000)-72.33-270.89)-(-229.21170315315857*627.503+(41.653-298.15*69.467/1000)-8.59)+1.9858775*298.15*log(1/24.46)/1000=
11.0435795929699 kcal/mol=46.2063370169861 kJ/mol =>
pKa=DG*log10(e)/RT=
11.0435795929699*log10(e)/(1.9858775*10**(-3)*298.15)
=8.1 (which is very bad -- it should be close to 4.75)
DG(formic/formate)=G(formate)+G(H+)-G(formic)-RT*ln(1/24.46)=
(-189.31741660871680*627.503+(14.845-298.15*56.927/1000)-72.59-270.89)-(-189.87745422616013*627.503+(23.495-298.15*59.329/1000)-8.90)+1.9858775*298.15*log(1/24.46)/1000=
7.01850845285312 kcal/mol=29.3654393667375 kJ/mol
pKa=5.1 (which is quite bad -- it should be close to 3.75)
Results: Isodesmic approach
From an older post:
"Assuming that we know that formic acid has a pKa of 3.75, then DG_solution=pKa*RT/log(e)=3.75*8.314*298.15/log10(e)/1000=21.404 kJ/mol. The reverse reaction is -21.404 kJ/mol."
That's about 5.116 kcal/mol
DG(acetate)+DG(formic)-(DG(acetic)+DG(formate))+DG(ref)=
((-228.64570824556665*627.503+(32.446-298.15*64.067/1000)-72.33)+(-189.87745422616013*627.503+(23.495-298.15*59.329/1000)-8.90))-((-229.21170315315857*627.503+(41.653-298.15*69.467/1000)-8.59)+(-189.31741660871680*627.503+(14.845-298.15*56.927/1000)-72.59))+5.116=
4.02507114014588 kcal/mol+5.116 kcal/mol =9.14107114014588 kcal/mol <=> pKa= 6.7
Which is better, but still not as good as here.
Using the E+zpe energies doesn't help much:
((-228.59851481032464*627.503+(32.446-298.15*64.067/1000)-72.33)+(-189.84410782745556*627.503+(23.495-298.15*59.329/1000)-8.90))-((-229.15085136038675*627.503+(41.653-298.15*69.467/1000)-8.59)+(-189.29766471881302*627.503+(14.845-298.15*56.927/1000)-72.59))+5.116=
9.10100587110175 kcal/mol <=> pKa=6.68
I really have no idea why the results are so bad when I had reasonable results with DFT/b3lyp/6-31++G** which should be worse than the E(CBS)+E(CC)+E(CV)+E(DK) approach for calculating electronic energies.
Solvation energies are a bit different and could explain some of the difference. Using the solvation energies from here I got:
((-228.64570824556665*627.503+(32.446-298.15*64.067/1000)-73.23)+(-189.87745422616013*627.503+(23.495-298.15*59.329/1000)-9.99))-((-229.21170315315857*627.503+(41.653-298.15*69.467/1000)-9.32)+(-189.31741660871680*627.503+(14.845-298.15*56.927/1000)-72.47))+5.116=
7.76107114008302 kcal/mol <=> pKa=5.69. Not 4.75, but closer.
Using rsolv 1.3 I get
((-228.64570824556665*627.503+(32.446-298.15*64.067/1000)-71.09)+(-189.87745422616013*627.503+(23.495-298.15*59.329/1000)-6.37))-((-229.21170315315857*627.503+(41.653-298.15*69.467/1000)-6.53)+(-189.31741660871680*627.503+(14.845-298.15*56.927/1000)-71.90))+5.116=
10.1610711401063 kcal/mol which is bad.
10.1610711401063 kcal/mol which is bad.
More thinking.
This paper says that the gas phase free energy for the deprotonation of acetic acid should be 341.1 kcal/mol
(-228.64570824556665*627.503+(32.446-298.15*64.067/1000)-6.82)-(-229.21170315315857*627.503+(41.653-298.15*69.467/1000))=340.75 kca/mol
We're within 1 kcal/mol
The same paper states 338.5 kcal/mol for formic acid:
(-189.31741660871680*627.503+(14.845-298.15*56.927/1000)-6.82)-(-189.87745422616013*627.503+(23.495-298.15*59.329/1000))=336.67 kca/mol
For our direct solution phase pKa calculation formic acid was off by about 1.9 kcal/mol which is similar to the error here.
