|
|
用AIMD跑分子动力学,想用python脚本处理,显示MD过程中键长的变化的程序MD-bond.py,但是我们复制脚本进去运行,一直提示缩进错误,请高手指点。以下是我复制的脚本内容
import numpy as np
# input atoms No.(b1,b2) bwtween which you want to calculate the bond
# input the MD-steps (Mnum) in your VASP MD calculation
b1= 12
b2= 27
Mnum=10000
b1=b1-1
b2=b2-1
a=open("XDATCAR",'r+')
b=open("bond.dat",'w+')
#skip line1-2
for i in range(2):
a.readline()
#get vectors
X=a.readline().strip().split()
Y=a.readline().strip().split()
Z=a.readline().strip().split()
X=[float(X) for i in range(3)]
Y=[float(Y) for i in range(3)]
Z=[float(Z) for i in range(3)]
#skip line-6
a.readline()
#get total atom number
num=a.readline().strip().split()
num=[int(num) for i in range(len(num))]
num=sum(num)
#define function of bond-length
def bond(m,n):
r=[n-m for i in range(3)]
r2=[r[0]*X+r[1]*Y+r[2]*Z for i in range(3)]
r3=np.sqrt(r2[0]**2+r2[1]**2+r2[2]**2)
return r3
#get bond-length of each step
def bond_s():
a.readline() #skip the line "Direct configuration="
coord=[]
for i in range(num):
r=a.readline().strip().split()
r=[float(r) for i in range(3)]
coord.append(r)
C_chose=[coord[b1],coord[b2]]
return C_chose
for j in range(Mnum):
m,n=bond_s()
r=bond(m,n)
b.write(str(r)+"n")
|
|