Coverage for /builds/ase/ase/ase/utils/forcecurve.py: 78.63%

1from collections import namedtuple 

2 

3import numpy as np 

4 

5from ase.geometry import find_mic 

6 

7 

8def fit_raw(energies, forces, positions, cell=None, pbc=None): 

9 """Calculates parameters for fitting images to a band, as for 

10 a NEB plot.""" 

11 energies = np.array(energies) - energies[0] 

12 n_images = len(energies) 

13 fit_energies = np.empty((n_images - 1) * 20 + 1) 

14 fit_path = np.empty((n_images - 1) * 20 + 1) 

15 

16 path = [0] 

17 for i in range(n_images - 1): 

18 dR = positions[i + 1] - positions[i] 

19 if cell is not None and pbc is not None: 

20 dR, _ = find_mic(dR, cell, pbc) 

21 path.append(path[i] + np.sqrt((dR**2).sum())) 

22 

23 lines = [] # tangent lines 

24 lastslope = None 

25 for i in range(n_images): 

26 if i == 0: 

27 direction = positions[i + 1] - positions[i] 

28 dpath = 0.5 * path[1] 

29 elif i == n_images - 1: 

30 direction = positions[-1] - positions[-2] 

31 dpath = 0.5 * (path[-1] - path[-2]) 

32 else: 

33 direction = positions[i + 1] - positions[i - 1] 

34 dpath = 0.25 * (path[i + 1] - path[i - 1]) 

35 

36 direction /= np.linalg.norm(direction) 

37 slope = -(forces[i] * direction).sum() 

38 x = np.linspace(path[i] - dpath, path[i] + dpath, 3) 

39 y = energies[i] + slope * (x - path[i]) 

40 lines.append((x, y)) 

41 

42 if i > 0: 

43 s0 = path[i - 1] 

44 s1 = path[i] 

45 x = np.linspace(s0, s1, 20, endpoint=False) 

46 c = np.linalg.solve( 

47 np.array( 

48 [ 

49 (1, s0, s0**2, s0**3), 

50 (1, s1, s1**2, s1**3), 

51 (0, 1, 2 * s0, 3 * s0**2), 

52 (0, 1, 2 * s1, 3 * s1**2), 

53 ] 

54 ), 

55 np.array([energies[i - 1], energies[i], lastslope, slope]), 

56 ) 

57 y = c[0] + x * (c[1] + x * (c[2] + x * c[3])) 

58 fit_path[(i - 1) * 20 : i * 20] = x 

59 fit_energies[(i - 1) * 20 : i * 20] = y 

60 

61 lastslope = slope 

62 

63 fit_path[-1] = path[-1] 

64 fit_energies[-1] = energies[-1] 

65 return ForceFit(path, energies, fit_path, fit_energies, lines) 

66 

67 

68class ForceFit( 

69 namedtuple( 

70 'ForceFit', ['path', 'energies', 'fit_path', 'fit_energies', 'lines'] 

71 ) 

72): 

73 """Data container to hold fitting parameters for force curves.""" 

74 

75 def plot(self, ax=None): 

76 import matplotlib.pyplot as plt 

77 

78 if ax is None: 

79 ax = plt.gca() 

80 

81 ax.plot(self.path, self.energies, 'o') 

82 for x, y in self.lines: 

83 ax.plot(x, y, '-g') 

84 ax.plot(self.fit_path, self.fit_energies, 'k-') 

85 ax.set_xlabel(r'path [Å]') 

86 ax.set_ylabel('energy [eV]') 

87 Ef = max(self.energies) - self.energies[0] 

88 Er = max(self.energies) - self.energies[-1] 

89 dE = self.energies[-1] - self.energies[0] 

90 ax.set_title( 

91 r'$E_\mathrm{{f}} \approx$ {:.3f} eV; ' 

92 r'$E_\mathrm{{r}} \approx$ {:.3f} eV; ' 

93 r'$\Delta E$ = {:.3f} eV'.format(Ef, Er, dE) 

94 ) 

95 return ax 

96 

97 

98def fit_images(images): 

99 """Fits a series of images with a smoothed line for producing a standard 

100 NEB plot. Returns a `ForceFit` data structure; the plot can be produced 

101 by calling the `plot` method of `ForceFit`.""" 

102 R = [atoms.positions for atoms in images] 

103 E = [atoms.get_potential_energy() for atoms in images] 

104 F = [atoms.get_forces() for atoms in images] # XXX force consistent??? 

105 A = images[0].cell 

106 pbc = images[0].pbc 

107 return fit_raw(E, F, R, A, pbc) 

108 

109 

110def force_curve(images, ax=None): 

111 """Plot energies and forces as a function of accumulated displacements. 

112 

113 This is for testing whether a calculator's forces are consistent with 

114 the energies on a set of geometries where energies and forces are 

115 available.""" 

116 

117 if ax is None: 

118 import matplotlib.pyplot as plt 

119 

120 ax = plt.gca() 

121 

122 nim = len(images) 

123 

124 accumulated_distances = [] 

125 accumulated_distance = 0.0 

126 

127 # XXX force_consistent=True will work with some calculators, 

128 # but won't work if images were loaded from a trajectory. 

129 energies = [atoms.get_potential_energy() for atoms in images] 

130 

131 for i in range(nim): 

132 atoms = images[i] 

133 f_ac = atoms.get_forces() 

134 

135 if i < nim - 1: 

136 rightpos = images[i + 1].positions 

137 else: 

138 rightpos = atoms.positions 

139 

140 if i > 0: 

141 leftpos = images[i - 1].positions 

142 else: 

143 leftpos = atoms.positions 

144 

145 disp_ac, _ = find_mic( 

146 rightpos - leftpos, cell=atoms.cell, pbc=atoms.pbc 

147 ) 

148 

149 def total_displacement(disp): 

150 disp_a = (disp**2).sum(axis=1) ** 0.5 

151 return sum(disp_a) 

152 

153 dE_fdotr = -0.5 * np.vdot(f_ac.ravel(), disp_ac.ravel()) 

154 

155 linescale = 0.45 

156 

157 disp = 0.5 * total_displacement(disp_ac) 

158 

159 if i == 0 or i == nim - 1: 

160 disp *= 2 

161 dE_fdotr *= 2 

162 

163 x1 = accumulated_distance - disp * linescale 

164 x2 = accumulated_distance + disp * linescale 

165 y1 = energies[i] - dE_fdotr * linescale 

166 y2 = energies[i] + dE_fdotr * linescale 

167 

168 ax.plot([x1, x2], [y1, y2], 'b-') 

169 ax.plot(accumulated_distance, energies[i], 'bo') 

170 ax.set_ylabel('Energy [eV]') 

171 ax.set_xlabel('Accumulative distance [Å]') 

172 accumulated_distances.append(accumulated_distance) 

173 accumulated_distance += total_displacement(rightpos - atoms.positions) 

174 

175 ax.plot(accumulated_distances, energies, ':', zorder=-1, color='k') 

176 return ax 

177 

178 

179def plotfromfile(*fnames): 

180 from ase.io import read 

181 

182 nplots = len(fnames) 

183 

184 for i, fname in enumerate(fnames): 

185 images = read(fname, ':') 

186 import matplotlib.pyplot as plt 

187 

188 plt.subplot(nplots, 1, 1 + i) 

189 force_curve(images) 

190 plt.show() 

191 

192 

193if __name__ == '__main__': 

194 import sys 

195 

196 fnames = sys.argv[1:] 

197 plotfromfile(*fnames)