A difficult case of xTB hessian

[xiki-tempula]

I'm testing the xTB routine and for most molecules (>95%). The xTB hessian behaves fine (give similar energy to PDB0-D3(BJ)/def2-TZVP). However there is one molecule that gives vastly different energy.

The way that I test the force field is that I generate 100 conformations of the molecule, and compare the MM energy of these 100 conformations.

The error seems to come from the angle part. However, I cannot diagnose the issue, I wonder if you mind help me on this? Thanks. The molecule is

HETATM    1  C0  lig     1      -1.960  21.550 -27.330  0.00  0.00           C  
HETATM    2  C1  lig     1      -1.290  22.410 -28.200  0.00  0.00           C  
HETATM    3  C2  lig     1      -1.960  22.950 -29.290  0.00  0.00           C  
HETATM    4  C3  lig     1      -3.300  22.620 -29.540  0.00  0.00           C  
HETATM    5  C4  lig     1      -3.970  21.760 -28.650  0.00  0.00           C  
HETATM    6  C5  lig     1      -3.290  21.230 -27.560  0.00  0.00           C  
HETATM    7 Cl6  lig     1      -5.640  21.350 -28.920  0.00  0.00          CL  
HETATM    8  C7  lig     1      -4.040  23.150 -30.710  0.00  0.00           C  
HETATM    9  O8  lig     1      -4.230  22.420 -31.670  0.00  0.00           O  
HETATM   10  N9  lig     1      -4.510  24.410 -30.650  0.00  0.00           N  
HETATM   11  C10 lig     1      -5.290  25.110 -31.580  0.00  0.00           C  
HETATM   12  C11 lig     1      -5.960  24.550 -32.680  0.00  0.00           C  
HETATM   13  C12 lig     1      -6.780  25.360 -33.460  0.00  0.00           C  
HETATM   14  N13 lig     1      -6.920  26.650 -33.180  0.00  0.00           N  
HETATM   15  C14 lig     1      -6.300  27.230 -32.160  0.00  0.00           C  
HETATM   16  C15 lig     1      -5.510  26.460 -31.310  0.00  0.00           C  
HETATM   17  N16 lig     1      -6.540  28.580 -31.880  0.00  0.00           N  
HETATM   18  C17 lig     1      -5.710  29.450 -31.230  0.00  0.00           C  
HETATM   19  O18 lig     1      -4.620  29.120 -30.810  0.00  0.00           O  
HETATM   20  C19 lig     1      -6.230  30.850 -31.030  0.00  0.00           C  
HETATM   21 Cl20 lig     1      -1.110  24.050 -30.330  0.00  0.00          CL  
HETATM   22  H21 lig     1      -6.760  31.250 -31.890  0.00  0.00           H  
HETATM   23  H22 lig     1      -1.440  21.130 -26.480  0.00  0.00           H  
HETATM   24  H23 lig     1      -0.260  22.660 -28.020  0.00  0.00           H  
HETATM   25  H24 lig     1      -3.810  20.560 -26.890  0.00  0.00           H  
HETATM   26  H25 lig     1      -4.320  24.980 -29.840  0.00  0.00           H  
HETATM   27  H26 lig     1      -5.840  23.500 -32.910  0.00  0.00           H  
HETATM   28  H27 lig     1      -7.310  24.940 -34.300  0.00  0.00           H  
HETATM   29  H28 lig     1      -5.060  26.910 -30.430  0.00  0.00           H  
HETATM   30  H29 lig     1      -7.400  29.010 -32.180  0.00  0.00           H  
HETATM   31  H30 lig     1      -6.890  30.870 -30.160  0.00  0.00           H  
HETATM   32  H31 lig     1      -5.460  32.770 -30.540  0.00  0.00           H  
HETATM   33  C32 lig     1      -5.140  31.760 -30.690  0.00  0.00           C  
HETATM   34  C33 lig     1      -4.010  31.220 -29.860  0.00  0.00           C  
HETATM   35  C34 lig     1      -3.800  31.510 -31.320  0.00  0.00           C  
HETATM   36  H35 lig     1      -3.550  31.930 -29.140  0.00  0.00           H  
HETATM   37  H36 lig     1      -4.090  30.140 -29.580  0.00  0.00           H  
HETATM   38  H37 lig     1      -3.190  32.410 -31.550  0.00  0.00           H  
HETATM   39  H38 lig     1      -3.730  30.620 -31.990  0.00  0.00           H  

[selimsami]

Could you please send all the relevant files?

[xiki-tempula]

So the file for reproducing the xTB topology is in folder xTB.

I don't have any QM software yet so I did a hacky way where I generated the QM hessian with psi4 (QM/lig_hessian.npy) and use the xTB routine to generate the new topology, I think you could regenerate the QM topology with any other QM engine. The topology file that I have generate is in QM/gas.top.

100 conformers are in MM/lig.gro and MM/traj.xtc, I used run.sh to generate the MM energy profile and plot.py to analysis the result. energy.pdf As is shown here that the QM one is similar to other force field and seems ok but the xTB one is very wrong. Archive.zip Thanks.

[selimsami]

; flexible dihedrals 15 17 18 20 3 -14476.443 2997.159 1898.391 -752.774 -1962.776 -2373.229

Something seem to have gone wrong with one of the dihedral fitting. You should probably have received a warning about this.

[xiki-tempula]

@selimsami Yes (scan_data_CN_H12C10N2O1_75efeccb473e9b543b61dc50f5d1fc16), it seems a bit weird that the same dihedral in the on the other one is fine where the only thing that have changed is the hessian.

[xiki-tempula]

@selimsami Wait let me do some further investigation on this.

[selimsami]

Starting geometry of the MM optimization is based on the QM optimized fragment. The difference between those is where I would start the investigation

[xiki-tempula]

Somehow I cannot reproduce the issue now.