使用NewtonX和Gaussian16产生激发态谐振分布（更新genHOD脚本）

ggdh

模拟吸收发射光谱的方法sob写了不少文章，可以先参考：
Gaussian中用TDDFT计算激发态和吸收、荧光、磷光光谱的方法
振动分辨的电子光谱的计算
使用Multiwfn绘制构象权重平均的光谱
NewtonX产生吸收发射光谱的原理是基于预先算好分子的正则振动模随机产生一批符合谐振分布（harmonic oscillator distribution）的分子结构。然后在这些结构基础之上计算垂直跃迁。最后将所有结构的垂直跃迁加和，得到相应的光谱图。
NewtonX官网提供了Tutorial 和 Docs，不过里面没有模拟荧光的教程，而且他的操作比较麻烦，这里我就来做一个简明版的教程。

前言：
这是本算例中计算的分子DMABN，

这个分子在极性溶剂中由于TICT现象，弛豫激发态中的NN二甲基平面垂直于苯平面，产生电荷分离态，空穴电子基本不重叠，此时S0-S1跃迁振子强度为0。
振子强度为0，对应于跃迁禁阻，荧光淬灭，此时，由于分子振动，打破了禁阻，从而发出荧光。
模拟这种条件下产生的荧光就我所知有如下几种方案：
1，使用gausisian16提供的freq=FCHT，具体做法参考振动分辨的电子光谱的计算，对于本例不是很合适，因为基态和激发态几何结构相差比较大。
2，基于谐振分布产生一批分子，分别计算垂直跃迁(这是本教程中使用的方法）。
3，基于激发态势能面上的绝热动力学产生一批结构，分别计算垂直跃迁。

正文
1.前提条件
电脑中安装了Gaussian09 and Gaussian16 安装方法参考http://bbs.keinsci.com/thread-5493-1-1.html (Gaussian09非必要）
NewtonX2.0 安装方法参考NewtonX2.0安装方法
openbabel

2. 使用Gaussian16优化激发态结构并计算频率
使用gv建立模型，使用如下关键词制作dmabn.gjf文件，用g16计算得到dmabn.log

1. #p b3lyp def2svp td(singlets) opt freq nosym

复制代码

3. 文件准备
依次运行下面的代码准备NewtonX需要的文件

1. obabel -i g09 dmabn.log -o xyz -O dmabn.xyz     #将log文件转换为xyz文件
   
2. xyz2nx < dmabn.xyz                                       #将xyz文件转换为NewtonX需要的geom文件
   
3. cp dmabn.log gaussian.log
   
4. mkdir JOB_AD                                                 #此目录名不可更改
   
5. echo 'def2svp' > JOB_AD/basis                         #计算垂直跃迁的基组在这里设置
   
6. cp dmabn.gjf JOB_AD/gaussian.com                  #将之前的输入文件拷贝过去作为计算垂直跃迁的模板输入文件
   
7. sed -i 's/freq//' ./JOB_AD/gaussian.com
   
8. sed -i 's/opt//' ./JOB_AD/gaussian.com              #删除模板输入文件中的opt和freq关键词

复制代码

4. 产生NewtonX输入文件
运行

1. nxinp

复制代码

命令，NewtonX会采用交互式引导产生输入文件，对新手非常友好，首先选1 ——  1. GENERATE INITIAL CONDITIONS
然后按照他的提示依次进行选择，如果对选项有不明白的地方，可以对照下面产生好的输入文件来完成操作，【回车跳过】表明设定该参数时可以直接回车使用默认值。关键词的详细说明可以参考手册第
运行完成后，当前目录下会产生一个名为initqp_input的文件

&dat

nact= 2                              #获取分子结构的方式，nact=2意思是Wigner distribution for the harmonic oscillator，【回车跳过】

numat = 21                         #体系原子数，【回车跳过】

npoints = 2000                    #一共生成多少个结构，此处需设定，个人建议设大一点，后面可以使用分割的方法只算一部分。

file_geom = geom                #分子坐标，【回车跳过】

iprog = 4                             #读取分子振动信息的来源，选4表明分子振动来自Gaussian的输出文件

file_nmodes = gaussian.log   #含有分子振动信息的Gaussian输出文件名

anh_f = 1                            #非谐矫正系数，如何设置可以参考谈谈谐振频率校正因子

temp= 300                          #设定温度，参考谈谈温度、压力、同位素设定对量子化学计算结果产生的影响

ics_flg = n                           #是否计算光电子能谱，【回车跳过】

chk_e = 1                            #是否计算垂直跃迁，这里默认是0，不可直接跳过，需要改成1

nis= 1                                 #初始态，1为基态，【回车跳过】

nfs= 2                                 #末态，2为S1态（Kasha规则了解一下），NewtonX会统计初态到末态，以及中间所有态的垂直跃迁能，如果算紫外吸收，此处应该根据情况设置20-50左右。

kvert = 1                             #能量限定规则，见下个参数，【回车跳过】

de =100                              #单位eV，如果产生的结构能量超过正负100eV，则认为不合理，剔除，这里设100eV相当于不做限制了，【回车跳过】

prog= 6.5                            #选择计算垂直跃迁的程序，6.5对应于Gaussian TDDFT

iseed = -1                            #设定随机数种子，【回车跳过，注意这里用-1可能会产生过大随机数导致无法运行，可以自己人工随便输一个】

lvprt = 1                              #打印级别，【回车跳过】

5，分割，计算，合并（如果只有一台机器学习使用，可以跳过5.1和5.3的步骤）
5.1 分割
上一步产生的输入文件指定产生2000个构象，这其实是过量的，实际上几百个大概就够了。
实际上一个节点算几百个TD任务也要很长时间，这就需要把任务分割成很多小份提到不同节点上。
输入：

1. split_initcond.pl

复制代码

命令，它会问你需要把任务分解成多少份，这里我们输入100。这样每个任务只需要算2000/100=20个结构的TD了。然后它会问是否需要提交到集群计算，这里需要准备一个作业调度系统的模板，然而各个作业调度系统设定不一样，这里就不单独对这一步说明，如果没有模板，选否即可。
此时当前文件夹下会产生一个名为INITIAL_CONDITIONS 的文件夹，进入后会发现里面有I1 I2 I3 I4 ... I100这样的子文件夹。
5.2 调用gaussian运算
分别进入每个文件夹，把

1. initcond.pl > initcond.log

复制代码

命令提交到作业调度系统上。如果某个文件夹下的任务顺利运行完毕，会自动产生final_output.1.2文件，否则，查看initcond.log 以及DEBUG文件夹中的.error中的信息判断，如果还找不到有用的错误信息，则需要把打印级别lvprt设为2。
5.3 合并
部分文件夹中的任务运行完毕后，就可以运行

1. merge_initcond.pl

复制代码

合并，它会问你需要合并几个文件夹，假设这里我们的I1到I25号文件夹中的任务都运行完毕（500个结构），我们就输入25。
然后会产生一个I_merged 文件夹。其中有个final_output.1.2的文件，其中记录了每个结构的坐标，跃迁能，振子强度等信息，下一步我们将利用这个文件产生光谱。

6. 产生光谱
在I_merged 文件夹下（如果上一步没有执行分割合并的步骤就在当前含有final_output.1.2文件的文件夹下）运行

1. nxinp

复制代码

命令，进入交互式问答环节，首先选5， （5. GENERATE TRAJECTORIES AND SPECTRUM）
然后按照它的提示依次进行选择。
会产生emission-rate.dat的一个文件，这个文件有4列数据，其中第一列是垂直跃迁能eV，第二列是垂直跃迁能nm，第三列是吸收发射的强度，
这里我们将这些数据导入excel表格中，用第二列和第三列作图。得到如下的结果：


7.无法正常运行的调试思路
中间会产生DEBUG和TEMP文件下，DEBUG中可能会找到出错的原因，而TEMP文件夹则是运行中产生的临时文件。根据这两个文件夹下的内容可以诊断出问题的原因。



后记
如何添加自己所需的属性？
默认的finaloutput文件中只给出了能量，振子强度等基本信息，如果用户需要考察激发态的其他信息，则需要自己手动修改源代码，前提是用户对perl语言和Fortran语言有一定的了解。
需要编辑3个文件，首先是：
bin/run_g09_initcond.pl
该脚本调用gaussian，gaussian运行完后从log文件中读取信息，并写入临时文件中，比如oos.2文件记录了第一激发态的振子强度。epot.3文件记录了第2激发态的能量。
如果需要添加自己需要的信息，在该脚本中添加从log文件中抓取其他信息的语句，比如把跃迁偶极抓取并写入tmd.n文件中（n代表第n-1激发态），把激发态偶极写入dip.n文件中。或者是调用Multwfn进行更多强大的分析，比如将空穴电子重叠写入overlap.n文件中。。。
然后是：
source/initcond/writeall.f90
该程序读取上一步中产生的各种临时文件中的信息。然后整合写入readall.out文件中。
修改该程序后，运行make编译后，从EXE_initcond文件夹中将可执行文件拷贝到bin下
最后是：
bin/initcond.pl文件
该程序整合上面两个程序，如果上面两个脚本做了修改，initcond.pl中也需相应的改动，比如把tdm.n拷贝成tdm供writeall程序读取。。。

后记2
如何添加自己所需的属性？（简单版）
输入nxinp后，有一步会问你是否chk_e，如果这里选1，那么就会调用gaussian计算tddft
这里不要选1，而是选0，那么之后运行initcond.pl时它不会调用gaussian，而仅仅产生结构。几秒内就可以完成，产生final_output文件，该文件中包含所有的结构信息。
那么我们拿到这些结构，用这些结构产生Gaussian输入文件，或者别的量化软件的输入文件，就可以为所欲为了。而不需要象后记1中的改原文件那么麻烦。
下面是具体步骤：
首先把finaloutput文件转换成xyz轨迹：

1. echo 0 \n 100 | finaloutput2dynmld.pl   #在有finaloutput文件的地方运行，把finaloutput文件转换成dyn.mld文件，dyn.mld文件就是xyz格式的轨迹文件

复制代码

然后把xyz轨迹，转换成gaussian的 gjf输入文件

1. obabel -i xyz dyn.mld -m -o gjf -O osc #该命令会产生以osc开头的后面跟数字的一系列文件，格式为gjf，不过没有关键词，也没有gjf后缀。

复制代码

如果想转成xyz文件，把-o gjf改成-o xyz就行了，修改文件名：

1. for f in osc[0-9]*; do mv $f `printf osc%05d.gjf ${f#osc}`; done

复制代码

把osc1，ocs2 这样的文件改名为osc00001.gjf， 0sc00002.gjf，这样文件名能够按后面数字大小排序，否则osc11会排在osc1和ocs2之间。


最后，自己写个脚本，把所有的osc开头的gjf 文件中添加上关键词，%nproc，%mem等条目。然后批量提交到服务器上就可以运算了。

后记3
genHOD脚本的介绍：genHOD将后记2中得操作集成到一个脚本中，（此脚本能自己产生合理的随机数，从而避免iseed=-1可能带来的问题）
使用条件：需要路径中安装有obabel以及newtonX中的initcond.pl和finaloutput2dynmld.pl这3个命令，文件的用法如下：

1. genHOD -n 200 -t 300 -s 100 -a 0.95 -o gjf dmabn.log

复制代码

其中dmabn.log是含有td freq计算的任务（当然基态的freq任务也可以） -n产生200个结构，-t设定温度300度，-s 文件编号从100开始 -a设定非谐矫正因子，-o设定输出文件格式。
运行完成之后，会产生一个叫做dmabn的文件夹，其中就有200个gjf文件，文件名是hod0001.gjf 到hod0200.gjf，对应于200个谐振分布的结构（hod的意思是harmonic oscillator distribution）后面需要自己用脚本修改这些gjf的关键词，提交计算后，用脚本批量提取log中的激发态信息进行分析。
更新genHOD：有时候 obabel 转gaussian的log文件会出错。这时候我们可以手动用别的方法把log文件转成相应的mol文件或者xyz文件。保证xyz文件或者mol文件和log文件放在一起。然后运行命令：

1. genHOD -n 200 -t 300 -s 100 -a 0.95 -o gjf dmabn.xyz

复制代码

即可，如果你使用的是xyz文件，那么就不需要调用obabel程序。

lonemen

太赞了！谢谢分享！

alonewolfyang

好东西，果断收藏

ABetaCarw

萌新求问，这个教程是专门针对如下情况的吗？
振子强度为0，对应于跃迁禁阻，荧光淬灭，此时，由于分子振动，打破了禁阻，从而发出荧光。

ggdh

ABetaCarw 发表于 2018-6-8 11:59
萌新求问，这个教程是专门针对如下情况的吗？
振子强度为0，对应于跃迁禁阻，荧光淬灭，此时，由于分子振 ...

荧光发光的情况都可以啊。
只是一般的荧光发光算个TDDFT就行了。没有必要用这招吧

旺旺雪饼

ggdh 发表于 2018-6-9 22:11
荧光发光的情况都可以啊。
只是一般的荧光发光算个TDDFT就行了。没有必要用这招吧

那用这招的目的是什么呢?

yakir

请问这样计算荧光和Gaussian的TDDFT优化激发态计算荧光原理一样吗？

sobereva

yakir 发表于 2018-6-30 23:04
请问这样计算荧光和Gaussian的TDDFT优化激发态计算荧光原理一样吗？

你说的那种简单而常用的做法没有体现出核的振动对光谱的影响

yakir

sobereva 发表于 2018-7-1 05:01
你说的那种简单而常用的做法没有体现出核的振动对光谱的影响

1.老师这个是不是就相当于您写的那个振动分辨的电子光谱？
2.关于荧光的基础知识比如跃迁禁阻、旋轨耦合、荧光猝灭等等老师有推荐的书吗？我想系统的了解一下
谢谢老师

sobereva

yakir 发表于 2018-7-1 14:47
1.老师这个是不是就相当于您写的那个振动分辨的电子光谱？
2.关于荧光的基础知识比如跃迁禁阻、旋轨耦合 ...

1 是
2
Photochemistry of organic compounds-From Concepts to Practice (Klan, Wirz, 2009)
Excited States and Photochemistry of Organic Molecules(1995)

yakir

sobereva 发表于 2018-7-1 18:45
1 是
2
Photochemistry of organic compounds-From Concepts to Practice (Klan, Wirz, 2009)

谢谢老师

Amiswen

“gaussian运行完后从log文件中读取信息,”出现了这种错怎么解，“Cannot open gaussian.log! at /share/apps/chemsoft/NX-2-B19/bin/lib/colib_perl.pm line 363, <MEM> line 1.”

ggdh

Amiswen 发表于 2018-7-17 20:28
“gaussian运行完后从log文件中读取信息,”出现了这种错怎么解，“Cannot open gaussian.log! at /share/ap ...

看错误提示，估计是你没有把算频率的文件重命名为gaussian.log把

Amiswen

ggdh 发表于 2018-7-18 21:31
看错误提示，估计是你没有把算频率的文件重命名为gaussian.log把

非常感谢！
另外，我发现用2.0版本时，emission-rate.dat里的第三列不是吸收发射强度，而是radiative decay rate，那么横坐标是能量，纵坐标对应radiative decay rate的图怎么理解呢？

ggdh

Amiswen 发表于 2018-11-10 12:44
非常感谢！
另外，我发现用2.0版本时，emission-rate.dat里的第三列不是吸收发射强度，而是radiative de ...

速率就是强度啊，发射速率越快，实验上观察到的强度就越强。