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

skl

在一篇文献中看到gromacs参数优化的内容，其中涉及壳层模型库仑势、极化能以及分子偶极矩的计算，我自己用python写了一下，不知道对不对，请大家帮我看看，有不对的地方还请指教（原文献一并附在附件中）
库仑势及极化能计算：

偶极矩的方差计算：

python实现：

1. #定义CaF2分子，按顺序分别为F-、F-_s、F-、F-_s、Li+、Li+_s,其中，_s表示壳层
   
2. atoms = [{'beta': 9.139, 'q': 5.590, 'site': np.array([0, 1.869, -0.361])}, \
   
3.          {'beta': 9.139, 'q': -6.590, 'site': np.array([-0.23324962,  0.16647118,  0.29554485])}, \
   
4.          {'beta': 9.139, 'q': 5.590, 'site': np.array([0, -1.869, -0.361])}, \
   
5.          {'beta': 9.139, 'q': -6.590, 'site': np.array([-0.23324962,  0.16647118,  0.29554485])}, \
   
6.          {'beta': 13.772, 'q': 20, 'site': np.array([0., 0.    , 0.325    ])}, \
   
7.          {'beta': 13.772, 'q': -22, 'site': np.array([0.15000612, 0.0927579 , 0.29707778])}]
   
8. #计算库仑势和极化能
   
9. def cacl_ener(self,atoms):
   
10.     k = 1/(4*math.pi*0.08854187818)
    
11.     v_sum = 0
    
12.     #计算库仑势
    
13.     for i in range(0,5):
    
14.         atom = atom[1:]
    
15.         for j in atom:
    
16.             bi = atoms[i]['beta']
    
17.             bj = atoms[j]['beta']
    
18.             qi = atoms[i]['q']
    
19.             #if qi < 0:
    
20.             #    qi *= -1
    
21.             qj = atoms[j]['q']
    
22.             #if qj <0:
    
23.             #    qj *= -1
    
24.             r = np.linalg.norm(atoms[i]['site'] - atoms[j]['site'])
    
25.             bij = (bi * bj) /(math.sqrt(bi*bi+bj*bj))            
    
26.             erf= math.erf(bij * r)
    
27.             #print(qi,qj,erf,r)
    
28.             e_coul = k * qi * qj *erf /r
    
29.             v_sum += e_coul
    
30.    
    
31.     #计算极化能 a_f、a_ca为F-和Ca2+的极化率
    
32.     r_f_1 = np.linalg.norm(atoms[0]['site'] - atoms[1]['site'])
    
33.     r_f_2 = np.linalg.norm(atoms[2]['site'] - atoms[3]['site'])
    
34.     r_ca= np.linalg.norm(atoms[4]['site'] - atoms[5]['site'])
    
35.     v_pol_f_1 = k * (atoms[1]['q']**2/(2*self.a_f)) * r_f_1**2
    
36.     v_pol_f_2 = k * (atoms[3]['q']**2/(2*self.a_f)) * r_f_2**2
    
37.     v_pol_ca = k * (atoms[5]['q']**2/(2*self.a_ca)) * r_ca**2
    
38.     v_sum = v_sum + v_pol_f_1+ v_pol_ca   + v_pol_f_2
    
39.     return(v_sum)
    
40. #计算偶极矩的方差，u_exp为实验测得的偶极矩
    
41. def calc_x2(self,atoms):
    
42.     x2 = 0
    
43.     for atom in atoms:
    
44.         u = atom['q'] * atom['site']
    
45.         u = np.linalg.norm(u)
    
46.         x2 += (u -self.u_exp) ** 2
    
47.     x2 = (1/6) * x2
    
48.     return x2
    
49.    

复制代码

1. 
   

复制代码

赵小壮

shell model 模型啊，我刚写了一个。库伦势怎么只计算倒空间的呀，erfc(\alpha r)/r部分呢，偶极矩的计算有很多种，你照着你这个计算也行