def yield_pair_hessians(self, index1, index2):
d_1 = <NUM_LIT:1>/self.distances[index1, index2]<EOL>d_3 = d_1**<NUM_LIT:3><EOL>if self.charges is not None:<EOL><INDENT>c1 = self.charges[index1]<EOL>c2 = self.charges[index2]<EOL>yield <NUM_LIT:2>*c1*c2*d_3, np.zeros((<NUM_LIT:3>, <NUM_LIT:3>))<EOL><DEDENT>if self.dipoles is not None:<EOL><INDENT>d_5 = d_1**<NUM_LIT:5><EOL>d_7 = d_1**<NUM_LIT:7><EOL>p1 = self.dipoles[index1]<EOL>p2 = self.dipoles[index2]<EOL>yield <NUM_LIT:12>*d_5*np.dot(p1, p2), np.zeros((<NUM_LIT:3>, <NUM_LIT:3>))<EOL>yield -<NUM_LIT>*d_7, np.outer(p1, p2) + np.outer(p2, p1)<EOL>if self.charges is not None:<EOL><INDENT>yield <NUM_LIT:12>*c1*d_5, np.zeros((<NUM_LIT:3>, <NUM_LIT:3>))<EOL>yield <NUM_LIT:12>*c2*d_5, np.zeros((<NUM_LIT:3>, <NUM_LIT:3>))<EOL><DEDENT><DEDENT>
Yields pairs ((s''(r_ij), grad_i (x) grad_i v(bar{r}_ij))