Coverage for /builds/ase/ase/ase/calculators/ff.py: 78.57%

1# fmt: off 

2 

3import numpy as np 

4 

5from ase.calculators.calculator import Calculator 

6from ase.utils import ff 

7 

8 

9class ForceField(Calculator): 

10 implemented_properties = ['energy', 'forces'] 

11 nolabel = True 

12 

13 def __init__(self, morses=None, bonds=None, angles=None, dihedrals=None, 

14 vdws=None, coulombs=None, **kwargs): 

15 Calculator.__init__(self, **kwargs) 

16 if (morses is None and 

17 bonds is None and 

18 angles is None and 

19 dihedrals is None and 

20 vdws is None and 

21 coulombs is None): 

22 raise ImportError( 

23 "At least one of morses, bonds, angles, dihedrals," 

24 "vdws or coulombs lists must be defined!") 

25 if morses is None: 

26 self.morses = [] 

27 else: 

28 self.morses = morses 

29 if bonds is None: 

30 self.bonds = [] 

31 else: 

32 self.bonds = bonds 

33 if angles is None: 

34 self.angles = [] 

35 else: 

36 self.angles = angles 

37 if dihedrals is None: 

38 self.dihedrals = [] 

39 else: 

40 self.dihedrals = dihedrals 

41 if vdws is None: 

42 self.vdws = [] 

43 else: 

44 self.vdws = vdws 

45 if coulombs is None: 

46 self.coulombs = [] 

47 else: 

48 self.coulombs = coulombs 

49 

50 def calculate(self, atoms, properties, system_changes): 

51 Calculator.calculate(self, atoms, properties, system_changes) 

52 if system_changes: 

53 for name in ['energy', 'forces', 'hessian']: 

54 self.results.pop(name, None) 

55 if 'energy' not in self.results: 

56 energy = 0.0 

57 for morse in self.morses: 

58 i, j, e = ff.get_morse_potential_value(atoms, morse) 

59 energy += e 

60 for bond in self.bonds: 

61 i, j, e = ff.get_bond_potential_value(atoms, bond) 

62 energy += e 

63 for angle in self.angles: 

64 i, j, k, e = ff.get_angle_potential_value(atoms, angle) 

65 energy += e 

66 for dihedral in self.dihedrals: 

67 i, j, k, l, e = ff.get_dihedral_potential_value( 

68 atoms, dihedral) 

69 energy += e 

70 for vdw in self.vdws: 

71 i, j, e = ff.get_vdw_potential_value(atoms, vdw) 

72 energy += e 

73 for coulomb in self.coulombs: 

74 i, j, e = ff.get_coulomb_potential_value(atoms, coulomb) 

75 energy += e 

76 self.results['energy'] = energy 

77 if 'forces' not in self.results: 

78 forces = np.zeros(3 * len(atoms)) 

79 for morse in self.morses: 

80 i, j, g = ff.get_morse_potential_gradient(atoms, morse) 

81 limits = get_limits([i, j]) 

82 for gb, ge, lb, le in limits: 

83 forces[gb:ge] -= g[lb:le] 

84 for bond in self.bonds: 

85 i, j, g = ff.get_bond_potential_gradient(atoms, bond) 

86 limits = get_limits([i, j]) 

87 for gb, ge, lb, le in limits: 

88 forces[gb:ge] -= g[lb:le] 

89 for angle in self.angles: 

90 i, j, k, g = ff.get_angle_potential_gradient(atoms, angle) 

91 limits = get_limits([i, j, k]) 

92 for gb, ge, lb, le in limits: 

93 forces[gb:ge] -= g[lb:le] 

94 for dihedral in self.dihedrals: 

95 i, j, k, l, g = ff.get_dihedral_potential_gradient( 

96 atoms, dihedral) 

97 limits = get_limits([i, j, k, l]) 

98 for gb, ge, lb, le in limits: 

99 forces[gb:ge] -= g[lb:le] 

100 for vdw in self.vdws: 

101 i, j, g = ff.get_vdw_potential_gradient(atoms, vdw) 

102 limits = get_limits([i, j]) 

103 for gb, ge, lb, le in limits: 

104 forces[gb:ge] -= g[lb:le] 

105 for coulomb in self.coulombs: 

106 i, j, g = ff.get_coulomb_potential_gradient(atoms, coulomb) 

107 limits = get_limits([i, j]) 

108 for gb, ge, lb, le in limits: 

109 forces[gb:ge] -= g[lb:le] 

110 self.results['forces'] = np.reshape(forces, (len(atoms), 3)) 

111 if 'hessian' not in self.results: 

112 hessian = np.zeros((3 * len(atoms), 3 * len(atoms))) 

113 for morse in self.morses: 

114 i, j, h = ff.get_morse_potential_hessian(atoms, morse) 

115 limits = get_limits([i, j]) 

116 for gb1, ge1, lb1, le1 in limits: 

117 for gb2, ge2, lb2, le2 in limits: 

118 hessian[gb1:ge1, gb2:ge2] += h[lb1:le1, lb2:le2] 

119 for bond in self.bonds: 

120 i, j, h = ff.get_bond_potential_hessian(atoms, bond) 

121 limits = get_limits([i, j]) 

122 for gb1, ge1, lb1, le1 in limits: 

123 for gb2, ge2, lb2, le2 in limits: 

124 hessian[gb1:ge1, gb2:ge2] += h[lb1:le1, lb2:le2] 

125 for angle in self.angles: 

126 i, j, k, h = ff.get_angle_potential_hessian(atoms, angle) 

127 limits = get_limits([i, j, k]) 

128 for gb1, ge1, lb1, le1 in limits: 

129 for gb2, ge2, lb2, le2 in limits: 

130 hessian[gb1:ge1, gb2:ge2] += h[lb1:le1, lb2:le2] 

131 for dihedral in self.dihedrals: 

132 i, j, k, l, h = ff.get_dihedral_potential_hessian( 

133 atoms, dihedral) 

134 limits = get_limits([i, j, k, l]) 

135 for gb1, ge1, lb1, le1 in limits: 

136 for gb2, ge2, lb2, le2 in limits: 

137 hessian[gb1:ge1, gb2:ge2] += h[lb1:le1, lb2:le2] 

138 for vdw in self.vdws: 

139 i, j, h = ff.get_vdw_potential_hessian(atoms, vdw) 

140 limits = get_limits([i, j]) 

141 for gb1, ge1, lb1, le1 in limits: 

142 for gb2, ge2, lb2, le2 in limits: 

143 hessian[gb1:ge1, gb2:ge2] += h[lb1:le1, lb2:le2] 

144 for coulomb in self.coulombs: 

145 i, j, h = ff.get_coulomb_potential_hessian(atoms, coulomb) 

146 limits = get_limits([i, j]) 

147 for gb1, ge1, lb1, le1 in limits: 

148 for gb2, ge2, lb2, le2 in limits: 

149 hessian[gb1:ge1, gb2:ge2] += h[lb1:le1, lb2:le2] 

150 self.results['hessian'] = hessian 

151 

152 def get_hessian(self, atoms=None): 

153 return self.get_property('hessian', atoms) 

154 

155 

156def get_limits(indices): 

157 gstarts = [] 

158 gstops = [] 

159 lstarts = [] 

160 lstops = [] 

161 for l, g in enumerate(indices): 

162 g3, l3 = 3 * g, 3 * l 

163 gstarts.append(g3) 

164 gstops.append(g3 + 3) 

165 lstarts.append(l3) 

166 lstops.append(l3 + 3) 

167 return zip(gstarts, gstops, lstarts, lstops)