计算化学公社

 找回密码 Forget password
 注册 Register
Views: 3359|回复 Reply: 6

[CP2K] 在VMD中绘制CP2K计算产生的振动矢量的方法

[复制链接 Copy URL]

327

帖子

9

威望

1995

eV
积分
2502

Level 5 (御坂)

发表于 Post on 2022-3-14 18:35:11 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 丁越 于 2022-3-14 22:49 编辑

在VMD中绘制CP2K计算产生的振动矢量的方法



一、前言

  sob大大之前写了《在VMD中绘制Gaussian计算的分子振动矢量的方法》(http://bbs.keinsci.com/thread-19287-1-1.html),于是在我在这个脚本的基础上魔改了一番,使得其能很方便绘制CP2K振动分析产生的分子振动矢量并且可以绘制红外光谱。另外,目前就我了解到的一些程序能支持CP2K相关的振动分析文件并可视化振动模式的有:

  * molden: 我仅仅用过windows版的,总之体验很不友好,首先就载入文件这一项让人很头疼,一时间找不到那个框框在哪里(它最小化缩到了桌面右下角,费半天劲才找到。。。)。

  * jmol:非常流行的一款可视化程序,完美支持.mol文件并可以观看振动动画。但是唯一有个缺点是振动矢量的箭头不美观,而且太小了。

  *PyVibMS: 暂且直接支持绘制由CP2K产生的振动分析输出文件,但是可以写脚本提取振动矢量信息,然后制作成这个程序的输入文件。

二、脚本使用

  1. 绘制振动矢量
   首先将结构文件载入vmd,可以是如上步你做结构优化的产生的xxx-pos-1.xyz文件、pdb、mol2等。接着,将该脚本cp2kNorm.tcl放到CP2K振动分析输出文件的目录下,然后将脚本中set filename "vib.out"中的vib.out改为你计算中所用的文件名,最后source cp2kNorm.tcl这个脚本。此时, 命令控制台上就输出了你做振动分析的原子数,产生了多少个正则模式,正则模式载入的情况,如下所示。另外,如果将脚本中第41行的if {0} 改为if {1},则会在屏幕上输出所有的正则模式

  1. vmd > source cp2kNorm.tcl                           
  2. The number of atoms for vibration analysis: 9      
  3. The total number of normal modes: 27               
  4. Now is loading Normal modes....                     
  5.     loading 1 of 27....                             
  6.     loading 2 of 27....                             
  7.     loading 3 of 27....                             
  8.     loading 4 of 27....                             
  9.     loading 5 of 27....                 
复制代码
  
  用法:norm [modeID 1] [scale 3] [radius 0.05]
  如上面可见,此结构一共27个正则模式,如果你想看第五个正则模式,那么就输入norm 5,如果什么也不加直接输入norm则默认第一个正则模式; norm后面还可以接两个额外参数, scale和radius,scale控制箭头长短,默认值3;radius控制箭头粗细,默认值0.05。绘制效果如下所示:

norm.png



  2. 绘制红外光谱图   
  虽然这个脚本支持了这个功能,但起初的目的只是为了自己写写代码看看红外光谱是怎么计算出来的,不写的话自己就会被一直蒙在鼓里。如果有人想绘制红外光谱,则可以直接利用强大的Multiwfn程序,现在已经支持CP2K振动分析的输出文件并绘制红外,里面还有丰富的可调节参数。
  用法: 直接输入大写的IR后就直接利用gnuplot将图绘制出来了,默认的FWHM是8cm^-1,基于洛伦兹函数展宽。如果想改为12 cm^-1, 则输入IR 12。输入IR后结果如下所示:
IR.png


  下面是脚本代码:
  1. #Draw vibration modes of CP2K using VMD and plot IR spectra with gnuplot

  2. set filename "vib.out"

  3. #get the number of atoms for vibration analysis
  4. set fileid [open $filename r]
  5. set natm -1
  6. while {[gets $fileid line] >=0} {
  7.                 if {[string first "ATOM  EL" $line] != -1} {break}
  8. }
  9. while {[string length $line] != 0} {
  10.         gets $fileid line
  11.         set natm [expr $natm+1]
  12. }
  13. seek $fileid 0 start
  14. puts "The number of atoms for vibratoin analysis: $natm"
  15. puts "The total number of normal modes: [expr $natm*3]"
  16. puts "Now is loading Normal modes...."

  17. #load Normal modes
  18. set nframe [expr 3*$natm]
  19. for {set i 1} {$i <= $nframe} {incr i} {
  20.         while {[gets $fileid line] >=0} {
  21.                 if {[string first "ATOM  EL" $line] != -1} {break}
  22.         }
  23.         set m1 [expr 3*$i-2]
  24.         set m2 [expr 3*$i-1]
  25.         set m3 [expr 3*$i]
  26.         for {set iatm 1} {$iatm <=$natm} {incr iatm} {
  27.                 gets $fileid line
  28.                 scan $line "%d %s %f %f %f %f %f %f %f %f %f" itmp jtmp \
  29.                 X($iatm,$m1) Y($iatm,$m1) Z($iatm,$m1) \
  30.                 X($iatm,$m2) Y($iatm,$m2) Z($iatm,$m2) \
  31.                 X($iatm,$m3) Y($iatm,$m3) Z($iatm,$m3)
  32.                 lappend indexatm $itmp
  33.         }
  34.         puts "[format "    loading %d of %d....\r" $i $nframe]"
  35. }

  36. #print Normal modes to screen by changing "if {0}" to "if {1}"
  37. if {0} {
  38.         for {set i 1} {$i <=$nframe} {incr i} {
  39.                 puts "${i}_th Normal modes:"
  40.                 for {set j 1} {$j <=$natm} {incr j} {
  41.                         puts "[format "%-1.2f \t %1.2f \t %1.2f \n" $X($j,$i) $Y($j,$i) $Z($j,$i)]"
  42.                 }
  43.         }
  44. }

  45. #load IR frequency and intensity into an array "a(freq)=intensities"
  46. seek $fileid 0 start
  47. for {set i 1} {$i <= $natm} {incr i} {
  48.         while {[gets $fileid line] >=0} {
  49.                 if {[string first "VIB|Frequency" $line] != -1} {break}
  50.         }        
  51.         lappend fq [lrange $line 2 4]        
  52.         while {[gets $fileid line] >=0} {
  53.                 if {[string first "VIB|Intensities" $line] != -1} {break}
  54.         }               
  55.         lappend int [lrange $line 1 3]
  56.         for {set t 0} {$t <3} {incr t} {
  57.                 foreach j [lindex $fq $t] k [lindex $int $t] {
  58.                         set a($j) $k
  59.                 }
  60.         }        
  61.         unset fq
  62.         unset int
  63. }                        
  64. close $fileid

  65. #Plot IR spectra by type command "IR", FWHM is 8 cm^-1 by default
  66. proc IR {{F 8}} {
  67.         #Perform Lorentzian function broadening for isolated IR intensities
  68.         global a
  69.         set pi [expr acos(-1.0)]
  70.         set nbin 10000
  71.         set minfq [lindex [lsort -real -increasing [array names a]] 0]
  72.         set maxfq [lindex [lsort -real -decreasing [array names a]] 0]
  73.         set minrg [expr $minfq-$F*12.0]
  74.         set maxrg [expr $maxfq+$F*12.0]
  75.         set tau [expr 1.0*($maxrg-$minrg)/$nbin]
  76.         
  77.         for {set i 0} {$i<$nbin} {incr i} {
  78.                 set fqval [expr $minrg+$i*$tau]
  79.                 set L($fqval) 0
  80.         }               
  81.                         
  82.         foreach j [lsort -real -increasing [array names a]] {
  83.                 for {set k 0} {$k<$nbin} {incr k} {
  84.                         set fqval [expr $minrg+$k*$tau]
  85.                         if {$fqval >= [expr $j-10.0*$F] && $fqval <= [expr $j+10.0*$F]} {
  86.                                 set L($fqval) [expr $L($fqval)+100.0*$F*$a($j)/(2*$pi*(($fqval-$j)**2+0.25*${F}**2))]
  87.                         }
  88.                 }
  89.         }
  90.         
  91.         #Export X-Y data set of curve to IRfile.dat
  92.         set IRfileid [open IRfile.dat w]
  93.         puts " IR data has been outputed into IRfile.dat, two columns are corresponding frequency and intensities,respectively"
  94.         foreach i [lsort -real -increasing [array names L]] {
  95.                 puts $IRfileid "$i \t $L($i)"
  96.         }
  97.         close $IRfileid
  98.         
  99.         #Generate plot script for gnuplot
  100.         set plotIRid [open plotIR.gp w]
  101.         puts $plotIRid "set grid"
  102.         puts $plotIRid "set xlabel 'Frequency (cm^-1)'"
  103.         puts $plotIRid "set ylabel 'Molar adsorption coefficient (L*mol^-1*cm^-1)'"
  104.         puts $plotIRid "set xrange \[0:$maxrg\]"
  105.         puts $plotIRid "plot 'IRfile.dat' u 1:2 w l lw 1 not,"
  106.         puts $plotIRid "pause -1"
  107.         close $plotIRid
  108.         
  109.         #plot IR spectra with gnuplot
  110.         exec gnuplot plotIR.gp
  111. }

  112. #Plot Normal modes
  113. proc norm {{imode 1} {scl 3} {rad 0.05}} {
  114.         draw delete all
  115.         draw color yellow
  116.         global natm indexatm X Y Z
  117.         #draw arrow
  118.         for {set i 1} {$i<=$natm} {incr i} {
  119.                 #get each atom coordinate as cylinder center
  120.                 set sel [atomselect top "serial [lindex $indexatm [expr $i-1]]"]
  121.                 $sel update
  122.                 set cenx [$sel get x]
  123.                 set ceny [$sel get y]
  124.                 set cenz [$sel get z]
  125.                 #Scale vector  
  126.                 set fragdx [expr $X($i,$imode)*$scl]
  127.                 set fragdy [expr $Y($i,$imode)*$scl]
  128.                 set fragdz [expr $Z($i,$imode)*$scl]
  129.                 #draw arrow
  130.                 set body 0.75
  131.                 set begx [expr $cenx-$fragdx/2]
  132.                 set begy [expr $ceny-$fragdy/2]
  133.                 set begz [expr $cenz-$fragdz/2]
  134.                 set endx [expr $cenx+$fragdx*$body-$fragdx/2]
  135.                 set endy [expr $ceny+$fragdy*$body-$fragdy/2]
  136.                 set endz [expr $cenz+$fragdz*$body-$fragdz/2]
  137.                 draw cylinder "$begx $begy $begz" "$endx $endy $endz" radius $rad filled yes resolution 20
  138.                 set endx2 [expr $cenx+$fragdx/2]
  139.                 set endy2 [expr $ceny+$fragdy/2]
  140.                 set endz2 [expr $cenz+$fragdz/2]
  141.                 draw cone "$endx $endy $endz" "$endx2 $endy2 $endz2" radius [expr $rad*2.5] resolution 20
  142.         }
  143. }

复制代码

  最后请教下坛子的各位老师:本来打算写一个能直接播放振动动画的,但是没想明白这方面该怎么实现,如果有思路的话还请指点下



cp2kNorm.tcl

4.63 KB, 下载次数 Times of downloads: 8

评分 Rate

参与人数
Participants 6
威望 +1 eV +25 收起 理由
Reason
djjj148 + 5 好物!
卡开发发 + 5 欢迎讨论
sobereva + 1
ChemG + 5 谢谢分享
hdhxx123 + 5 好物!
wolfli369 + 5

查看全部评分 View all ratings

自由发挥,野蛮生长

4万

帖子

99

威望

4万

eV
积分
89975

管理员

公社社长+计算化学玩家

发表于 Post on 2022-3-15 14:18:01 | 显示全部楼层 Show all
molden倒也能用,但是界面设计太过时了,显示效果也很丑。但凡能有其它像样的程序能替代,我估计大多数人都不想用那玩意

振动动画不是很好实现,一个可行的办法是根据振动矢量和当前结构产生振动过程的各个点的坐标作为一个xyz轨迹文件,自动让VMD载入,并让让VMD播放。
另一个做法是魔改一下我的OfakeG,弄个CP2KfakeG,把CP2K振动分析输出文件弄成Gaussian的格式,让GaussView显示动画
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班,内容介绍以及往届资料购买请点击链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的最佳途径。培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

327

帖子

9

威望

1995

eV
积分
2502

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2022-3-15 14:34:19 | 显示全部楼层 Show all
sobereva 发表于 2022-3-15 14:18
molden倒也能用,但是界面设计太过时了,显示效果也很丑。但凡能有其它像样的程序能替代,我估计大多数人都 ...

谢谢sob老师,既然不好搞的话看振动动画就用jmol就行,这个比molden好用太多
自由发挥,野蛮生长

3015

帖子

3

威望

1万

eV
积分
14751

Level 6 (一方通行)

从头算界孔乙己

发表于 Post on 2022-3-15 18:10:22 | 显示全部楼层 Show all
丁越 发表于 2022-3-15 14:34
谢谢sob老师,既然不好搞的话看振动动画就用jmol就行,这个比molden好用太多

其实产生轨迹可以用ase试试,利用neb那个模块的插点可以形成轨迹,正反两个过程拼一下再转*.xyz轨迹就行了,当然方法不是唯一的。显示箭头的方法除了vmd,ase也可以,将力或者速度做替换理论上用ase显示,其他还可以使用vesta,不过这个就得自己多写一些东西了。总而言之,自己动手总比伸手强。
近期不及时回复。欢迎无偏见非商业的学术讨论,但是看家本领和课题组的传统艺能别人会毫无保留告诉你?

327

帖子

9

威望

1995

eV
积分
2502

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2022-3-15 18:52:35 | 显示全部楼层 Show all
卡开发发 发表于 2022-3-15 18:10
其实产生轨迹可以用ase试试,利用neb那个模块的插点可以形成轨迹,正反两个过程拼一下再转*.xyz轨迹就行 ...

谢谢卡开发发老师提供思路,学习python和ASE之类的还在N个月后的计划中。。目前还不能给自己揽一大堆要学的,怕是现在手头上需要学的都学不懂,不过过段时间可以回头再搞

评分 Rate

参与人数
Participants 1
eV +1 收起 理由
Reason
卡开发发 + 1 欢迎讨论

查看全部评分 View all ratings

自由发挥,野蛮生长

287

帖子

8

威望

1664

eV
积分
2111

Level 5 (御坂)

发表于 Post on 2022-3-16 01:07:04 | 显示全部楼层 Show all
有兴趣的话可以看看能不能把cp2k的支持加入PyVibMS 我给你merge进去

327

帖子

9

威望

1995

eV
积分
2502

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2022-3-16 08:34:05 | 显示全部楼层 Show all
本帖最后由 丁越 于 2022-3-16 08:39 编辑
smutao 发表于 2022-3-16 01:07
有兴趣的话可以看看能不能把cp2k的支持加入PyVibMS 我给你merge进去

太好了老师,希望PyVibMS更能拓展一下功能,搞一个CP2K的GUI
老师的意思是让我写一下吗(不知道我有没有误解你的话)?要真这样的话那可就太有点为难我了,我还是菜鸟初学者,编程只会点shell和TCL
自由发挥,野蛮生长

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-7 02:25 , Processed in 0.678075 second(s), 26 queries .

快速回复 返回顶部 返回列表 Return to list