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

1import numpy as np 

2 

3from ase.calculators.calculator import Calculator 

4from ase.utils import ff 

5 

6 

7class ForceField(Calculator): 

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

9 nolabel = True 

10 

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

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

13 Calculator.__init__(self, **kwargs) 

14 if (morses is None and 

15 bonds is None and 

16 angles is None and 

17 dihedrals is None and 

18 vdws is None and 

19 coulombs is None): 

20 raise ImportError( 

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

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

23 if morses is None: 

24 self.morses = [] 

25 else: 

26 self.morses = morses 

27 if bonds is None: 

28 self.bonds = [] 

29 else: 

30 self.bonds = bonds 

31 if angles is None: 

32 self.angles = [] 

33 else: 

34 self.angles = angles 

35 if dihedrals is None: 

36 self.dihedrals = [] 

37 else: 

38 self.dihedrals = dihedrals 

39 if vdws is None: 

40 self.vdws = [] 

41 else: 

42 self.vdws = vdws 

43 if coulombs is None: 

44 self.coulombs = [] 

45 else: 

46 self.coulombs = coulombs 

47 

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

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

50 if system_changes: 

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

52 self.results.pop(name, None) 

53 if 'energy' not in self.results: 

54 energy = 0.0 

55 for morse in self.morses: 

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

57 energy += e 

58 for bond in self.bonds: 

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

60 energy += e 

61 for angle in self.angles: 

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

63 energy += e 

64 for dihedral in self.dihedrals: 

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

66 atoms, dihedral) 

67 energy += e 

68 for vdw in self.vdws: 

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

70 energy += e 

71 for coulomb in self.coulombs: 

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

73 energy += e 

74 self.results['energy'] = energy 

75 if 'forces' not in self.results: 

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

77 for morse in self.morses: 

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

79 limits = get_limits([i, j]) 

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

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

82 for bond in self.bonds: 

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

84 limits = get_limits([i, j]) 

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

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

87 for angle in self.angles: 

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

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

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

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

92 for dihedral in self.dihedrals: 

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

94 atoms, dihedral) 

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

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

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

98 for vdw in self.vdws: 

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

100 limits = get_limits([i, j]) 

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

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

103 for coulomb in self.coulombs: 

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

105 limits = get_limits([i, j]) 

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

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

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

109 if 'hessian' not in self.results: 

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

111 for morse in self.morses: 

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

113 limits = get_limits([i, j]) 

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

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

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

117 for bond in self.bonds: 

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

119 limits = get_limits([i, j]) 

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

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

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

123 for angle in self.angles: 

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

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

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

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

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

129 for dihedral in self.dihedrals: 

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

131 atoms, dihedral) 

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

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

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

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

136 for vdw in self.vdws: 

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

138 limits = get_limits([i, j]) 

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

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

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

142 for coulomb in self.coulombs: 

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

144 limits = get_limits([i, j]) 

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

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

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

148 self.results['hessian'] = hessian 

149 

150 def get_hessian(self, atoms=None): 

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

152 

153 

154def get_limits(indices): 

155 gstarts = [] 

156 gstops = [] 

157 lstarts = [] 

158 lstops = [] 

159 for l, g in enumerate(indices): 

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

161 gstarts.append(g3) 

162 gstops.append(g3 + 3) 

163 lstarts.append(l3) 

164 lstops.append(l3 + 3) 

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