跑完锂离子电池电解液体系的分子动力学后，怎么去分析锂离子在溶液中配位结构？

Shuang-Jie

大家好，请教一下，用gromacs跑完锂离子电池电解液体系的分子动力学后，怎么去分析锂离子在溶液中配位结构呢？去知道锂离子的溶剂化结构有哪些，以及不同结构的占比是什么样的呢？之前请教有前辈说写TCL脚本可以实现，大家能否给个类似的参照呢？


比如说文献的描述：
MD simulations were conducted to reveal the solvationstructure of Li+ in different mixtures of two solutions. In the 1mol L−1 LiTFSI-EC+DMC solution, the three most dominantsolvation structures are Li(EC)2(DMC)(TFSI)2, Li(EC)-(TFSI)4, and Li(EC)4(TFSI), with contents of 17.39%,8.30%, and 7.90%, respectively. Overall, the coordinationnumbers of Li+ with EC, DMC, and TFSI− are 2.37, 0.76, and2.77, respectively (Figure S4 and Table S2), while in 1 mol L−1LiTFSI-EC+EGD solution, the three most dominant solvationstructures are Li(EC)2(EGD)(TFSI)2 (12.74%), Li-(EC)3(EGD)(TFSI)2 (9.14%), and 1Li(EC)3(EGD)(TFSI)(8.76%); the respective coordination numbers with EC, EGD,and TFSI− are 2.28, 1.08, and 2.55 (Figure S7 and Table S3).；

slxc920113

1. 
   
2. 
   
3. set outfile [open Log.txt w]
   
4. set frames [molinfo top get numframes]
   
5. #resname Li 改成自己需要的中心原子
   
6. set center_atoms [atomselect top "resname Li"]
   
7. set center_ids [$center_atoms get resid]
   
8. set length [llength $center_ids]
   
9. # 10000 改成自己需要的采样数目
   
10. set step [expr $frames * $length / 10000 ]
    
11. $center_atoms delete
    
12. for { set frame 1 } { $frame < $frames } { incr frame} {
    
13.     for { set id 0 } { $id < $length } {incr id} {
    
14.         set current_resid [lindex $center_ids $id]
    
15.         # resname EC EMC 改成自己需要的配位分子名称，O改成自己需要的配位原子
    
16.         set nearby_atoms [atomselect top "resname EC EMC and name "O.*" and within 2.8 of resid $current_resid" frame $frame]
    
17.         set nearby_resids [lsort -unique [$nearby_atoms get resid]]
    
18.         $nearby_atoms delete
    
19.         set resid_list [join [lsort -unique [concat $nearby_resids $current_resid]] " "]
    
20.         set complete_atoms [atomselect top "resid $resid_list" frame $frame]
    
21.         puts $outfile [$complete_atoms get resname]
    
22.         puts $outfile [$complete_atoms get residue]
    
23.         # 会导出大量的溶剂化结构的pdb文件，建议在文件夹下执行脚本，不要在桌面，不需要导出具体的结构可以删除下面这一行
    
24.         $complete_atoms writepdb "${frame}_${id}.pdb"
    
25.         $complete_atoms delete
    
26.     }
    
27. }
    
28. close $outfile
    
29. 
    

复制代码

需要你自己根据体系修改里面的内容，比如name Li，resname 溶剂的名称，配位原子，配位半径 等等。
用VMD打开gro文件，然后载入xtc轨迹之后，在控制台执行脚本。
最后会统计所有Li的环境，写入Log.txt文件，因为有些时候一个分子会有多个配位原子，如DME，为了避免同一个分子多位点配位的时候被误识别成多个分子，所以输出的内容一行是resname，一行是对应的resid，用于下一步的去重。
自己再写一个python小脚本去计这些团簇出现的次数和频率即可。

1. 
   
2. out = open("Log.txt").read()
   
3. lines = out.strip().split('\n')
   
4. collect = {}
   
5. for i in range(0,len(lines),2):
   
6.     hash = {}
   
7.     resname = lines[i].split()
   
8.     resid = lines[i+1].split()
   
9.     if not len(resname) or not len(resid):
   
10.         continue
    
11.     thisid = resid[0]
    
12.     hash[resname[0]] = 1
    
13.     for i in range(1,len(resid)):
    
14.         if not resid[i] == thisid:
    
15.             name = resname[i]
    
16.             if name in hash:
    
17.                 hash[name] += 1
    
18.             else:
    
19.                 hash[name] = 1
    
20.         thisid = resid[i]
    
21.     keys = sorted(hash.keys())
    
22.     cluster = ''.join(map(lambda k:'%s(%s)'%(k,hash[k]),keys))
    
23.     if cluster in collect:
    
24.         collect[cluster] += 1
    
25.     else:
    
26.         collect[cluster] = 1
    
27. file = open('statistic.txt','w+')
    
28. items = list(collect.items())
    
29. items = sorted(items,key=lambda x:x[1],reverse=True)
    
30. for (k,v) in items:
    
31.     print(k,v)
    
32. file.writelines(list(map(lambda x:'%s:%s\n'%x,items)))
    
33. file.close()
    

复制代码

xjw

您好，请问您成功解决了吗？我用楼下代码一直不成功，Log.txt文件没有内容

aff

slxc920113 发表于 2023-6-16 13:32
需要你自己根据体系修改里面的内容，比如name Li，resname 溶剂的名称，配位原子，配位半径 等等。
用VM ...

请问大佬 \"O..\" \"O.\" \"F..\" \"F.\" 这个部分是什么含义呢？ 还想请问下 2.8是rdf后得到的第一溶剂层的距离吗？

slxc920113

aff 发表于 2024-7-4 14:45
请问大佬 \"O..\" \"O.\" \"F..\" \"F.\" 这个部分是什么含义呢？ 还想请问下 2.8是rdf后得到的第一溶剂 ...

也可以写成"O.*"，表示的原子名称是O开头的，后面跟着n个字符，\"是因为要做字符串的转义，select语句外面已经有引号了。通常的习惯是原子名称=元素名称+序号。如果你的拓扑文件里面不算这种规则，那你要根据你的实际情况来。

aff

谢谢大佬回复！
我的体系是EC,EMC,NaPF6.  根据您的代码修改了下
set outfile [open Log.txt w]
set frames [molinfo top get numframes]
set atomsLi [atomselect top "name Na1"]
set ids [$atomsNa get resid]
set length [llength $ids]
$atomsNa delete
for { set frame 1 } { $frame < $frames } { incr frame} {
    for { set id 0 } { $id < $length } {incr id} {
        set atoms [atomselect top "resid [lindex $ids $id] or same fragment as {resname EC EMC PF6 and name "O.*" "F.*" and within 3 of resid [lindex $ids $id] }" frame $frame]
        puts $outfile [$atoms get resname]
        puts $outfile [$atoms get resid]
        $atoms delete
    }
}
close $outfile

但在close $outfile 之前出现了 extra characters after close-quote。 并且在vmd的目录下也没有出现Log.txt。可以麻烦大佬看看是为什么输出不了Log.txt文件吗？

slxc920113

aff 发表于 2024-7-4 16:46
谢谢大佬回复！
我的体系是EC,EMC,NaPF6.  根据您的代码修改了下
set outfile [open Log.txt w]

\"O.*\"，需要加\转义，因为外面还有一层引号。

aff

slxc920113 发表于 2024-7-4 22:35
\"O.*\"，需要加\转义，因为外面还有一层引号。

非常感谢大佬的回复！！！！！！  对我的体系其实只需要使用大佬的代码，改resname就可以了。 还需要请教一个问题，代码运行没有出问题，但是Log.txt文件却不知道在哪，请问这个文件存放在vmd的pwd的目录下吗？

slxc920113

aff 发表于 2024-7-5 09:38
非常感谢大佬的回复！！！！！！  对我的体系其实只需要使用大佬的代码，改resname就可以了。 还需要请教 ...

VMD的启动路径，具体看你的实际情况，如果是linux下应该就是当前工作目录，如果是windows，如果你是从快捷方式打开的，可以看你的快捷方式里面的起始目录。如果直接进入到VMD的根目录下双击exe打开的，启动路径就是这个根目录。

PS：有个小技巧，你可以把快捷方式中起始目录内容删掉，这个时候无论你把快捷方式拷贝到哪里双击打开，哪里就是VMD的工作目录，和linux下一个效果。

Lemon97z

slxc920113 发表于 2023-6-16 13:32
需要你自己根据体系修改里面的内容，比如name Li，resname 溶剂的名称，配位原子，配位半径 等等。
用VM ...

您好 根据您的脚本 我现在得到了log.txt 但是单个溶剂重复出现23次  这是正常的吗？

slxc920113

Lemon97z 发表于 2024-7-26 16:57
您好 根据您的脚本 我现在得到了log.txt 但是单个溶剂重复出现23次  这是正常的吗？

上面写得很清楚，因为有些时候一个分子会有多个配位原子，如DME，所以输出的内容一行是周围原子的resname，一行是对应的resid，用于下一步的去重。
自己再写一个python小脚本去计这些团簇出现的次数和频率即可。

hesr

slxc920113 发表于 2024-7-5 13:01
VMD的启动路径，具体看你的实际情况，如果是linux下应该就是当前工作目录，如果是windows，如果你是从快 ...

您好，跪求您加我15017930973，我尝试多次都失败了，请求您指导我一下

Hanna110

slxc920113 发表于 2024-7-4 14:52
也可以写成"O.*"，表示的原子名称是O开头的，后面跟着n个字符，\"是因为要做字符串的转义，select语句外 ...

O..和O.有什么区别呢？