alan-arnold
commented
1 year ago Well this is weird, and more than a little embarrassing. I really should look more closely at the beginnng of the output files, instead of jumping to the business end.
It turns out that the 1st two runs of the 1st set of three ie cluster growth using --xtbiff, and the 1st restart, which were both done in a particular working directory, correctly have charge =0 for Q6.
After these 2 runs, I changed to a new working directory, on a NVMe SSD drive, to see if run times might improve. (They didn't :)
Although they used exactly the same Q6.xyz and h2o.xyz structure files as these 1st two runs, these later 10 runs incorrectly use charge=1, which apparently was read from a .CHRG file in this new working directory. I cannot recall having created this file, but presumably I must have at some stage for another previous project. In my defence, a file starting with a period is hidden in MacOS (under which Ubuntu is running in my case), making this a bit tricky to find!!
As an aside, in all 12 runs, uhf is reported as zero, which is not consistent with charge=0 for two runs and charge=1 for the others, so it looks like neither crest nor xtb checks for this inconsistency. Perhaps something for a future feature request.
And I'm not sure if this charge difference accounts for the inconsistent Esolv values. But, just in case, I've deleted the offending hidden .CHRG file and will do another set of 3 runs, using --xtbiff to create a cluster (much quicker than aISS), followed by 2 restarts, and see if that gives Esolv consistent with Set 1.
I'll report back ...
cplett
github-actions[bot]
Hello again Christoph. You might recall that I've been using crest/xtb to investigate at the solvation energies and host-guest binding properties of a range of aquated cucurbiturils.
In the particular series of calculations described below, I've been trying to calculate the solvation energy of cucurbit[6]uril (Q6), with limited success.
All runs were performed with the crest 2.12 and xtb 6.5.1 codes downloaded on 24 Dec 2022 and compiled with v2021.6.0 ifort/icc from the 2022.2.0 Intel oneAPI toolkits, and the precompiled version of xtbiff 1.1, under Ubuntu 22.04. This is the same setup that you helped me get working in https://github.com/crest-lab/crest/issues/163
I performed four sets of three 5-ps QCG gsolv calculations. In each set, the 1st run grows a 100-water cluster around Q6 then the 2nd and 3rd runs were restarts, using the same grow cluster as the 1st run. In all runs, the --nofix option was used to permit relaxation of the Q6 solute in the solvation process.
For each of these 4 sets of 3 calculations, 4 final clusters are produced by default, from which Esolv is calculated, followed by frequency calulations to evaluate Gsolv and Hsolv of the Q6 solute
The data below are the means and standard deviations of the thermochemical properties reported by these 12 frequency calculations for the gas-phase solute, and 4 solute- and 4 solvent-clusters, of the 3 replicate runs for each of the 4 grow sets
In the xTB-IFF-grown cluster, 4 waters are placed inside the Q6 cavity. In each of the 3 aISS-grown clusters, 2 additional waters are in the cavity. In all clusters, most of the remaining waters are tightly around the two 6-oxygen portals, with the remainder spread over the outside of the more hydrophobic solute walls. Although the 6-O portals are completely solvated, the exterior walls are not quite completely covered with waters.
Despite these reservations, the intra-set and inter-set reproducibilty of these 12 frequency calculations of 100-water solvations is remarkable.
Because the Q6 solute is quite rigid, a single low-energy conformer results in all calculations (a 2nd slightly-higher conformer is observed in one of the 3 runs of the 1st set), and the thermochemical properties of the solute- and solvent-clusters for the 4 sets seem to be essentially statistically identical.
What is surprising is that the 4th set, which does NOT explicitly call for alpb solvation, is identical to the other three, which do explicitly call for alpb in the command line.
1st set of 3 calculations (xTB-IFF used to grow a 1st cluster, with alpb implicit solvation):
crest Q6.xyz --T 24 --qcg h2o.xyz --gsolv --nsolv 100 --alpb h2o --mdlen 5 --xtbiff --nofix -keepdirfollowed by 2 runs ofcrest Q6.xyz --T 24 --qcg h2o.xyz --gsolv --alpb h2o --mdlen 5 --xtbiff --nofix -keepdir2nd set of 3 calculations (aISS used to grow a 2nd cluster, with alpb implicit solvation):
crest Q6.xyz --T 24 --qcg h2o.xyz --gsolv --nsolv 100 --alpb h2o --mdlen 5 --nofix -keepdirfollowed by 2 runs ofcrest Q6.xyz --T 24 --qcg h2o.xyz --gsolv --alpb h2o --mdlen 5 --nofix -keepdir3rd set of 3 calculations (aISS again used to grow a 3rd cluster, again with alpb implicit solvation):
crest Q6.xyz --T 24 --qcg h2o.xyz --gsolv --nsolv 100 --alpb h2o --mdlen 5 --nofix -keepdirfollowed by 2 runs ofcrest Q6.xyz --T 24 --qcg h2o.xyz --gsolv --alpb h2o --mdlen 5 --nofix -keepdir4th set of 3 calculations (aISS used to grow a 4th cluster, but NOT with alpb implicit solvation):
crest Q6.xyz --T 24 --qcg h2o.xyz --gsolv h2o --nsolv 100 --mdlen 5 --nofix -keepdirfollowed by 2 runs ofcrest Q6.xyz --T 24 --qcg h2o.xyz --gsolv h2o --mdlen 5 --nofix -keepdirThat's the good news.
The not-so-good news is that Esolv (and consequently, Gsolv and Hsolv) is not reproducible at all. Indeed, Esolv has a spread of 88 kcal/mol across the 12 runs!!
1st set of 3 calculations (xTB-IFF used to grow a 1st cluster, with alpb implicit solvation): Esolv = -144.72(25), -160.39(25), -216.44(50) kcal/mol
2nd set of 3 calculations (aISS used to grow a 2nd cluster, with alpb implicit solvation): Esolv = -131.28(50), -143.34(25), -138.01(25) kcal/mol
3rd set of 3 calculations (aISS again used to grow a 3rd cluster, again with alpb implicit solvation): Esolv = -136.87(76), -128.45(75), -131.62(25) kcal/mol
4th set of 3 calculations (aISS used to grow a 4th cluster, but NOT with alpb implicit solvation): Esolv = -169.76(26), -169.64(50), -187.22(50) kcal/mol
The numbers in parentheses are the number of structures within a 3.00 kcal/mol window in the final crest_rotamers_5.xyz
So, I have a few questions that I hope you can help me with:
Would you expect the Esolv results to be more reproducible, given the highly consistent values of the thermochemical properties found in the frequency calculations and the rigidity of the Q6 solute? Don't the frequency calculations suggest that the 12x4 100-water+Q6 clusters have very similar structures, and therefore should have similar Esolv's??
Why is the number of structures within a 3.00 kcal/mol window in the final crest_rotamers_5.xyz ensembles always close to a multiple of 25? Given the enormous conformational space of these 100-water clusters I would have expected to see much more diversity in the number of structures within the 3kcal/mol window - perhaps always more than 25, by why not 37 or 42 or 131 ...?
Why are the thermochemical properties from the frequency calculations for the 4th set, which does not call for alpb solvation in the command line, so similar to the other 3 sets which do use alpb? Is alpb used by default, even if you don't call for it?
When xtbiff was used to grow a 100-water cluster (Set 1) the grow time was 5h 42m. When aISS was used (Sets 2, 3 and 4) the grow times were 25h 22m, 24h 44m and 26h 53m, respectively. Would you expect the new aISS algorithm to be so much slower than xtbiff to grow these clusters??
As usual, I've used -keepdir so a lot of intermediate files are available if you need me to check something.