Inconsistent results for the identical system with different order of coordinates

[FanGuozheng]

Describe the bug
For an identical system, when only switching orders of coordinates, the DFTB+ results with solvation effect, including charges, energies are different. This only exists when using GeneralizedBorn. If we turn off GeneralizedBorn or switch to COSMO or Vacuum, the DFTB+ results are the same, and which should be the same.

To Reproduce

DFTB input:
Geometry = genFormat {
<<< geo.gen
}
Hamiltonian = DFTB {
SCC = Yes
SCCTolerance = 1.0E-5
ShellResolvedSCC = No
MaxSCCIterations = 500
Mixer = Broyden {
MixingParameter = 0.2
}
Filling = Fermi {
Temperature [K] = 310.15
}
MaxAngularMomentum = {
H = "s"
C = "p"
N = "p"
O = "p"
}
HubbardDerivs = {
N = -0.1535
C = -0.1492
O = -0.1575
H = -0.1857
}
ThirdOrderFull = Yes
HCorrection = Damping {
Exponent = 4.00
}
Solvation = GeneralizedBorn { # GFN2-xTB/GBSA(water)
Solvent = fromConstants {
Epsilon = 80.2
MolecularMass [amu] = 18.0
Density [kg/l] = 1.0
}
Temperature = 0.0009821802
FreeEnergyShift [kcal/mol] = 1.16556316
BornScale = 1.55243817
BornOffset = 2.462811043694508E-02
Radii = vanDerWaalsRadiiD3 {}
Descreening = Values {
H = 0.71893869
C = 0.74298311
N = 0.90261230
O = 0.75369019
}
SASA {
ProbeRadius = 1.843075777670416
Radii = vanDerWaalsRadiiD3 {}
SurfaceTension = Values {
H = -3.34983060E-01
C = -7.47690650E-01
N = -2.31291292E+00
O = 9.17979110E-01
}
}
HBondCorr = Yes
HBondStrength = Values {
H = -7.172800544988973E-02
C = -2.548469535762511E-03
N = -1.976849501504001E-02
O = -8.462476828587280E-03
}
}
SlaterKosterFiles = Type2FileNames {
Prefix = "../3ob-3-1/"
Separator = "-"
Suffix = ".skf"
}
Solver = MAGMA {}
Charge = -1
}

Options = {
WriteChargesAsText = Yes
}

First geometry input:
38 C
C O N H
1 1 2.9822000000E+01 -5.2810000000E+00 -1.2905000000E+01
2 1 3.0345000000E+01 -5.0570000000E+00 -1.1492000000E+01
3 1 2.9876000000E+01 -3.7380000000E+00 -1.0890000000E+01
4 1 3.0048000000E+01 -3.7880000000E+00 -9.3770000000E+00
5 1 3.0641000000E+01 -2.5090000000E+00 -8.7970000000E+00
6 1 3.0065000000E+01 -2.2350000000E+00 -7.4200000000E+00
7 2 2.9035000000E+01 -2.7470000000E+00 -6.9990000000E+00
8 3 3.0751000000E+01 -1.3950000000E+00 -6.6690000000E+00
9 2 3.0551000000E+01 -1.3740000000E+00 -5.4510000000E+00
10 1 2.9305000000E+01 -3.9800000000E+00 -1.3503000000E+01
11 1 2.7835000000E+01 -4.1240000000E+00 -1.3816000000E+01
12 2 2.6977000000E+01 -3.7640000000E+00 -1.3019000000E+01
13 3 2.7571000000E+01 -4.6270000000E+00 -1.5027000000E+01
14 1 2.6416000000E+01 -4.5080000000E+00 -1.5727000000E+01
15 1 2.5566000000E+01 -5.6240000000E+00 -1.5775000000E+01
16 1 2.4365000000E+01 -5.5660000000E+00 -1.6482000000E+01
17 1 2.4019000000E+01 -4.3890000000E+00 -1.7147000000E+01
18 1 2.4874000000E+01 -3.2750000000E+00 -1.7109000000E+01
19 1 2.6077000000E+01 -3.3310000000E+00 -1.6405000000E+01
20 4 2.9082300000E+01 -6.0869000000E+00 -1.2886300000E+01
21 4 3.0634900000E+01 -5.6745000000E+00 -1.3518500000E+01
22 4 3.1436200000E+01 -5.0833000000E+00 -1.1489600000E+01
23 4 3.0026300000E+01 -5.9002000000E+00 -1.0875500000E+01
24 4 2.8822600000E+01 -3.5199000000E+00 -1.1073600000E+01
25 4 3.0484300000E+01 -2.9221000000E+00 -1.1284900000E+01
26 4 3.0710700000E+01 -4.6018000000E+00 -9.0740000000E+00
27 4 2.9074900000E+01 -4.0708000000E+00 -8.9702000000E+00
28 4 3.0367400000E+01 -1.6483000000E+00 -9.4096000000E+00
29 4 3.1731800000E+01 -2.5520000000E+00 -8.7841000000E+00
30 4 3.1578100000E+01 -9.4740000000E-01 -7.0362000000E+00
31 4 2.9825100000E+01 -3.8755000000E+00 -1.4458200000E+01
32 4 2.9515000000E+01 -2.9982000000E+00 -1.3092600000E+01
33 4 2.8350400000E+01 -4.9973000000E+00 -1.5549900000E+01
34 4 2.5830000000E+01 -6.5356000000E+00 -1.5257600000E+01
35 4 2.3708400000E+01 -6.4231000000E+00 -1.6512600000E+01
36 4 2.3089500000E+01 -4.3345000000E+00 -1.7694400000E+01
37 4 2.4598500000E+01 -2.3692000000E+00 -1.7628700000E+01
38 4 2.6730000000E+01 -2.4705000000E+00 -1.6378600000E+01

second geometry input (move first 10 atoms in above geometry to the end in this geometry):
38 C
C O N H
1 1 2.7835000000E+01 -4.1240000000E+00 -1.3816000000E+01
2 2 2.6977000000E+01 -3.7640000000E+00 -1.3019000000E+01
3 3 2.7571000000E+01 -4.6270000000E+00 -1.5027000000E+01
4 1 2.6416000000E+01 -4.5080000000E+00 -1.5727000000E+01
5 1 2.5566000000E+01 -5.6240000000E+00 -1.5775000000E+01
6 1 2.4365000000E+01 -5.5660000000E+00 -1.6482000000E+01
7 1 2.4019000000E+01 -4.3890000000E+00 -1.7147000000E+01
8 1 2.4874000000E+01 -3.2750000000E+00 -1.7109000000E+01
9 1 2.6077000000E+01 -3.3310000000E+00 -1.6405000000E+01
10 4 2.9082300000E+01 -6.0869000000E+00 -1.2886300000E+01
11 4 3.0634900000E+01 -5.6745000000E+00 -1.3518500000E+01
12 4 3.1436200000E+01 -5.0833000000E+00 -1.1489600000E+01
13 4 3.0026300000E+01 -5.9002000000E+00 -1.0875500000E+01
14 4 2.8822600000E+01 -3.5199000000E+00 -1.1073600000E+01
15 4 3.0484300000E+01 -2.9221000000E+00 -1.1284900000E+01
16 4 3.0710700000E+01 -4.6018000000E+00 -9.0740000000E+00
17 4 2.9074900000E+01 -4.0708000000E+00 -8.9702000000E+00
18 4 3.0367400000E+01 -1.6483000000E+00 -9.4096000000E+00
19 4 3.1731800000E+01 -2.5520000000E+00 -8.7841000000E+00
20 4 3.1578100000E+01 -9.4740000000E-01 -7.0362000000E+00
21 4 2.9825100000E+01 -3.8755000000E+00 -1.4458200000E+01
22 4 2.9515000000E+01 -2.9982000000E+00 -1.3092600000E+01
23 4 2.8350400000E+01 -4.9973000000E+00 -1.5549900000E+01
24 4 2.5830000000E+01 -6.5356000000E+00 -1.5257600000E+01
25 4 2.3708400000E+01 -6.4231000000E+00 -1.6512600000E+01
26 4 2.3089500000E+01 -4.3345000000E+00 -1.7694400000E+01
27 4 2.4598500000E+01 -2.3692000000E+00 -1.7628700000E+01
28 4 2.6730000000E+01 -2.4705000000E+00 -1.6378600000E+01
29 1 2.9822000000E+01 -5.2810000000E+00 -1.2905000000E+01
30 1 3.0345000000E+01 -5.0570000000E+00 -1.1492000000E+01
31 1 2.9876000000E+01 -3.7380000000E+00 -1.0890000000E+01
32 1 3.0048000000E+01 -3.7880000000E+00 -9.3770000000E+00
33 1 3.0641000000E+01 -2.5090000000E+00 -8.7970000000E+00
34 1 3.0065000000E+01 -2.2350000000E+00 -7.4200000000E+00
35 2 2.9035000000E+01 -2.7470000000E+00 -6.9990000000E+00
36 3 3.0751000000E+01 -1.3950000000E+00 -6.6690000000E+00
37 2 3.0551000000E+01 -1.3740000000E+00 -5.4510000000E+00
38 1 2.9305000000E+01 -3.9800000000E+00 -1.3503000000E+01

Error output:
The Mulliken charge for the same atom will have 1E-1 e level difference, to minimize the text here, the following part show the total energy difference.

First geometry:
Total energy: -45.8655611501 H -1248.0654 eV
Second geometry:
Total energy: -45.8708466853 H -1248.2092 eV

Expected behaviour
Get the same results for the identical system with different order of coordinates when using GeneralizedBorn.

Additional context

[tsihyoung]

The problem lies in the incorrect nNeighbour, iNeighbour and neighDist2 returned by subroutine updateNeighbourList.

By default, symm = .false., thus lpIAtom2: do iAtom2 = 1, iAtom2End only loops to iAtom1. This will underestimate the number of neighbours for atoms with index larger than 1.

To fix it, I suggest changing symm to .true..

[tsihyoung]

Also, the problem is not related to GB model itself, but is rooted in SASA.

Another problem, although not related to the issue, is that the OpenMP directive is not correct in subroutine getSASA. It should be !$omp shared(sasa, dsdr) instead of !$omp reduction(+: sasa, dsdr) as these two arrays are not undergoing summation inside the loop.

[aradi]

@FanGuozheng thanks, I can reproduce it, I'll try to fix it. @tsihyoung Thanks for the hints, they are very valuable. As almost all parts of DFTB+ work with the asymmetric neighbor list, I'll probably try to fix the SASA instead.

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs.

[stale]

This stale issue has been automatically closed.

[bhourahine]

@aradi is this still a live issue?

[aradi]

Sure, it is still an issue.