分子动力学模拟gro格式转换为 car

分子模拟中，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