Coverage for gpaw/test/pw/test_stresstest.py: 100%

1import numpy as np 

2 

3from gpaw.grid_descriptor import GridDescriptor 

4from gpaw.spline import Spline 

5import gpaw.mpi as mpi 

6from gpaw.pw.descriptor import PWDescriptor 

7from gpaw.pw.lfc import PWLFC 

8 

9 

10def test_stress(): 

11 x = 2.0 

12 rc = 3.5 

13 r = np.linspace(0, rc, 100) 

14 

15 n = 40 

16 a = 8.0 

17 cell_cv = np.array([[a, 0.5, -1], [0, a, 2], [-1, 0, a + 1]]) 

18 gd = GridDescriptor((n, n, n), cell_cv, comm=mpi.serial_comm) 

19 

20 a_R = gd.empty() 

21 z = np.linspace(0, n, n, endpoint=False) 

22 a_R[:] = 2 + np.sin(2 * np.pi * z / n) 

23 

24 spos_ac = np.array([(0.15, 0.45, 0.95)]) 

25 

26 pd = PWDescriptor(45, gd) 

27 a_G = pd.fft(a_R) 

28 

29 s = Spline.from_data(0, rc, 2 * x**1.5 / np.pi * np.exp(-x * r**2)) 

30 p = Spline.from_data(1, rc, 2 * x**1.5 / np.pi * np.exp(-x * r**2)) 

31 d = Spline.from_data(2, rc, 2 * x**1.5 / np.pi * np.exp(-x * r**2)) 

32 

33 lfc = PWLFC([[s, p, d]], pd) 

34 lfc.set_positions(spos_ac) 

35 b_LG = pd.zeros(9) 

36 lfc.add(b_LG, {0: np.eye(9)}) 

37 e1 = pd.integrate(a_G, b_LG) 

38 assert abs(lfc.integrate(a_G)[0] - e1).max() < 1e-11 

39 

40 s1 = [] 

41 for i in range(9): 

42 x = [0, 0, 0, 0, 0, 0, 0, 0, 0] 

43 x[i] = 1 

44 s1.append(lfc.stress_tensor_contribution(a_G, {0: x}) - 

45 np.eye(3) * e1[i]) 

46 

47 x = 1e-6 

48 for dist in [[[x, 0, 0], [0, 0, 0], [0, 0, 0]], 

49 [[0, 0, 0], [0, x, 0], [0, 0, 0]], 

50 [[0, 0, 0], [0, 0, 0], [0, 0, x]], 

51 [[0, x, 0], [x, 0, 0], [0, 0, 0]], 

52 [[0, 0, x], [0, 0, 0], [x, 0, 0]], 

53 [[0, 0, 0], [0, 0, x], [0, x, 0]]]: 

54 e = dist + np.eye(3) 

55 c_cv = np.dot(cell_cv, e) 

56 gd = GridDescriptor((n, n, n), c_cv, comm=mpi.serial_comm) 

57 pd = PWDescriptor(45, gd) 

58 aa_G = a_G / np.linalg.det(e) 

59 lfc = PWLFC([[s, p, d]], pd) 

60 lfc.set_positions(spos_ac) 

61 b_LG = pd.zeros(9) 

62 lfc.add(b_LG, {0: np.eye(9)}) 

63 e2 = pd.integrate(aa_G, b_LG) 

64 s2 = (e2 - e1) / x 

65 error = (np.array(s1) * dist).sum(1).sum(1) / x - s2 

66 print(s2, abs(error).max()) 

67 assert abs(error).max() < 2e-6