r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import numpy as np

3from ase.calculators.calculator import Calculator

4from ase.utils import ff

7class ForceField(Calculator):

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

9 nolabel = True

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

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

150 def get_hessian(self, atoms=None):

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

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)