def hessian_component(self, index1, index2):
result = np.zeros((<NUM_LIT:3>, <NUM_LIT:3>), float)<EOL>if index1 == index2:<EOL><INDENT>for index3 in range(self.numc):<EOL><INDENT>if self.scaling[index1, index3] > <NUM_LIT:0>:<EOL><INDENT>d_1 = <NUM_LIT:1>/self.distances[index1, index3]<EOL>for (se, ve), (sg, vg), (sh, vh) in zip(<EOL>self.yield_pair_energies(index1, index3),<EOL>self.yield_pair_gradients(index1, index3),<EOL>self.yield_pair_hessians(index1, index3)<EOL>):<EOL><INDENT>result += (<EOL>+sh*self.dirouters[index1, index3]*ve<EOL>+sg*(np.identity(<NUM_LIT:3>, float) - self.dirouters[index1, index3])*ve*d_1<EOL>+sg*np.outer(self.directions[index1, index3],  vg)<EOL>+sg*np.outer(vg, self.directions[index1, index3])<EOL>+se*vh<EOL>)*self.scaling[index1, index3]<EOL><DEDENT><DEDENT><DEDENT><DEDENT>elif self.scaling[index1, index2] > <NUM_LIT:0>:<EOL><INDENT>d_1 = <NUM_LIT:1>/self.distances[index1, index2]<EOL>for (se, ve), (sg, vg), (sh, vh) in zip(<EOL>self.yield_pair_energies(index1, index2),<EOL>self.yield_pair_gradients(index1, index2),<EOL>self.yield_pair_hessians(index1, index2)<EOL>):<EOL><INDENT>result -= (<EOL>+sh*self.dirouters[index1, index2]*ve<EOL>+sg*(np.identity(<NUM_LIT:3>, float) - self.dirouters[index1, index2])*ve*d_1<EOL>+sg*np.outer(self.directions[index1, index2],  vg)<EOL>+sg*np.outer(vg, self.directions[index1, index2])<EOL>+se*vh<EOL>)*self.scaling[index1, index2]<EOL><DEDENT><DEDENT>return result<EOL>
Compute the hessian of the energy for one atom pair