本帖最后由 丁越 于 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}，则会在屏幕上输出所有的正则模式
- vmd > source cp2kNorm.tcl
- The number of atoms for vibration analysis: 9
- The total number of normal modes: 27
- Now is loading Normal modes....
- loading 1 of 27....
- loading 2 of 27....
- loading 3 of 27....
- loading 4 of 27....
- loading 5 of 27....
用法：norm [modeID 1] [scale 3] [radius 0.05] 如上面可见，此结构一共27个正则模式，如果你想看第五个正则模式，那么就输入norm 5，如果什么也不加直接输入norm则默认第一个正则模式; norm后面还可以接两个额外参数, scale和radius，scale控制箭头长短，默认值3；radius控制箭头粗细，默认值0.05。绘制效果如下所示：
2. 绘制红外光谱图
虽然这个脚本支持了这个功能，但起初的目的只是为了自己写写代码看看红外光谱是怎么计算出来的，不写的话自己就会被一直蒙在鼓里。如果有人想绘制红外光谱，则可以直接利用强大的Multiwfn程序，现在已经支持CP2K振动分析的输出文件并绘制红外，里面还有丰富的可调节参数。
用法： 直接输入大写的IR后就直接利用gnuplot将图绘制出来了，默认的FWHM是8cm^-1，基于洛伦兹函数展宽。如果想改为12 cm^-1, 则输入IR 12。输入IR后结果如下所示：
下面是脚本代码： - #Draw vibration modes of CP2K using VMD and plot IR spectra with gnuplot
- set filename "vib.out"
- #get the number of atoms for vibration analysis
- set fileid [open $filename r]
- set natm -1
- while {[gets $fileid line] >=0} {
- if {[string first "ATOM EL" $line] != -1} {break}
- }
- while {[string length $line] != 0} {
- gets $fileid line
- set natm [expr $natm+1]
- }
- seek $fileid 0 start
- puts "The number of atoms for vibratoin analysis: $natm"
- puts "The total number of normal modes: [expr $natm*3]"
- puts "Now is loading Normal modes...."
- #load Normal modes
- set nframe [expr 3*$natm]
- for {set i 1} {$i <= $nframe} {incr i} {
- while {[gets $fileid line] >=0} {
- if {[string first "ATOM EL" $line] != -1} {break}
- }
- set m1 [expr 3*$i-2]
- set m2 [expr 3*$i-1]
- set m3 [expr 3*$i]
- for {set iatm 1} {$iatm <=$natm} {incr iatm} {
- gets $fileid line
- scan $line "%d %s %f %f %f %f %f %f %f %f %f" itmp jtmp \
- X($iatm,$m1) Y($iatm,$m1) Z($iatm,$m1) \
- X($iatm,$m2) Y($iatm,$m2) Z($iatm,$m2) \
- X($iatm,$m3) Y($iatm,$m3) Z($iatm,$m3)
- lappend indexatm $itmp
- }
- puts "[format " loading %d of %d....\r" $i $nframe]"
- }
- #print Normal modes to screen by changing "if {0}" to "if {1}"
- if {0} {
- for {set i 1} {$i <=$nframe} {incr i} {
- puts "${i}_th Normal modes:"
- for {set j 1} {$j <=$natm} {incr j} {
- puts "[format "%-1.2f \t %1.2f \t %1.2f \n" $X($j,$i) $Y($j,$i) $Z($j,$i)]"
- }
- }
- }
- #load IR frequency and intensity into an array "a(freq)=intensities"
- seek $fileid 0 start
- for {set i 1} {$i <= $natm} {incr i} {
- while {[gets $fileid line] >=0} {
- if {[string first "VIB|Frequency" $line] != -1} {break}
- }
- lappend fq [lrange $line 2 4]
- while {[gets $fileid line] >=0} {
- if {[string first "VIB|Intensities" $line] != -1} {break}
- }
- lappend int [lrange $line 1 3]
- for {set t 0} {$t <3} {incr t} {
- foreach j [lindex $fq $t] k [lindex $int $t] {
- set a($j) $k
- }
- }
- unset fq
- unset int
- }
- close $fileid
- #Plot IR spectra by type command "IR", FWHM is 8 cm^-1 by default
- proc IR {{F 8}} {
- #Perform Lorentzian function broadening for isolated IR intensities
- global a
- set pi [expr acos(-1.0)]
- set nbin 10000
- set minfq [lindex [lsort -real -increasing [array names a]] 0]
- set maxfq [lindex [lsort -real -decreasing [array names a]] 0]
- set minrg [expr $minfq-$F*12.0]
- set maxrg [expr $maxfq+$F*12.0]
- set tau [expr 1.0*($maxrg-$minrg)/$nbin]
-
- for {set i 0} {$i<$nbin} {incr i} {
- set fqval [expr $minrg+$i*$tau]
- set L($fqval) 0
- }
-
- foreach j [lsort -real -increasing [array names a]] {
- for {set k 0} {$k<$nbin} {incr k} {
- set fqval [expr $minrg+$k*$tau]
- if {$fqval >= [expr $j-10.0*$F] && $fqval <= [expr $j+10.0*$F]} {
- set L($fqval) [expr $L($fqval)+100.0*$F*$a($j)/(2*$pi*(($fqval-$j)**2+0.25*${F}**2))]
- }
- }
- }
-
- #Export X-Y data set of curve to IRfile.dat
- set IRfileid [open IRfile.dat w]
- puts " IR data has been outputed into IRfile.dat, two columns are corresponding frequency and intensities,respectively"
- foreach i [lsort -real -increasing [array names L]] {
- puts $IRfileid "$i \t $L($i)"
- }
- close $IRfileid
-
- #Generate plot script for gnuplot
- set plotIRid [open plotIR.gp w]
- puts $plotIRid "set grid"
- puts $plotIRid "set xlabel 'Frequency (cm^-1)'"
- puts $plotIRid "set ylabel 'Molar adsorption coefficient (L*mol^-1*cm^-1)'"
- puts $plotIRid "set xrange \[0:$maxrg\]"
- puts $plotIRid "plot 'IRfile.dat' u 1:2 w l lw 1 not,"
- puts $plotIRid "pause -1"
- close $plotIRid
-
- #plot IR spectra with gnuplot
- exec gnuplot plotIR.gp
- }
- #Plot Normal modes
- proc norm {{imode 1} {scl 3} {rad 0.05}} {
- draw delete all
- draw color yellow
- global natm indexatm X Y Z
- #draw arrow
- for {set i 1} {$i<=$natm} {incr i} {
- #get each atom coordinate as cylinder center
- set sel [atomselect top "serial [lindex $indexatm [expr $i-1]]"]
- $sel update
- set cenx [$sel get x]
- set ceny [$sel get y]
- set cenz [$sel get z]
- #Scale vector
- set fragdx [expr $X($i,$imode)*$scl]
- set fragdy [expr $Y($i,$imode)*$scl]
- set fragdz [expr $Z($i,$imode)*$scl]
- #draw arrow
- set body 0.75
- set begx [expr $cenx-$fragdx/2]
- set begy [expr $ceny-$fragdy/2]
- set begz [expr $cenz-$fragdz/2]
- set endx [expr $cenx+$fragdx*$body-$fragdx/2]
- set endy [expr $ceny+$fragdy*$body-$fragdy/2]
- set endz [expr $cenz+$fragdz*$body-$fragdz/2]
- draw cylinder "$begx $begy $begz" "$endx $endy $endz" radius $rad filled yes resolution 20
- set endx2 [expr $cenx+$fragdx/2]
- set endy2 [expr $ceny+$fragdy/2]
- set endz2 [expr $cenz+$fragdz/2]
- draw cone "$endx $endy $endz" "$endx2 $endy2 $endz2" radius [expr $rad*2.5] resolution 20
- }
- }
| 
发表于 Post on 2022-3-14 18:35:11
收藏 Add to favorites
发表于 Post on 2022-3-15 14:18:01
楼主 Author
真没想到您居然是做实验的，您这shell写的比我可熟练太多了，我最近在用QEtoolkit.sh呢
特别喜欢代码运行或者程序得到结果所带来的快感。说实话学习计算化学特别感谢sob老师以及论坛的各路大神们，对新手入门帮助非常大。唉，人这一辈子转眼间几年十几年就过去了，我觉得找点自己喜欢干的事，自己能够乐在其中就足够了。