def _calc_subsystem_energies(self, ewald_matrix):
"""Modify the give matrix that consists of the eral and reciprocal sums
so that each entry x_ij is the full Ewald sum energy of a system
consisting of atoms i and j.
"""
q = self.q
# Create the self-term array where q1[i,j] is qi**2 + qj**2, except for
# the diagonal, where it is qi**2
q1 = q[None, :]**2 + q[:, None]**2
diag = np.diag(q1)/2
np.fill_diagonal(q1, diag)
q1_prefactor = -self.a/self.sqrt_pi
# Create the charge correction array where q2[i,j] is (qi + qj)**2,
# except for the diagonal where it is qi**2
q2 = q[None, :] + q[:, None]
q2 **= 2
diag = np.diag(q2)/4
np.fill_diagonal(q2, diag)
q2_prefactor = -np.pi/(2*self.volume*self.a_squared)
correction_matrix = q1_prefactor*q1 + q2_prefactor*q2
# Add the terms coming from x_ii and x_jj to the off-diagonal along
# with the corrections
n_atoms = self.n_atoms
final_matrix = np.zeros((n_atoms, n_atoms))
for i in range(n_atoms):
for j in range(n_atoms):
if i == j:
final_matrix[i, j] = ewald_matrix[i, j]
else:
pair_term = 2*ewald_matrix[i, j]
self_term_ii = ewald_matrix[i, i]
self_term_jj = ewald_matrix[j, j]
energy_total = pair_term + self_term_ii + self_term_jj
final_matrix[i, j] = energy_total
final_matrix += correction_matrix
return final_matrix
评论列表
文章目录
