Coverage for gpaw/test/test_d3_extension.py: 10%

1import pytest 

2from gpaw.mpi import world 

3from gpaw.new.ase_interface import GPAW 

4from gpaw import restart 

5from gpaw.new.extensions import D3 

6import numpy as np 

7from ase import Atoms 

8 

9 

10@pytest.mark.parametrize('parallel', [(1, 1), (1, 2), (2, 2), (2, 1)]) 

11@pytest.mark.parametrize('mode', [{'name': 'pw', 'ecut': 300}, 'lcao']) 

12def test_d3_extensions(mode, parallel, in_tmp_dir, gpaw_new, dftd3): 

13 if not gpaw_new: 

14 pytest.skip('Only GPAW new.') 

15 

16 from ase.calculators.dftd3 import PureDFTD3 

17 domain, band = parallel 

18 if world.size < domain * band: 

19 pytest.skip('Not enough cores for this test.') 

20 if world.size > domain * band * 2: 

21 pytest.skip('Too many cores for this test.') 

22 

23 # 1. Create a calculation with a particular list of extensions. 

24 def get_atoms(): 

25 from ase.build import molecule 

26 atoms = molecule('H2') 

27 atoms[0].position[1] += 1 

28 atoms.center(vacuum=2) 

29 atoms.set_pbc((True, True, True)) 

30 return atoms 

31 

32 def D3ref(atoms): 

33 atoms = atoms.copy() 

34 atoms.calc = PureDFTD3(xc='PBE') 

35 return atoms.get_potential_energy(), atoms.get_forces() 

36 

37 atoms = get_atoms() 

38 

39 def get_calc(atoms): 

40 # To test multiple extensions, create two sprigs which add 

41 # up to k=ktot, which is what is tested in this test 

42 calc = GPAW(extensions=[D3(xc='PBE')], 

43 symmetry='off', 

44 parallel={'band': band, 'domain': domain}, 

45 kpts=(2, 1, 1), 

46 mode=mode) 

47 atoms.calc = calc 

48 return calc 

49 

50 calc = get_calc(atoms) 

51 

52 E, F = atoms.get_potential_energy(), atoms.get_forces() 

53 D3_E, D3_F = D3ref(atoms) 

54 

55 # Write the GPW file for the restart test later on (4.) 

56 print('Wrote the potential energy', E) 

57 calc.write('calc.gpw') 

58 

59 # 2. Test that moving the atoms works after an SFC convergence 

60 atoms.positions[0, 2] -= 0.1 

61 movedE, movedF = atoms.get_potential_energy(), atoms.get_forces() 

62 

63 movedD3_E, movedD3_F = D3ref(atoms) 

64 # Reset atoms to their original positions 

65 atoms.positions[0, 2] += 0.1 

66 

67 # 3. Calculate a reference result without extensions 

68 calc = GPAW(mode=mode, 

69 kpts=(2, 1, 1), 

70 symmetry='off') 

71 atoms.calc = calc 

72 

73 E0, F0 = atoms.get_potential_energy(), atoms.get_forces() 

74 

75 # Manually evaluate the spring energy, and compare forces 

76 assert E == pytest.approx(E0 + D3_E) 

77 assert F == pytest.approx(F0 + D3_F) 

78 

79 # Evaluate the reference energy and forces also for the moved atoms 

80 atoms.positions[0, 2] -= 0.1 

81 movedE0, movedF0 = atoms.get_potential_energy(), atoms.get_forces() 

82 assert movedE == pytest.approx(movedE0 + movedD3_E) 

83 assert movedF == pytest.approx(movedF0 + movedD3_F) 

84 

85 # 4. Test restarting from a file 

86 atoms, calc = restart('calc.gpw', Class=GPAW) 

87 # Make sure the cached energies and forces are correct 

88 # without a new calculation 

89 assert E == pytest.approx(atoms.get_potential_energy()) 

90 assert F == pytest.approx(atoms.get_forces()) 

91 

92 if mode == 'lcao': 

93 # See issue #1369 

94 return 

95 

96 # Make sure the recalculated energies are forces are correct 

97 atoms.set_positions(atoms.get_positions() + 1e-10) 

98 assert E == pytest.approx(atoms.get_potential_energy(), abs=1e-5) 

99 assert F == pytest.approx(atoms.get_forces(), abs=1e-5) 

100 

101 # 5. Test full blown relaxation. 

102 from ase.optimize import BFGS 

103 atoms = get_atoms() 

104 calc = get_calc(atoms) 

105 relax = BFGS(atoms) 

106 relax.run() 

107 nsteps = relax.nsteps 

108 assert atoms.get_distance(0, 1) == pytest.approx(0.76915, abs=1e-2) 

109 Egs = atoms.get_potential_energy() 

110 L = atoms.get_distance(0, 1) 

111 

112 # 6. Test restarting from a relaxation. 

113 atoms = get_atoms() 

114 calc = get_calc(atoms) 

115 relax = BFGS(atoms, restart='relax_restart') 

116 for _, _ in zip(relax.irun(), range(3)): 

117 pass 

118 calc.write('restart_relax.gpw') 

119 atoms, calc = restart('restart_relax.gpw', Class=GPAW) 

120 relax = BFGS(atoms, restart='relax_restart') 

121 relax.run() 

122 

123 assert relax.nsteps + 3 == nsteps 

124 assert atoms.get_distance(0, 1) == pytest.approx(L, abs=1e-2) 

125 assert atoms.get_potential_energy() == pytest.approx(Egs, abs=1e-4) 

126 

127 

128@pytest.mark.parametrize('parallel', [(1, 1), (1, 2), (2, 2), (2, 1)]) 

129def test_d3_stress(parallel, in_tmp_dir, dftd3): 

130 from ase.calculators.dftd3 import DFTD3 

131 from ase.optimize import CellAwareBFGS 

132 from ase.build import bulk 

133 from ase.filters import FrechetCellFilter 

134 from gpaw.new.ase_interface import GPAW 

135 

136 domain, band = parallel 

137 if world.size < domain * band: 

138 pytest.skip('Not enough cores for this test.') 

139 if world.size > domain * band * 2: 

140 pytest.skip('Too many cores for this test.') 

141 

142 def get_atoms(): 

143 atoms = bulk('C', a=3.5) 

144 atoms.set_cell(atoms.get_cell(), 

145 scale_atoms=True) 

146 return atoms 

147 

148 kwargs = dict(xc='LDA', 

149 parallel={'band': band, 'domain': domain}, 

150 kpts=(2, 2, 2), txt='relax', 

151 convergence={'density': 1e-5}, 

152 mode=dict(name='pw', ecut=300)) 

153 

154 def get_calc(x): 

155 return GPAW(**kwargs, **x) 

156 

157 # 1. Old fashioned D3 calculation 

158 atoms = get_atoms() 

159 atoms.calc = DFTD3(xc='PBE', dft=get_calc({})) 

160 relax = CellAwareBFGS(FrechetCellFilter(atoms, exp_cell_factor=1), 

161 restart='restart_oldfashioned') 

162 relax.run(smax=0.001) 

163 atoms_old_ref = atoms.copy() 

164 E_ref = atoms.get_potential_energy() 

165 

166 # 2. New style D3 calculation 

167 atoms = get_atoms() 

168 atoms.calc = get_calc(dict(extensions=[D3(xc='PBE')])) 

169 relax = CellAwareBFGS(FrechetCellFilter(atoms, exp_cell_factor=1), 

170 restart='restart_new') 

171 relax.run(smax=0.001) 

172 nsteps = relax.nsteps 

173 

174 assert np.allclose(atoms.cell, atoms_old_ref.cell) 

175 assert np.allclose(atoms.get_scaled_positions(), 

176 atoms_old_ref.get_scaled_positions()) 

177 assert E_ref == pytest.approx(atoms.get_potential_energy(), abs=1e-4) 

178 

179 # 3. Restarting geometry relaxation of new style D3 calculation 

180 atoms = get_atoms() 

181 atoms.calc = get_calc(dict(extensions=[D3(xc='PBE')])) 

182 relax = CellAwareBFGS(FrechetCellFilter(atoms, exp_cell_factor=1), 

183 restart='restart_cont') 

184 relax.smax = 1e-4 

185 for _, _ in zip(relax.irun(), range(3)): 

186 pass 

187 atoms.calc.write('restart_cell_relax.gpw') 

188 atoms, calc = restart('restart_cell_relax.gpw', Class=GPAW) 

189 relax = CellAwareBFGS(FrechetCellFilter(atoms, exp_cell_factor=1), 

190 restart='restart_cont') 

191 relax.run(smax=0.001) 

192 

193 assert relax.nsteps + 3 == nsteps 

194 

195 assert np.allclose(atoms.cell, atoms_old_ref.cell, rtol=1e-5, atol=1e-5) 

196 assert np.allclose(atoms.get_scaled_positions(), 

197 atoms_old_ref.get_scaled_positions()) 

198 assert E_ref == pytest.approx(atoms.get_potential_energy(), abs=1e-4) 

199 

200 

201def test_d3_isolated_atom(dftd3): 

202 atoms = Atoms('He') 

203 atoms.center(vacuum=3) 

204 calc = GPAW(xc='PBE', 

205 extensions=[D3(xc='PBE')], 

206 mode='pw') 

207 atoms.calc = calc 

208 atoms.get_potential_energy() 

209 assert np.allclose(atoms.get_forces(), 0, atol=1e-5) 

210 print(calc.dft.d3)