分子模拟中,Materials Studio 跑分子动力学很慢,效率是gromacs的1/20。为了提高效率,可先在gromacs上运行分子动力学,获得平衡结构。再将gro文件转换为car文件并导入MS软件,进行分析,譬如XRD。car文件可保持周期性体系在ms中打开,也可在其它分子模拟软件中打开。(openbabel软件可以轻松实现gro2cif,但是没有电荷。)
car文件主要参数是坐标和电荷。
根据rtp文件中的atoms,编写一个charge()函数。这个可用python解决,代码如下:
#!/bin/env python
# -*- encoding: utf-8 -*-
f1 = open('atoms.txt','r',encoding='utf-8') #打开txt文本进行读取
f2 = open('charge.py','w')
f2.write("def chargefun(list0):\n")
f2.write(" if list0=='': print('There is at least one mistake!');\n")
while True: #循环,读取文本里面的所有内容
line = f1.readline() #一行一行读取
if not line: #如果没有内容,则退出循环
break
list = line.strip().split()
a=list[0]
b=list[2]
f2.write(" elif list0=='%s': charge='%s';\n"%(a,b))
f2.write(" else: print('%s There is at least one mistake!'%list0);\n")
f2.write(" return charge\n")
f1.close()
f2.close()
将charge()复制到将要编写的gro2car文件中,也可以通过from...import...引用。完成之后的代码如下:
#!/bin/env python
# -*- encoding: utf-8 -*-
import linecache
cartype="{0:<8}{1:<16}{2:15}{3:12}{4:5}{5:7}{6:8}{7:4}{8:12}" #car格式
'''
本例中有2个残基,30000+原子数,当原子编号大于9999时,gro文件的排版会有点小问题。
原子命名要求以元素名作为首字母,如果元素名多于1个字母,请自行修改代码。
'''
filename = 'as-c6' #文件名
asname = 'ASC6' #第一个残基名称
asmonum = 44 #第一个残基分子数
asatnum = 103 #第一个残基原子数
xyname = 'XYL' #第二个残基名称
xymonum = 1508 #第二个残基分子数
xyatnum = 18 #第二个残基原子数
f=open(filename+'.car','w')
f.write('!BIOSYM archive 3\n')
f.write('PBC=ON\n\n')
line0=linecache.getline(filename+'.gro',asmonum*asatnum+xymonum*xyatnum+3)
list1=line0.strip().split()
f.write('PBC ')
f.write("%.4f"%(float(list1[0])*10)+' ')
f.write("%.4f"%(float(list1[1])*10)+' ')
f.write("%.4f"%(float(list1[2])*10)+' ')
f.write(' 90 90 90 (P1)\n')
def chargefun(list0):
if list0=='': print('There is at least one mistake!');
elif list0=='C00': charge='0.192834809';
elif list0=='C01': charge='-0.215146234';
elif list0=='C02': charge='-0.215146234';
elif list0=='C03': charge='0.192834809';
elif list0=='C04': charge