search

crest-lab / crest

CREST - A program for the automated exploration of low-energy molecular chemical space.
https://crest-lab.github.io/crest-docs/
GNU Lesser General Public License v3.0
198 stars 42 forks source link

Problem with pKa conformational average calculation #133

Closed ik160 closed 2 years ago

ik160 commented 2 years ago

Hello, I have installed CREST and xTB for pKa calculation. Following the method in the paper(https://pubs.acs.org/doi/10.1021/acs.jpca.1c03463) and its supporting information, I got the correct pKa(below).

       ==============================================
       |                                            |
       |                 C R E S T                  |
       |                                            |
       |  Conformer-Rotamer Ensemble Sampling Tool  |
       |          based on the GFN methods          |
       |             P.Pracht, S.Grimme             |
       |          Universitaet Bonn, MCTC           |
       ==============================================
       Version 2.12,   Thu 19. Mai 16:32:32 CEST 2022
  Using the xTB program. Compatible with xTB version 6.4.0

   Cite work conducted with this code as

   • P.Pracht, F.Bohle, S.Grimme, PCCP, 2020, 22, 7169-7192.
   • S.Grimme, JCTC, 2019, 15, 2847-2862.

   and for works involving QCG as

   • S.Spicher, C.Plett, P.Pracht, A.Hansen, S.Grimme,
     JCTC, 2022, 18 (5), 3174-3189.

   with help from:
   C.Bannwarth, F.Bohle, S.Ehlert, S.Grimme,
   C.Plett, P.Pracht, S.Spicher

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

 Command line input:
 > crest acid.xyz --pka 7

 -------------------------
 xTB Geometry Optimization
 -------------------------
 Geometry successfully optimized.

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++                             Calculation of the acid                                  ++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Optimizing acid geometry ...  done.

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++                             Calculation of the base                                  ++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Optimizing base geometry ...  done.

 Atom 3 (O ) identified as reactive center.
 Gsolv(xtb) A / B    : -2.094533565900000E-002 -0.105175485811000
 dE(xtb,uncorr.)     :  0.191486145461003
 dG(xtb,uncorr.)     :  0.180767612980003
 Energy correction for acid/base reaction:
 Ecorr =    0.2470493940 Eh

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++                                 pKa CALCULATION                                      ++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

 For the reaction AH + H₂O --->  A⁻ + H₃O⁺
 (Note: H₂O/H₃O⁺ is not included in ΔG)
 T     =          298.15 K
 G(AH) =    -21.63908339 Eh
 G(A⁻) =    -21.45831577 Eh
 ΔG       =      0.18076761 Eh,  113.43 kcal/mol
 **ΔG+Ecorr =      0.42781701 Eh,  268.46 kcal/mol**

 polynomial free energy relationship (FER):
  pKa   = c0 + c1*kdiss + c2*kdiss² + ... + c_n*kdiss^n
  with kdiss = ΔG(aq)/ln(10)RT and GFN2 parameters (c0-c3) fitted on PKA74
       ΔG(aq)=  268.46 kcal/mol
       c0 =  -1855.02527700
       c1 =     26.07598200
       c2 =     -0.12496355
       c3 =      0.00020571

    ___________________________________
   |                                   |
   | calculated pKa =             4.79 |
   |___________________________________|

 -----------------
 Wall Time Summary
 -----------------
--------------------
Overall wall time  : 0h : 0m : 0s

 CREST terminated normally.

But in the ensemble case(attached xyz files), CREST outputted too small pKa(below). It seems that ΔG does not include "Ecorr".

Could you help?


       ==============================================
       |                                            |
       |                 C R E S T                  |
       |                                            |
       |  Conformer-Rotamer Ensemble Sampling Tool  |
       |          based on the GFN methods          |
       |             P.Pracht, S.Grimme             |
       |          Universitaet Bonn, MCTC           |
       ==============================================
       Version 2.12,   Thu 19. Mai 16:32:32 CEST 2022
  Using the xTB program. Compatible with xTB version 6.4.0

   Cite work conducted with this code as

   • P.Pracht, F.Bohle, S.Grimme, PCCP, 2020, 22, 7169-7192.
   • S.Grimme, JCTC, 2019, 15, 2847-2862.

   and for works involving QCG as

   • S.Spicher, C.Plett, P.Pracht, A.Hansen, S.Grimme,
     JCTC, 2022, 18 (5), 3174-3189.

   with help from:
   C.Bannwarth, F.Bohle, S.Ehlert, S.Grimme,
   C.Plett, P.Pracht, S.Spicher

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

 Command line input:
 > crest --pka --pkaensemble acid_conformers.xyz base_conformers.xyz

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++                      Read ensemble data for pKa calculation                          ++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Note: In order to work properly the ensembles must be in the xyz
 format and contain the final free energies in solution (in Eh) as
 comment line for each structure.

 Ensemble read for acid: acid_conformers.xyz with 18 structures
 Ensemble read for base: base_conformers.xyz with 12 structures
 Calculating population average from read-in (free) energies ...
 |G(AH)| =    -21.69880855 Eh
 |G(A⁻)| =    -21.50825600 Eh

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++                                 pKa CALCULATION                                      ++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

 For the reaction AH + H₂O --->  A⁻ + H₃O⁺
 (Note: H₂O/H₃O⁺ is not included in ΔG)
 T     =          298.15 K
 G(AH) =    -21.69880855 Eh
 G(A⁻) =    -21.50825600 Eh
 ΔG       =      0.19055254 Eh,  119.57 kcal/mol

 polynomial free energy relationship (FER):
  pKa   = c0 + c1*kdiss + c2*kdiss² + ... + c_n*kdiss^n
  with kdiss = ΔG(aq)/ln(10)RT and GFN2 parameters (c0-c3) fitted on PKA74
       ΔG(aq)=  119.57 kcal/mol
       c0 =  -1855.02527700
       c1 =     26.07598200
       c2 =     -0.12496355
       c3 =      0.00020571

    ___________________________________
   |                                   |
   | calculated pKa =          -391.00 |
   |                                   |
   | min.pKa (conf) =          -427.15 |
   | max.pKa (conf) =          -352.69 |
   |___________________________________|

 -----------------
 Wall Time Summary
 -----------------
--------------------

Overall wall time  : 0h : 0m : 0s

 CREST terminated normally.

base_conformers.txt acid.txt acid_conformers.txt

pprcht commented 2 years ago

Hi, The ensemble function was intended for DFT calculations. It may be possible to enable the "Ecorr" term, but I will have to look that up over the weekend.

pprcht commented 2 years ago

Sorry for the delayed answer. I've looked it up and there is already a function for what you want to to. You have to execute the command crest acid.xyz --pkaprep acid_conformers.xyz base_conformers.xyz which will calculate the correction term for all base structures with reference to the lowest acid structure and generate two files G_acid_ensemble.xyz and G_base_ensemble.xyz. These two can then be used with the pkaensemble command as you did previously, i.e., crest --pka --pkaensemble G_acid_conformers.xyz G_base_conformers.xyz

Hope this helps. I will try to prepare some more detailed instructions for the online documentation soon.

ik160 commented 2 years ago

Thank you for your great help! This is what my want. With pkaprep, I did calculate min, max, average pKa. Thanks again.

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++                                 pKa CALCULATION                                      ++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

 For the reaction AH + H₂O --->  A⁻ + H₃O⁺
 (Note: H₂O/H₃O⁺ is not included in ΔG)
 T     =          298.15 K
 G(AH) =    -21.63824324 Eh
 G(A⁻) =    -21.20941820 Eh
 ΔG       =      0.42882504 Eh,  269.09 kcal/mol

 polynomial free energy relationship (FER):
  pKa   = c0 + c1*kdiss + c2*kdiss² + ... + c_n*kdiss^n
  with kdiss = ΔG(aq)/ln(10)RT and GFN2 parameters (c0-c3) fitted on PKA74
       ΔG(aq)=  269.09 kcal/mol
       c0 =  -1855.02527700
       c1 =     26.07598200
       c2 =     -0.12496355
       c3 =      0.00020571

    ___________________________________
   |                                   |
   | calculated pKa =             5.16 |
   |                                   |
   | min.pKa (conf) =             2.19 |
   | max.pKa (conf) =             7.32 |
   |___________________________________|

 -----------------
 Wall Time Summary
 -----------------
--------------------
Overall wall time  : 0h : 0m : 0s

 CREST terminated normally.