Inconsistent energy and force when atoms are near the periodic boundries

[zyt0y]

Describe the bug
In principle, the total energy and force should have translational invariance under the periodic boundtry condition. When I translate a periodic system (water on graphene), the total energy and maximal atomic force are not always the same.

To Reproduce
Setup:
python: 3.6.3
ase: 3.17.0
xtb-python: 20.1
xtb: 6.3.0

The structure for water trimer on graphene is built first. Then I translate the system using a equally spaced distance.

for i in range(20):
    atoms.translate([0, -i*0.1, 0])

Energy and force are calculated using following command:

atoms.set_calculator(XTB(method="GFN0-xTB"))
#atoms.set_calculator(XTB(method="GFN1-xTB"))
atoms.get_forces()

Energy and maximal atomic force are plotted for structrues at different displacements and with different methods. Reference calculations are performed using DFTB+. The carbon atoms near the cell boundaries are plotted in blue colors for better visual comparisons.

Trajectory files in ASE format pbc.zip are also provided here.

Expected behaviour
GFN0-xTB and GFN1-xTB should havetranslational invariance under the periodic boundtry condition and show similar curve as periodic DFTB+.

Additional context
I am running xtb using xtb-python interface. The problem should lie on xtb part, I guess.

[awvwgk]

Could be an issue with the Wigner–Seitz cell (WSC) used for the cyclic cluster model (CCM) which is generating the periodic boundary conditions for the Hamiltonian. The current WSC generator can miss certain crucial interacting pairs if they are not included in a 3×3×3 supercell, i.e. if they are just in the 2, 0, 0 cell.

First, we should check if using the full real-space cutoff rather than the CCM is showing the same behaviour.

[zyt0y]

get_realspace_cutoff is defined in https://github.com/grimme-lab/xtb/blob/master/src/pbc.f90#L31 . Should I look into all places where get_realspace_cutoff is called?

[awvwgk]

Sorry, I'm currently short of time and haven't investigated further yet. It would certainly help to nail down the energy expression that is not translational invariant. The Hamiltonian was a wild guess, after checking the GFN1-xTB Hamiltonian is done in real-space, while the GFN0-xTB Hamiltonian is done with the CCM, so this might not be the issue.

The get_realspace_cutoff should only have very limited use in the code, I moved most of the boundary condition handling to https://github.com/grimme-lab/xtb/blob/master/src/type/latticepoint.f90. Since you don't deform the cell, the lattice point generator should always generate the same amount of lattice points for all cases.

[zyt0y]

I will try to figure out how PBC is working in xtb and see if I can indenitify where the problem is.

[zulissi]

Hi - we've also been noticing inconsistencies with PBC for molecules on surfaces in GFN0. Any luck tracking down the issue @tangzeyuan ?

[zyt0y]

Hi - we've also been noticing inconsistencies with PBC for molecules on surfaces in GFN0. Any luck tracking down the issue @tangzeyuan ?

No clues yet.

[awvwgk]

I haven't had time to investigate yet, either. Maybe I can clear enough projects by the end of this month to find time to work on the PBC implementation again. But can't promise.

[manuelmelle]

Same problem found:

I am running some preliminary calculations on a layered system, CoPS3, and found also some strange behaviour with the FIRE optimization of the periodic system with GFN0.

Atoms near the PBC boundaries show larger displacements which also could indicate that the boundary atoms do not interact properly with their neighbours through the PBC.

In this context, I would like to know if **it is possible to run fixed PBC optimizations or to select which cell parameters are optimized.

ps: great method and software BTW!

[zulissi]

FWIW, I found that using XTB for single-points in ASE and doing the optimization there (controlling cell parameters and free atoms etc) worked better for me than the internal optimizers. It still had the issues with PBC effects though.

[manuelmelle]

Thanks for the advice, I´ll try ASE and try to have a look at the source if I have some time.

From what you say, there seems to be a clear problem with the forces as mentioned, that depends on the system.
I computed other systems without trouble, but the one I am trying now is really problematic.

[zulissi]

Inspired by the suggestion that it might be a dispersion problem from grimme-lab/xtb-python#34, I removed the dispersion calculation and recompiled and it appears to be consistent now (at least on the same test structures he posted there).

FWIW, it would also be nice to have a settings to turn off the dispersion. We're using this as a base calculator to fit to DFT data that doesn't have dispersion, so regardless of PBC issues it's helpful to have it off.

[awvwgk]

Thanks for reminding me on this, we recently were able to track this issue down to the generation of the lattice points (many thanks to @Thomas3R for spotting this).

I'll work on a patch for the 6.4.0 release this week.

[awvwgk]

Sorry for not following up on this issue for some time. I decided for a change in strategy for the periodic boundary conditions in xtb, because trying to implement everything here from scratch turned out quite time consuming and error prone.

Instead, I'm now working on a unified library implementation, which will be available in dftb+, xtb, qcxms, ... as common backend for the xTB Hamiltonians. I think this is the better way to go instead of creating monolithic projects. The next step on the periodic boundary conditions is therefore on getting a full k-point periodic xTB Hamiltonian into the dftb+ project.