Unusually High Electrostatic Energy for Ionic Systems w/ SAPT2+3?

Hi all,

I wanted to investigate the interactions between charged groups with SAPT. I followed the input from the manual. The output showed very large values for electrostatics (-118.19630204 kcal/mol). I’m not too sure if this makes physical sense. I minimised my structures first with MP2 and a double zeta basis set. Can anyone shed some light on why I am getting these values?

Thanks in advance,
Eric

Input

molecule arg_carb {
     1 1
     N         2.871509         0.425563         -0.032645
     C         1.544924         0.290062         -0.085864
     N         0.971332         -0.894907         0.114337
     N         0.766447         1.346429         -0.360572
     C         3.781042         -0.679543         0.207392
     H         3.252449         1.359164         -0.140984
     H         3.701853         -1.438759         -0.580422
     H         3.592498         -1.147313         1.181541
     H         4.798827         -0.286649         0.204851
     H         1.18766         2.268526         -0.326197
     H         -0.236255         1.274701         -0.118454
     H         1.538184         -1.72862         0.208982
     H         -0.040825         -0.9988         -0.075165
     --
     -1 1
     O         -1.814201         -1.04752         -0.474678
     O         -1.962211         0.994214         0.431762
     C         -2.478162         -0.08179         0.006249
     C         -3.977669         -0.229866         0.113582
     H         -4.361611         -0.983381         -0.582074
     H         -4.482006         0.727497         -0.05895
     H         -4.225293         -0.552688         1.134205


     units angstrom
}

set globals {
    basis          aug-cc-pvdz
    guess          sad
    scf_type       df
}

set sapt {
    print          1
    nat_orbs_t2    true
    freeze_core    true
    aio_cphf      true
    aio_df_ints   true
}

energy('sapt2+3')

Output

  Special recipe for scaled SAPT0 (see Manual):
    Electrostatics sSAPT0        -188.35778907 [mEh]    -118.19630204 [kcal/mol]    -494.53337521 [kJ/mol]
    Exchange sSAPT0                44.07183788 [mEh]      27.65549695 [kcal/mol]     115.71061035 [kJ/mol]
    Induction sSAPT0              -36.27183851 [mEh]     -22.76092325 [kcal/mol]     -95.23171200 [kJ/mol]
    Dispersion sSAPT0             -11.32677616 [mEh]      -7.10765964 [kcal/mol]     -29.73845080 [kJ/mol]
  Total sSAPT0                   -191.88456586 [mEh]    -120.40938798 [kcal/mol]    -503.79292767 [kJ/mol]

You are surprised a complex of a cation and an anion has very high electrostatics?
Check with a supermolecular method if the interaction energy is OK

1 Like

Hi Hokru,

Thanks for your fast response. No that isn’t what has surprised me. The magnitude of the interaction seems unexpectedly large to me, but I think that is because I misunderstand the SAPT method…

I have compared the results to a supermolecular method and there is good agreement between energies.

DLPNO-CCSD(T)/Def2-TZVPP    SAPT2+3

A -244.6823392 -243.5952203
B -228.1812019 -227.272618
AB -473.0537938 -471.0490881

So what does the Electrostatics sSAPT0 energy actually refer to? Apologies for my misunderstanding.

No idea what these numbers are supposed to tell me.
Roughly, the SAPT electrostatics is the Coulomb energy (think classical charge-cloud electrostatics) between the unperturbed monomer densities. The references in the psi4 manual go into specific details.

sSAPT0 uses empirical scaling factors to improve upon normal SAPT0.

Thanks for your help!

To add to this, it looks like you are calculating the interaction energy between a guanidinium group from an arginine side chain and a carboxylate group from a glutamate or aspartate side chain. If so, the electrostatic interaction between these charged groups is likely influenced by the electric field from the rest of the protein and/or solvent. Your arginine zeta carbon and the carboxylate group are about 3.5 Å apart. The Coulomb energy between two classical point charges at this distance is -95 kcal/mol, so the strength of the interaction you are seeing is on the right order of magnitude for a vacuum calculation.

However, this vacuum calculation likely doesn’t represent the strength of the actual interaction in the protein. Estimates for the dielectric constant inside of proteins vary widely and are strongly dependent on what side chains are nearby, but Figure 7 of this 2013 paper estimates that the dielectric constant near charged amino acids is ~20. This means that the strength of the electrostatic interaction for our classical point charges drops from 95 kcal/mol to ~5 kcal/mol.

If you want to account for this kind of dielectric screening in your QM energy, you’ll have to include many more atoms of the protein in your SAPT calculation, which can get expensive fast. You may have to use a less rigorous version of SAPT (e.g. sSAPT0 instead of SAPT2+3) depending on your system size and available resources. I think it is better to use cheaper QM with a good representation of your system than to use better QM with a poor representation of your system, but there are others who disagree.

Of course, if you aren’t working with a protein or if you want to look at vacuum energies, then don’t worry about what I said.

1 Like