Rebo2Scr seems to calculate randomly wrong forces when using OpenMP Parallel

[threich]

When calculating the forces in an amorphous C slab using Rebo2Scr within the python ase-interface and several threads, the result seems to be random. When using the following script (the script and the structure are attached in the zip file below)

from atomistica import Rebo2Scr
import ase.io

a = ase.io.read('struct.extxyz')
calc = Rebo2Scr()
a.set_calculator(calc)

for i in range(30):
    print(a.get_forces() ** 2).sum(1).max()

the result obtained using OMP_NUM_THREADS=1 is in all print statements 9.83344680117e-05.
However, when using e.g. OMP_NUM_THREADS=20 the results are very different, e.g.:
0.0787229654406
9.83344680117e-05
9.83344680117e-05
0.204342136593
0.00617123803722
9.83344680117e-05
0.293423808598
0.00141455572593
9.83344680117e-05
9.83344680117e-05
0.0624515880291
0.0556148664835
9.83344680117e-05
0.00617123803722
9.83344680117e-05
9.83344680117e-05
9.83344680117e-05
9.83344680117e-05
0.0952089154343
0.10312023489
9.83344680117e-05
0.134743528249
0.0787578185209
9.83344680117e-05
1.17992497141
0.634367892818
9.83344680117e-05
9.83344680117e-05
0.0156597251095
0.0891190587402

The scattering seems to increase with the number of atoms and the number of threads.
The tested versions include v0.4.1, 0.6.0 and 0.7.0 using ase 3.16.2.
MWE.zip

The problem vanishes when using Rebo2 instead of Rebo2Scr.

[threich]

The problem can also be observed using a small cluster system, but less frequently. The script

from atomistica import Rebo2Scr
import numpy as np
from ase import Atoms

positions = [[3.59005750, 4.30036973, 9.46763746],
             [5.54900821, 5.54722436, 8.41101491],
             [5.98411364, 3.34799370, 9.65966376],
             [5.77213711, 4.93952228, 11.23127506],
             [2.97728181, 5.48820437, 10.29817917],
             [4.54303387, 3.15626810, 10.00535256],
             [6.19751149, 4.90048844, 9.65999527],
             [4.22467444, 4.63103841, 8.04930501]]

a = Atoms('C8', positions=positions)
a.center(vacuum=4)
a.pbc = False
calc = Rebo2Scr()
a.set_calculator(calc)

forces_ref = a.get_forces()
fmax_ref = np.sqrt((forces_ref**2).sum(axis=1).max())

for i in range(1000000):
    forces = a.get_forces()
    fmax = np.sqrt((forces**2).sum(axis=1).max())
    if abs(fmax-fmax_ref) > 0.01:
        print(i, fmax, fmax_ref)
        print(forces, forces_ref)
        break

delivers e.g.

(293, 5.5817569662132813, 5.5985440310988368)
(array([[-2.68364478, -0.8806511 , 1.0967975 ],
[-2.35801924, -4.38387914, 1.34969855],
[-3.42333517, 2.14237654, 1.57051134],
[ 2.45025748, 1.02652469, -3.2583111 ],
[ 0.58812198, -3.92490632, -3.03422628],
[ 1.25960847, 3.88804171, -2.98347877],
[ 2.14322211, -1.5944522 , 1.62993681],
[ 2.02378915, 3.72694582, 3.62907194]]),
array([[-2.59503897, -0.92833751, 1.10798925],
[-2.37908816, -4.40293054, 1.36957525],
[-3.43565141, 2.15032581, 1.5694922 ],
[ 2.45025748, 1.02652469, -3.2583111 ],
[ 0.58812198, -3.92490632, -3.03422628],
[ 1.33981683, 4.08392998, -3.03838738],
[ 1.99426236, -1.73263258, 1.63765885],
[ 2.03731989, 3.72802647, 3.6462092 ]]))

[pastewka]

Which compiler was that? Did you try both gcc and Intel?

[threich]

The reported tests were done using atomistica compiled with Intel compilers 18.0.2 20180210.
Just now I tried a GCC 7.3.0. compiled version, but the issue persists.

[pastewka]

Is this resolved with the latest commits that switch of the splines?

[threich]

No. Initially, the problem disappeared with the commits from 6.2., but the fix of the cutoffs (12.2) reintroduced the issue.

[pastewka]

@threich I just fixed a bug in the neighor list construction. Can you check whether this problem still occurs with the current master? It may be related to #30

[threich]

Unfortunately, I can still reproduce the error with the current master using the MWE.zip-example and e.g. OMP_NUM_THREADS=20.

[pastewka]

Okay, thanks for checking.