计算化学公社

标题: 跑完锂离子电池电解液体系的分子动力学后，怎么去分析锂离子在溶液中配位结构？ [打印本页]
作者
Author: Shuang-Jie    时间: 2023-6-16 11:18
标题: 跑完锂离子电池电解液体系的分子动力学后，怎么去分析锂离子在溶液中配位结构？
大家好，请教一下，用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).；


作者
Author: slxc920113    时间: 2023-6-16 13:32
本帖最后由 slxc920113 于 2025-11-4 11:32 编辑


  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

复制代码



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

  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()
复制代码


作者
Author: xjw    时间: 2024-4-22 13:39
您好，请问您成功解决了吗？我用楼下代码一直不成功，Log.txt文件没有内容
作者
Author: aff    时间: 2024-7-4 14:45
slxc920113 发表于 2023-6-16 13:32
需要你自己根据体系修改里面的内容，比如name Li，resname 溶剂的名称，配位原子，配位半径 等等。
用VM ...

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

也可以写成"O.*"，表示的原子名称是O开头的，后面跟着n个字符，\"是因为要做字符串的转义，select语句外面已经有引号了。通常的习惯是原子名称=元素名称+序号。如果你的拓扑文件里面不算这种规则，那你要根据你的实际情况来。
作者
Author: aff    时间: 2024-7-4 16:46
谢谢大佬回复！
我的体系是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文件吗？
作者
Author: slxc920113    时间: 2024-7-4 22:35
本帖最后由 slxc920113 于 2024-7-4 22:37 编辑
aff 发表于 2024-7-4 16:46
谢谢大佬回复！
我的体系是EC,EMC,NaPF6.  根据您的代码修改了下
set outfile [open Log.txt w]

\"O.*\"，需要加\转义，因为外面还有一层引号。
作者
Author: aff    时间: 2024-7-5 09:38
slxc920113 发表于 2024-7-4 22:35
\"O.*\"，需要加\转义，因为外面还有一层引号。

非常感谢大佬的回复！！！！！！  对我的体系其实只需要使用大佬的代码，改resname就可以了。 还需要请教一个问题，代码运行没有出问题，但是Log.txt文件却不知道在哪，请问这个文件存放在vmd的pwd的目录下吗？
作者
Author: slxc920113    时间: 2024-7-5 13:01
本帖最后由 slxc920113 于 2024-7-5 13:02 编辑
aff 发表于 2024-7-5 09:38
非常感谢大佬的回复！！！！！！  对我的体系其实只需要使用大佬的代码，改resname就可以了。 还需要请教 ...

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

PS：有个小技巧，你可以把快捷方式中起始目录内容删掉，这个时候无论你把快捷方式拷贝到哪里双击打开，哪里就是VMD的工作目录，和linux下一个效果。
作者
Author: Lemon97z    时间: 2024-7-26 16:57
本帖最后由 Lemon97z 于 2024-7-26 18:14 编辑
slxc920113 发表于 2023-6-16 13:32
需要你自己根据体系修改里面的内容，比如name Li，resname 溶剂的名称，配位原子，配位半径 等等。
用VM ...

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

上面写得很清楚，因为有些时候一个分子会有多个配位原子，如DME，所以输出的内容一行是周围原子的resname，一行是对应的resid，用于下一步的去重。
自己再写一个python小脚本去计这些团簇出现的次数和频率即可。
作者
Author: hesr    时间: 2024-10-6 02:04
slxc920113 发表于 2024-7-5 13:01
VMD的启动路径，具体看你的实际情况，如果是linux下应该就是当前工作目录，如果是windows，如果你是从快 ...

您好，跪求您加我15017930973，我尝试多次都失败了，请求您指导我一下
作者
Author: Hanna110    时间: 2025-9-29 20:02
slxc920113 发表于 2024-7-4 14:52
也可以写成"O.*"，表示的原子名称是O开头的，后面跟着n个字符，\"是因为要做字符串的转义，select语句外 ...

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



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3