请大家帮我看一下计算库伦势、极化能以及偶极矩的算法对不对 - 计算化学公社

#定义CaF2分子，按顺序分别为F-、F-_s、F-、F-_s、Li+、Li+_s,其中，_s表示壳层
atoms = [{'beta': 9.139, 'q': 5.590, 'site': np.array([0, 1.869, -0.361])}, \
{'beta': 9.139, 'q': -6.590, 'site': np.array([-0.23324962, 0.16647118, 0.29554485])}, \
{'beta': 9.139, 'q': 5.590, 'site': np.array([0, -1.869, -0.361])}, \
{'beta': 9.139, 'q': -6.590, 'site': np.array([-0.23324962, 0.16647118, 0.29554485])}, \
{'beta': 13.772, 'q': 20, 'site': np.array([0., 0. , 0.325 ])}, \
{'beta': 13.772, 'q': -22, 'site': np.array([0.15000612, 0.0927579 , 0.29707778])}]
#计算库仑势和极化能
def cacl_ener(self,atoms):
k = 1/(4*math.pi*0.08854187818)
v_sum = 0
#计算库仑势
for i in range(0,5):
atom = atom[1:]
for j in atom:
bi = atoms[i]['beta']
bj = atoms[j]['beta']
qi = atoms[i]['q']
#if qi < 0:
# qi *= -1
qj = atoms[j]['q']
#if qj <0:
# qj *= -1
r = np.linalg.norm(atoms[i]['site'] - atoms[j]['site'])
bij = (bi * bj) /(math.sqrt(bi*bi+bj*bj))
erf= math.erf(bij * r)
#print(qi,qj,erf,r)
e_coul = k * qi * qj *erf /r
v_sum += e_coul
#计算极化能 a_f、a_ca为F-和Ca2+的极化率
r_f_1 = np.linalg.norm(atoms[0]['site'] - atoms[1]['site'])
r_f_2 = np.linalg.norm(atoms[2]['site'] - atoms[3]['site'])
r_ca= np.linalg.norm(atoms[4]['site'] - atoms[5]['site'])
v_pol_f_1 = k * (atoms[1]['q']**2/(2*self.a_f)) * r_f_1**2
v_pol_f_2 = k * (atoms[3]['q']**2/(2*self.a_f)) * r_f_2**2
v_pol_ca = k * (atoms[5]['q']**2/(2*self.a_ca)) * r_ca**2
v_sum = v_sum + v_pol_f_1+ v_pol_ca + v_pol_f_2
return(v_sum)
#计算偶极矩的方差，u_exp为实验测得的偶极矩
def calc_x2(self,atoms):
x2 = 0
for atom in atoms:
u = atom['q'] * atom['site']
u = np.linalg.norm(u)
x2 += (u -self.u_exp) ** 2
x2 = (1/6) * x2
return x2

复制代码