QCG G_solv issue with monatomic solutes (noble gases)

[alan-arnold]

Hi. I've been playing with QCG to look at the aquation of the noble gases, as a prelude to QCG simulations of noble gas binding to cucurbiturils in the presence of explicit waters.

I'm using the precompiled Linux versions of xtb 6.5.1. with crest 2.12 and xtbiff 1.1

Below are two snippets from the end of the output files from QCG runs of Kr in clusters of 50 and 100 waters, respectively.

The --nopreopt option is necessary to prevent an "Initial geometry optimization failed!" error, also mentioned in issue #149

Notice the translational entropy of the Kr in the two runs (-254.63 and -509.37 cal/mol/K) under the "Solute Gas properties" heading. Why should STRA(Kr) be twice as large with 100 waters compared to 50? Naively, I would have expected the solute-only entropy components to be independent of the number of waters in the simulation.

The Sackur-Tetrode equation gives 39.1 cal/mol/K for the molar translational entropy of Kr80 at 298K and 1atm but I'm not sure whether this translates (if you pardon the pun) into a simulation in water :)

So, I worry whether G_solv is calculated correctly here? The experimental value for the solvation free energy of Kr in water, mentioned in https://pure.rug.nl/ws/portalfiles/portal/64489770/1.451846.pdf is 6.95 kJ/mol (1.66 kcal/mol), very much different to these calculations.

As an aside, there's a line in the output (not shown below) which says "Single point computation with GBSA model", but these runs are using ALPB. Is this an oversight?

crest Kr.xyz --T 24 --qcg h2o.xyz --gsolv --nsolv 50 --alpb h2o --mdlen 10 -keepdir --nopreopt > Kr_crest_10_qcg50_1.out

..............................................................	quantum cluster growth: ESOLV
9.70 kcal/mol

.........................................

......................................... | Frequency evaluation | .........................................

SOLUTE MOLECULE Starting reoptimizations + Frequency computation of ensemble 1 jobs to do.

done. SOLUTE CLUSTER Starting reoptimizations + Frequency computation of ensemble 4 jobs to do.

done. SOLVENT CLUSTER Starting reoptimizations + Frequency computation of ensemble 4 jobs to do.

done.

Solute Gas properties H(T) SVIB SROT STRA G(T) [kcal/mol] [ cal/mol/K ] [kcal/mol] ........................................................ 0.00 0.00 0.00 -254.63 0.00

Solute cluster properties H(T) SVIB SROT STRA G(T) [kcal/mol] [ cal/mol/K ] [kcal/mol] ........................................................ 1 833.01 475.65 40.24 46.51 665.34 2 833.94 474.64 40.24 46.51 666.55 3 833.73 471.66 40.22 46.51 667.24 4 833.72 479.92 40.19 46.51 664.78

Solvent cluster properties H(T) SVIB SROT STRA G(T) [kcal/mol] [ cal/mol/K ] [kcal/mol] ........................................................ 1 831.49 458.26 40.16 46.25 669.10 2 832.40 461.03 40.05 46.25 669.21 3 832.01 453.32 40.02 46.25 671.13 4 832.10 467.79 40.10 46.25 666.89

_____ Evaluation ____

........................................................ Gsolv and Hsolv ref. state: [1 M gas/solution] G_solv (incl.RRHO) = -69.75 kcal/mol H_solv (incl.RRHO) = 11.80 kcal/mol ........................................................

........................................................ Gsolv and Hsolv ref. state: [1 M gas/solution] G_solv (incl.RRHO) = -71.65 kcal/mol H_solv (incl.RRHO) = 9.91 kcal/mol ........................................................

........................................................ Solvation free energies with scaled translational and rotational degrees of freedom: Gsolv (scaling)

 0.58 (0.05)    <<
-3.22 (0.10)    <<
-7.03 (0.15)    <<

-10.83 (0.20) << -14.63 (0.25) << -18.43 (0.30) << -22.23 (0.35) << -26.03 (0.40) << -29.83 (0.45) << -33.63 (0.50) << -37.44 (0.55) << -41.24 (0.60) << -45.04 (0.65) << -48.84 (0.70) << -52.64 (0.75) << -56.44 (0.80) << -60.24 (0.85) << -64.04 (0.90) << -67.84 (0.95) << -71.65 (1.00) << ........................................................

..................................................	Gsolv with SCALED RRHO contributions: 0.75		[1 bar gas/ 1 M solution]
G_solv (incl.RRHO)+dV(T)= -52.64 kcal/mol

..................................................

AND

crest Kr.xyz --T 24 --qcg h2o.xyz --gsolv --nsolv 100 --alpb h2o --mdlen 10 -keepdir --nopreopt > Kr_crest_10_qcg100_1.out &

.........................................	quantum cluster growth: ESOLV
6.51 kcal/mol

.........................................

......................................... | Frequency evaluation | .........................................

SOLUTE MOLECULE Starting reoptimizations + Frequency computation of ensemble 1 jobs to do.

done. SOLUTE CLUSTER Starting reoptimizations + Frequency computation of ensemble 4 jobs to do.

done. SOLVENT CLUSTER Starting reoptimizations + Frequency computation of ensemble 4 jobs to do.

done.

Solute Gas properties H(T) SVIB SROT STRA G(T) [kcal/mol] [ cal/mol/K ] [kcal/mol] ........................................................ 0.00 0.00 0.00 -509.37 0.00

Solute cluster properties H(T) SVIB SROT STRA G(T) [kcal/mol] [ cal/mol/K ] [kcal/mol] ........................................................ 1 1669.93 937.97 43.52 48.45 1362.86 2 1669.75 942.42 43.50 48.45 1361.35 3 1670.19 953.07 43.46 48.45 1358.62 4 1669.62 947.90 43.45 48.45 1359.61

Solvent cluster properties H(T) SVIB SROT STRA G(T) [kcal/mol] [ cal/mol/K ] [kcal/mol] ........................................................ 1 1668.12 925.65 43.48 48.31 1364.77 2 1668.23 930.05 43.46 48.31 1363.58 3 1668.74 941.00 43.41 48.31 1360.83 4 1668.16 935.61 43.41 48.31 1361.86

_____ Evaluation ____

........................................................ Gsolv and Hsolv ref. state: [1 M gas/solution] G_solv (incl.RRHO) = -147.42 kcal/mol H_solv (incl.RRHO) = 8.50 kcal/mol ........................................................

........................................................ Gsolv and Hsolv ref. state: [1 M gas/solution] G_solv (incl.RRHO) = -149.32 kcal/mol H_solv (incl.RRHO) = 6.60 kcal/mol ........................................................

........................................................ Solvation free energies with scaled translational and rotational degrees of freedom: Gsolv (scaling)

-4.99 (0.05)    <<

-12.58 (0.10) << -20.18 (0.15) << -27.78 (0.20) << -35.37 (0.25) << -42.97 (0.30) << -50.57 (0.35) << -58.16 (0.40) << -65.76 (0.45) << -73.35 (0.50) << -80.95 (0.55) << -88.55 (0.60) << -96.14 (0.65) << -103.74 (0.70) << -111.34 (0.75) << -118.93 (0.80) << -126.53 (0.85) << -134.12 (0.90) << -141.72 (0.95) << -149.32 (1.00) << ........................................................

..................................................	Gsolv with SCALED RRHO contributions: 0.75		[1 bar gas/ 1 M solution]
G_solv (incl.RRHO)+dV(T)= -111.34 kcal/mol

..................................................

[cplett]

Hi, The solute properties like the translational entropy should indeed be independent of the solvent cluster as it is computed in the gas phase. When I tested this with the newest versions compiled from the source code, the translational entropy ends up always being 39.17 cal/mol/K for the solute in the gas phase. This also leads to a positive solvation free energy for Kr (tested with 40 water molecules). I'm wondering whether this is a problem related to the different versions. It might be helpful if you could check the molecule in the qcg_tmp/tmp_freq/tmp_gas1 directory. There, the gas-phase frequency calculations should be performed and only the solute coordinates should be present.

Regarding the GBSA/ALPB printout: actually, if the --alpb flag is used, ALPB is also employed. It is just a wrong printout. Thank you for pointing this out.

[alan-arnold]

Thanks for checking this. It looks like the problem is in my version of xtb

• xtb version 6.5.1 (579679a) compiled by 'ehlert@majestix' on 2022-07-11

xtb_freq.out in qcg_tmp/tmp_freq/tmp_gas1 didn't converge:

[ERROR] Program stopped due to fatal error -2- Geometry optimization failed -1- xtb_geoopt: Geometry optimization did not converge

Any chance of getting hold of the most recent x64 Linux binaries? Unfortunately I don't have access to Intel compilers for fast run-times. A 50-waters Kr run takes 3-4 hr on Ubuntu 22.04 under Parallels 18 on my old 12-core Mac Pro so every second counts :) In my hands, conda-forge makes are much slower on macOS.

[cplett]

This seems to be related to macOS as xtb 6.5.1 on ubtunut works for the example. I will check on this and also include some checks for single atoms in the QCG. Could you please provide me the solute.coord file in the qcg_tmp/tmp_freq/tmp_gas1 directory, the Kr.xyz file, and maybe also the complete xtb_freq.out file? This would help me to spot the problem.

[alan-arnold]

Thanks for adding the single-atom exception checking into the QCG code. I've added the files you requested below, after appending .txt to their filenames because this editor doesn't support .coord, .xyz or .out :( solute.coord copy.txt xtb_freq copy.out.txt Kr copy.xyz.txt

[cplett]

Thank you for the input. Somehow, I couldn't reproduce the error yet. Could you please try, if a xtb Kr.xyz --freq --gfn2 does work for you?

[alan-arnold]

xtb Kr.xyz --freq --gfn2 seems to work OK: Kr_freq_gfn2_1.out.txt

[cplett]

I noticed that I made a mistake. The --freq is actually not a valid option. I meant to write --hess. Sorry for that. However, I guess the latest patch of the source code #159 might solve your problem. If you have the gfortran compiler, you can try to compile the current source code of crest with gfortran and use this version.

[alan-arnold]

xtb Kr.xyz --hess --gfn2 seems to work OK and S_trans is calculated correctly to be 39.17 cal/mol/K. Kr_hess_gfn2_1.out.txt

I'll try to compile the latest version. ... Success! A Meson build of the latest xtb and crest codes, using gfortran and using openblas produces the correct S_trans and G_solv that looks reasonable. Kr_crest_10_qcg20_2.out.txt