Crest version 3 using QCG, issues with water as solvent

[andrewtarzia]

Running QCD on a neutral molecule with water and GFN-FF and ALPB:water seems to fail with the final lines of the output looking like:

  Size       E        De       Detot      Density     Eatom    av. R  Rlast   Volume    Opt
            [Eh]    [kcal]     [kcal]     [u/Å^3]     [kcal]   [bohr] [bohr] [bohr^3]
   Wall Potential too small, increasing size by 5 %
   New scaling factor 0.84
   Wall Potential too small, increasing size by 5 %
   New scaling factor 0.88
   Wall Potential too small, increasing size by 5 %
   New scaling factor 0.93
   Wall Potential too small, increasing size by 5 %
   New scaling factor 0.97
   Wall Potential too small, increasing size by 5 %
   New scaling factor 1.00
error while reading input coordinates

Looking into qcd_temp, it seemed like it could not write the topology for water. Although running the example (https://crest-lab.github.io/crest-docs/page/examples/qcg/example_1.html) with the same crest and xtb binaries caused no issues (it is completing steps at least). However, the energies I obtained from the preoptimisation differ from those in the output file of the example:

  =========================================
  |            Preoptimization            |
  =========================================
  Total Energy of solute:     -15.6531328 Eh
  Total energy of solvent:     -0.3276561 Eh

Crest version:
Version 3.0, Sat Apr 6 18:06:37 UTC 2024 commit (d321183) compiled by 'runner@fv-az778-216'

xtb version:
* xtb version 6.6.1 (8d0f1dd) compiled by 'conda@1efc2f54142f' on 2023-08-01

Any guesses as to what the issue is?

[andrewtarzia]

An update on this, I have rerun the QCG example 1 with this command using the continuous version of crest (Version 3.0.1, Fri May 3 19:27:58 UTC 2024 commit (7c6807c) compiled by 'runner@fv-az976-986'):

crest bac.xyz --qcg water.xyz --nsolv 100 --gfnff --T 12 --alpb water --nofix --xnam xtb

And the output is uploaded at the bottom.

When I do the same run but on my test, larger, neutral system - it gives the same result. It is possible that a newer version of xtb may fix this issue? (I am using xtb version 6.6.1 (8d0f1dd) compiled by 'conda@1efc2f54142f' on 2023-08-01)

       Version 3.0.1, Fri May  3 19:27:58 UTC 2024
       commit (7c6807c) compiled by 'runner@fv-az976-986'
 
   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.
   • P.Pracht, S.Grimme, C.Bannwarth, F.Bohle, S.Ehlert,
     G.Feldmann, J.Gorges, M.Müller, T.Neudecker, C.Plett,
     S.Spicher, P.Steinbach, P.Wesołowski, F.Zeller,
     J. Chem. Phys., 2024, 160, 114110.

   for works involving QCG cite

   • S.Spicher, C.Plett, P.Pracht, A.Hansen, S.Grimme,
     JCTC, 2022, 18 (5), 3174-3189.
   • C.Plett, S. Grimme,
     Angew. Chem. Int. Ed. 2023, 62, e202214477.

   for works involving MECP screening cite

   • P.Pracht, C.Bannwarth, JCTC, 2022, 18 (10), 6370-6385.
 
   Original code
     P.Pracht, S.Grimme, Universität Bonn, MCTC
   with help from (alphabetical order):
     C.Bannwarth, F.Bohle, S.Ehlert, G.Feldmann, J.Gorges,
     S.Grimme, C.Plett, P.Pracht, S.Spicher, P.Steinbach,
     P.Wesolowski, F.Zeller
 
   Online documentation is available at
   https://crest-lab.github.io/crest-docs/
 
   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. See the
   GNU Lesser General Public License (LGPL) for more details.

 Command line input:
 $ /home/atarzia/workingspace/squish-transform/env/bin/crest bac.xyz --qcg water.xyz --nsolv 100 --gfnff --T 12 --alpb water --nofix --xnam /home/atarzia/workingspace/squish-transform/env/bin/xtb

  --gfnff : Use of GFN-FF requested.
  --T 12 (CPUs/Threads selected)
  --alpb water : implicit solvation
  -xnam :
     xtb executable was set to: "/home/atarzia/workingspace/squish-transform/env/bin/xtb"

> Setting up backup calculator ... done.
 ----------------
 Calculation info
 ----------------
> User-defined calculation level:
 : GFN-FF calculation via GFNFF lib
 :   Molecular charge    : 0
 :   Solvation model     : alpb
 :   Solvent             : water
 :   Read dipoles?       : yes
 :   Weight              : 1.00000
 
 
  ========================================
  |           ----------------           |
  |                 Q C G                |
  |           ----------------           |
  |        Quantum Cluster Growth        |
  |       University of Bonn, MCTC       |
  ========================================
   S. Grimme, S. Spicher, C. Plett.
 
   Cite work conducted with this code as

   S. Spicher, C. Plett, P. Pracht, A. Hansen, S. Grimme, JCTC, 2022, 18, 3174-3189.
 
 
   The use of the aISS algorithm is requested (recommend).
   This requires xtb version 6.6.0 or newer.
   xTB-IFF can still be used with the --xtbiff flag.
 
 
  =========================================
  |   quantum cluster growth: INPUT       |
  =========================================
 
  QCG: Only Cluster Generation
 
  input parameters     
  solute                 : bac.xyz
  charge                 : 0
  uhf                    : 0
  solvent                : water.xyz
  # of solvents to add   : 100
  Cluster generated that are above 10 % populated 
  # of CPUs used         : 12
  Solvation model        : water               
  xtb opt level          : normal
  System temperature [K] : 298.1
  RRHO scaling factor    : 0.75
 
 
 Solute geometry
  molecular radius (Bohr**1):   11.20
  molecular area   (Bohr**2): 2554.19
  molecular volume (Bohr**3): 5887.65
 Solvent geometry
  molecular radius (Bohr**1):    3.88
  molecular area   (Bohr**2):  194.90
  molecular volume (Bohr**3):  244.27
 
  radius of solute    :    18.06
  radius of solvent   :     6.25
 
  =========================================
  |            Preoptimization            |
  =========================================
  Total Energy of solute:     -15.6531328 Eh
  Total energy of solvent:     -0.3276561 Eh
 
  ________________________________________________________________________
 
  __________________     Solute Cluster Generation   _____________________
 
  ________________________________________________________________________
 
 
  =========================================
  |   quantum cluster growth: GROW        |
  =========================================
 
 
  Water as solvent recognized, adjusting scaling factor for outer wall pot to 0.80
 
 Solute:
       unit ellipsoid axis a,b,c     :   0.428   0.289   0.283
 Solvent:
       unit ellipsoid axis a,b,c     :   0.384   0.323   0.292
 
  solvent anisotropy            :     1.130
  solute anisotropy             :     1.198
  roff inner wall               :     6.997
  solute max dist               :    39.988
  solvent max dist              :     7.298
  inner unit axis               :     0.529     0.240     0.231
  inner ellipsoid/Bohr          :    36.943    16.790    16.158
  scaling factor outer ellipsoid:     0.800
  outer ellipsoid/Bohr          :    22.848    15.403    15.110
 
  Size       E        De       Detot      Density     Eatom    av. R  Rlast   Volume    Opt
            [Eh]    [kcal]     [kcal]     [u/Å^3]     [kcal]   [bohr] [bohr] [bohr^3]
   Wall Potential too small, increasing size by 5 %
   New scaling factor 0.84
    1   -15.995210   -9.05      -9.05       1.127     -1.659      0.0   0.0    5952.5   normal
    2   -16.338408   -9.75     -18.80       1.129     -1.668      8.7   8.4    6121.7   normal
    3   -16.679637   -8.52     -27.32       1.126     -1.676      8.3  12.7    6316.7   normal
    4   -17.020016   -7.98     -35.30       1.123     -1.685      9.5   7.5    6514.5   normal
    5   -17.365124  -10.95     -46.25       1.123     -1.695      9.6   7.8    6695.3   normal
   Wall Potential too small, increasing size by 5 %
   New scaling factor 0.88
   Wall Potential too small, increasing size by 5 %
   New scaling factor 0.93
   Wall Potential too small, increasing size by 5 %
   New scaling factor 0.97
    6   -17.719856  -16.99     -63.24       1.128     -1.705      9.0   7.1    6840.3   normal
    7   -18.061862   -9.00     -72.25       1.131     -1.714      8.6  11.9    7005.7   normal
    8   -18.406147  -10.43     -82.68       1.133     -1.724      9.3   7.9    7167.4   normal
    9   -18.757864  -15.10     -97.78       1.138     -1.734      9.3   8.7    7316.6   normal
   Wall Potential too small, increasing size by 5 %
   New scaling factor 1.00
   10   -19.100032   -9.11    -106.89       1.138     -1.744      9.3   9.5    7491.4   normal
   11   -19.445383  -11.10    -117.99       1.139     -1.753      9.3   8.8    7661.1   normal

^ It is still running.

[andrewtarzia]

I am not sure how to read these outputs of xtb.out in solute_properties in qcg_tmp, attached

solute_properties: xtb.txt
solvent_properties: xtb_solv.txt

[pprcht]

Hm, that's odd, so it is system (size?) dependent maybe? I can confirm that hexane+water also runs without issue.
Unfortunately the output doesn't provide much insight.
If the failing run has a qcg_tmp/tmp_grow subdirectory, it might be worth looking into what xtb says in its output.

[pprcht]

Ah, ok, I think I have an I idea what is happening. I took the coordinates for your system from the xtb outputs you linked and tried the calculation myself. The error you observe is reproducible:
It seems like there is an SCC convergence issue in the docking computation, which apparently (deliberately) is calculated with some xTB level of theory, despite GFN-FF being selected via the cmd. This means it is system dependent for your system.

As for how to resolve the issue, I need to check what's viable. Sorry for the delay.

[andrewtarzia]

I guess I am not surprised considering the "complexity" of the system - however, they behave nicely with GFN2 (including QCG calculations with GFN2 on the exact same system, which is curious).

No urgency and thank you very much for looking into this! For now, I am using GFN2 and everything seems to be behaving. I wanted to benchmark GFN-FF to see if it was accurate enough because of the huge saving. Especially because the QCG calculations are simply to setup reasonable structures for future calculations.