计算化学公社
标题: 基于fch中的Hessian矩阵计算振动频率的简单程序Hess2freq [打印本页]
作者Author: sobereva 时间: 2016-5-27 07:30
标题: 基于fch中的Hessian矩阵计算振动频率的简单程序Hess2freq
基于fch中的Hessian矩阵计算振动频率的简单程序Hess2freq
Hess2freq: A simple program to calculate vibrational frequencies based on Hessian matrix in fch
文/Sobereva @北京科音 2016-May-27
写了个教学目的的程序Hess2freq，可以读取fch中的Hessian矩阵来计算谐振频率和正则坐标。给了非常易读的源代码，而且文档里把用到的公式和原理都详细介绍了，对想了解振动分析原理和细节的人应该会很有益。
下载地址： http://sobereva.com/soft/Hess2freq.rar
输出例子：
Hess2freq: Load Hessian from .fch file and then compute harmonic frequencies
Programmed by Sobereva (Sobereva@sina.com)
Release date: 2016-May-27
The number of atoms: 4
The number of vibrational modes: 6
Atomic masses:
12.000000 1.007825 1.007825 15.994915
************ Force constant matrix (i.e mass-weighted Hessian) ************
1 2 3 4 5
1 0.918452D-01 0.427913D-01 0.170607D-05 0.235189D-01 -0.105091D-02
2 0.427913D-01 0.504781D-02 -0.520948D-05 -0.159298D-01 -0.181706D-01
3 0.170607D-05 -0.520948D-05 0.534475D-02 -0.499779D-05 -0.157484D-05
4 0.235189D-01 -0.159298D-01 -0.499779D-05 0.332464D-02 -0.708577D-02
5 -0.105091D-02 -0.181706D-01 -0.157484D-05 -0.708577D-02 0.510611D-01
6 -0.192394D-04 0.167567D-04 0.216766D-02 0.189037D-04 -0.886049D-06
7 -0.852345D-01 -0.332194D-01 0.279538D-04 -0.285394D-01 -0.259396D-02
8 -0.478064D-01 0.343795D-01 0.187013D-04 0.414582D-01 -0.178675D-01
9 0.506723D-04 0.711252D-05 -0.114972D-01 -0.112089D-05 0.343570D-05
10 -0.640613D-01 -0.247270D-01 -0.724008D-05 -0.140419D-01 0.334003D-02
11 -0.248003D-01 -0.844094D-02 0.213252D-06 0.516973D-02 0.740654D-02
12 -0.936792D-05 -0.147931D-05 -0.228756D-02 -0.134884D-06 0.724064D-06
6 7 8 9 10
1 -0.192394D-04 -0.852345D-01 -0.478064D-01 0.506723D-04 -0.640613D-01
2 0.167567D-04 -0.332194D-01 0.343795D-01 0.711252D-05 -0.247270D-01
3 0.216766D-02 0.279538D-04 0.187013D-04 -0.114972D-01 -0.724008D-05
4 0.189037D-04 -0.285394D-01 0.414582D-01 -0.112089D-05 -0.140419D-01
5 -0.886049D-06 -0.259396D-02 -0.178675D-01 0.343570D-05 0.334003D-02
6 0.865363D-03 0.273490D-04 -0.423721D-04 -0.466045D-02 0.505424D-05
7 0.273490D-04 0.336434D+00 0.848079D-01 -0.147836D-03 -0.345946D-02
8 -0.423721D-04 0.848079D-01 -0.715275D-01 -0.133212D-04 0.971329D-02
9 -0.466045D-02 -0.147836D-03 -0.133212D-04 0.247107D-01 -0.649994D-05
10 0.505424D-05 -0.345946D-02 0.971329D-02 -0.649994D-05 0.598807D-01
11 -0.365557D-05 0.813642D-02 -0.733868D-02 -0.367918D-05 0.181410D-01
12 -0.924923D-03 0.603159D-05 -0.221842D-05 0.492551D-02 0.663398D-05
11 12
1 -0.248003D-01 -0.936792D-05
2 -0.844094D-02 -0.147931D-05
3 0.213252D-06 -0.228756D-02
4 0.516973D-02 -0.134884D-06
5 0.740654D-02 0.724064D-06
6 -0.365557D-05 -0.924923D-03
7 0.813642D-02 0.603159D-05
8 -0.733868D-02 -0.221842D-05
9 -0.367918D-05 0.492551D-02
10 0.181410D-01 0.663398D-05
11 0.729420D-02 0.165643D-05
12 0.165643D-05 0.977187D-03
Diagonalization passed
************ Normal coordinates (columns) ************
1 2 3 4 5
1 0.0666 0.0130 0.0000 0.0513 0.4179
2 -0.1089 0.1194 -0.0000 -0.0833 0.1751
3 -0.0000 0.0001 0.1312 -0.0000 0.0001
4 -0.3936 -0.9639 0.0004 -0.1564 0.1945
5 0.0303 -0.0534 0.0001 0.9406 -0.2801
6 0.0004 0.0001 0.1834 0.0000 -0.0000
7 -0.1713 0.0346 -0.0005 0.0914 0.7221
8 0.8931 -0.2116 -0.0001 -0.2670 0.0714
9 -0.0001 -0.0004 -0.9730 0.0000 -0.0002
10 -0.0144 0.0488 -0.0000 -0.0344 -0.3713
11 0.0235 -0.0729 0.0000 0.0200 -0.1182
12 0.0000 -0.0000 -0.0487 -0.0000 -0.0000
6
1 -0.1015
2 -0.0311
3 0.0000
4 -0.0689
5 -0.0080
6 0.0001
7 0.9710
8 0.2015
9 -0.0005
10 0.0193
11 0.0111
12 0.0000
The frequencies (cm-1) corresponding to overall translation and rotation:
-0.00997 0.01937 -0.05532 -12.87321 -15.76915 -19.42877
Harmonic vibrational frequencies:
Mode 1: -0.56590E+14 Hz -1887.64948 cm-1
Mode 2: 0.23453E+14 Hz 782.31259 cm-1
Mode 3: 0.27533E+14 Hz 918.38630 cm-1
Mode 4: 0.38543E+14 Hz 1285.65606 cm-1
Mode 5: 0.58055E+14 Hz 1936.50274 cm-1
Mode 6: 0.96313E+14 Hz 3212.64652 cm-1
作者Author: 我本是个娃娃 时间: 2016-5-27 08:49
沙发
作者Author: 978142355 时间: 2016-5-27 10:16
大赞sob老师。问一个问题，Lapack库在windows上是如何配置到VS2010中的？看了许多教程，都未成功，因为自己一直编写小程序从来没用什么数学库，借此机会问问老师这个该如何配置？我在网站上下载了windows版本的lapack库，打开文件夹里面好多的东西，不知道哪些应该放在VS中，求老师指点。（PS：没有看到*.lib文件，只找到了include文件夹，所以不知该如何是好）
作者Author: jliang 时间: 2016-5-27 10:18
Sob老师的物理功底真深！相比之下让学物理的人都感到惭愧！
作者Author: sobereva 时间: 2016-5-27 16:12
blas库在CVF/ivf下的使用：
新建一个工程，类型选静态库。把blas库所有.f都添加进source file里，然后build，就得到了.lib文件。
lapack在CVF/ivf下的使用
把lapack压缩包内的SRC内的.f都加入工程，把INSTALL里的辅助子程序slamch.f和dlamch.f加入子程序。把BLAS库文件加入工程。build即可。
blas在linux下编译：
解压后，把make.inc里面编译器改成自己的。然后make即可。
lapack在linux下编译：
解压后，把make.inc.example改名为make.inc，把FORTRAN、LOADER都改为ifort，把OPTS改成-O3，把BLASLIB改成blas.a的路径。然后make即可。
作者Author: 978142355 时间: 2016-5-27 20:11
大赞sob老师，谢谢，十分感谢，总能在这里找到十分正确又细致的答案，向sob学习，向sob看齐。
作者Author: 让你变成回忆 时间: 2017-9-11 19:39
sob老师，您好！
目前在学习Fortran编程，正在跟着您写的Hess2freq程序练习。
练习过程中发现您的程序中 loclabel子程序部分自己不怎么理得清关系，您可以简单说一下是怎么在某个文件中找到我要的label的呢？
谢谢sob老师。
作者Author: sobereva 时间: 2017-9-11 20:29
loclabel就是定位到含有某一串字符的那一行的子程序
从文本文件中读取信息的时候，得先定位到合适的位置，才能开始读取，loclabel就是用来干这个的。每一部分数据输出之前肯定会有一些提示字符能用来定位。
作者Author: 让你变成回忆 时间: 2017-9-11 20:51
subroutine loclabel(fileid,label,ifound,irewind)
integer fileid,ierror
integer,optional :: ifound,irewind
character*200 c200
CHARACTER(LEN=*) label
if ((.not.present(irewind)).or.(present(irewind).and.irewind==1)) rewind(fileid)
do while(.true.)
read(fileid,"(a)",iostat=ierror) c200
if (index(c200,label)/=0) then
backspace(fileid)
if (present(ifound)) ifound=1 !Found result
return
end if
if (ierror/=0) exit
end do
if (present(ifound)) ifound=0
end subroutine
这是sob老师您写的代码，我有几个不明白的地方，恳请您提示下：
1、这里面 ifound、irewind意义是什么，为什么不能直接通过index找到某个特定字符就行。
2、另外我发现，由于具体的字符出现的位置不同，您在主程序中写的也不同。比如 需要原子数的时候，您的程序中并没有read一个空行，我看fchk文件这是由于原子数就出现在查找字符串（原子数）的同一行，但导入质量矩阵等其他的时候，有一个read的空行，是由于这在fchk文件中换行了的原因吗？
因为最近一直在跟着sob老师您写的程序练习编程（如Hess2freq，optDFTW等），这部分内容不是很了解，所以冒昧请教您。 感谢~
作者Author: sobereva 时间: 2017-9-11 21:00
ifound判断有没有找到相应字符串，irewind用来决定是否做rewind。loclabel不是为了hess2freq专门写的，而是我写的各种程序中都会用到的一个子程序，比较通用。
loclabel将读取位置定位到含有某个字符串的行，如果要读的数据就在那一行，就不需要先read一下跳到下一行。
作者Author: 让你变成回忆 时间: 2017-9-11 21:02
对，我就是看您的很多程序中都有这样一样loclabel的子程序，所以想了解清楚点一下具体的过程。
大致懂了。
感谢sob老师的回答。谢谢~
作者Author: 让你变成回忆 时间: 2018-4-20 20:53
sob老师您好，再次请教您编程问题：
在您的程序中用到了showmatgau这个子程序。为什么不单独列成子程序来使用呢？而是在主程序中加入contains 后面接该子程序。
我尝试把这个子程序单独作为模板来使用，因为我觉得这样输出矩阵超级方便。但是我把这部分内容单独写为一个子程序，然后再去调用，就会出现问题，而向您那样放入主程序中，并在之前加入contains，就没有任何问题。
翻阅了一些书籍，没有找到相关的内容。恳求您再次解答。谢谢sob老师。
作者Author: sobereva 时间: 2018-4-21 07:57
只是为了当前这个小程序方便才这样写。
showmatgau本是Multiwfn程序里的子程序，是独立的。
作者Author: tiandikuoyuan 时间: 2021-7-24 21:37
老师，请问软件输出的Normal coordinates与Gaussian 输出文件中的normal coordinates是否是完全一致的呀，我发现部分数据正负相反，有些数据一致；比如说示例文件ethanol.out中 19 20 21
A A A
Frequencies -- 3276.7816 3289.0510 4114.6169
Red. masses -- 1.1023 1.1047 1.0669
Frc consts -- 6.9734 7.0412 10.6423
IR Inten -- 48.8354 53.7301 41.2114
Raman Activ -- 63.3543 47.1572 95.0383
Depolar (P) -- 0.7457 0.7500 0.3130
Depolar (U) -- 0.8544 0.8571 0.4767
Atom AN X Y Z X Y Z X Y Z
1 6 -0.06 -0.07 0.00 0.00 0.00 0.09 0.00 0.00 0.00
2 1 0.01 0.26 -0.39 0.03 0.41 -0.56 0.00 0.00 0.00
3 1 0.64 0.36 0.00 0.00 0.00 0.01 0.00 0.00 0.00
4 1 0.01 0.26 0.39 -0.03 -0.41 -0.56 0.00 0.00 0.00
5 6 0.00 -0.01 0.00 0.00 0.00 0.02 0.00 0.00 0.00
6 1 0.00 0.03 -0.04 0.00 0.06 -0.09 0.00 0.00 0.00
7 1 0.00 0.03 0.04 0.00 -0.06 -0.09 0.00 0.00 0.00
8 8 0.00 0.00 0.00 0.00 0.00 0.00 0.05 -0.04 0.00
9 1 0.00 0.00 0.00 0.00 0.00 0.00 -0.79 0.61 0.00
Hess2freq软件输出的数据如下：
19 20 21
1 -0.0553 -0.0000 0.0000
2 -0.0742 -0.0000 -0.0003
3 0.0000 -0.0926 -0.0000
4 0.0100 -0.0300 0.0003
5 0.2638 -0.4114 -0.0002
6 -0.3889 0.5599 0.0002
7 0.6428 0.0000 0.0012
8 0.3618 0.0000 0.0011
9 0.0000 -0.0144 -0.0000
10 0.0100 0.0300 0.0003
11 0.2638 0.4114 -0.0002
12 0.3889 0.5599 -0.0002
13 -0.0011 -0.0000 0.0005
14 -0.0061 -0.0000 0.0019
15 0.0000 -0.0156 0.0000
16 0.0012 -0.0002 0.0026
17 0.0307 -0.0635 0.0004
18 -0.0420 0.0883 0.0007
19 0.0012 0.0002 0.0026
20 0.0307 0.0635 0.0004
21 0.0420 0.0883 -0.0007
22 0.0004 0.0000 -0.0506
23 0.0003 -0.0000 0.0371
24 -0.0000 0.0003 -0.0000
25 -0.0001 -0.0000 0.7910
26 0.0013 0.0000 -0.6086
27 0.0000 0.0006 0.0000
我想通过这个软件将fchk文件转换成orca可识别的hess文件，发现三个文件中Normal coordinates都不完全一致，请问下是哪里存在问题吗？
Gaussian 水的数据：
1 2 3
A1 A1 B2
Frequencies -- 1616.6354 3787.5994 3892.8460
Red. masses -- 1.0838 1.0441 1.0826
Frc consts -- 1.6689 8.8250 9.6661
IR Inten -- 80.7528 4.4711 44.1094
Atom AN X Y Z X Y Z X Y Z
1 8 -0.00 -0.00 -0.07 0.00 0.00 0.05 -0.00 -0.07 0.00
2 1 0.00 0.42 0.57 0.00 0.59 -0.39 0.00 0.56 -0.43
3 1 0.00 -0.42 0.57 -0.00 -0.59 -0.39 0.00 0.56 0.43
Hess2freq软件输出的数据如下：
************ Normal coordinates (columns) ************
1 2 3
1 -0.0000 -0.0000 0.0000
2 0.0000 -0.0000 -0.0706
3 -0.0712 -0.0492 0.0000
4 -0.0000 -0.0000 -0.0000
5 0.4220 -0.5886 0.5605
6 0.5651 0.3903 -0.4282
7 0.0000 0.0000 -0.0000
8 -0.4220 0.5886 0.5605
9 0.5651 0.3903 0.4282
orca 计算得到的out文件：
------------
NORMAL MODES
------------
These modes are the Cartesian displacements weighted by the diagonal matrix
M(i,i)=1/sqrt(m[i]) where m[i] is the mass of the displaced atom
Thus, these vectors are normalized but *not* orthogonal
6 7 8
0 0.000000 0.000000 0.000000
1 0.000000 0.000000 -0.070639
2 0.071235 0.049153 0.000000
3 0.000000 0.000000 0.000000
4 -0.421752 0.588752 0.560593
5 -0.565321 -0.390079 -0.428068
6 0.000000 0.000000 0.000000
7 0.421752 -0.588754 0.560591
8 -0.565321 -0.390080 0.428067
作者Author: wzkchem5 时间: 2021-7-24 21:50
normal coordinates只有方向是有物理意义的，把整个normal coordinate矢量取负号，不改变任何物理可观测量。即使都是高斯做的计算，你把结构稍微扰动一下再做opt+freq，normal coordinate符号也会有变化。
作者Author: tiandikuoyuan 时间: 2021-7-25 10:23
明白了，谢谢
作者Author: Kalinite 时间: 2022-6-8 16:57
sob老师您好！
说明文档中提到：“Gaussian输出最终的振动频率前，会先从3N个频率中自动判断出整体平动和转动模式来，然后将它们对力常数矩阵的贡献投影掉使得振动频率更精确。”这里是如何判断整体平动转动模式的呢，是否也是和这个程序中一样去掉最小的六个。另外这里的投影可以给一下具体的数学表达吗？最近在做一些工作，想把一些笛卡尔坐标系下的矢量转换到内坐标下，对如何删除平动转动自由度感到很是困惑。
作者Author: wzkchem5 时间: 2022-6-8 18:08
我不知道sob老师的做法是什么，但是一般量化程序的做法并不是对角化以后去掉最小的6个频率，而是先推导出平动和转动模式的解析形式，再把这些模式从Hessian里投影掉，然后对角化剩下的3N-6维的Hessian。这样的好处是假如有振动频率极其小的振动模式，不至于误把这个振动扔掉。
作者Author: Kalinite 时间: 2022-6-8 20:49
请问您有关于“推导出平动和转动模式的解析形式”这方面的资料吗
作者Author: wzkchem5 时间: 2022-6-8 20:56
没有，这个就是经典力学的内容，找一本力学教材看看刚体平动和转动的简正模怎么算就知道了
作者Author: sobereva 时间: 2022-6-9 07:38
Gaussian有个官方文档Vibrational Analysis in Gaussian，里面有细节
作者Author: WaterMirror 时间: 2022-6-9 08:24
Gaussian官方文档Vibrational Analysis in Gaussian
作者Author: Kalinite 时间: 2022-6-9 13:27
感谢大家的帮助！
作者Author: 王二葛 时间: 2022-6-10 11:41
推荐 ajz34 的讲解和代码
- 频率分析 (1)：从 fchk 文件得到分子频率与简正模式
- https://ajz34.readthedocs.io/zh_CN/latest/QC_Notes/Freq_Series/freq_1.html
另推荐他写的 python 实现双杂化泛函梯度教程，其中的量化入门部分很友好
- xDH 在 Python 下实现的简易教程
- https://py-xdh.readthedocs.io/
作者Author: 风起~ 时间: 2022-6-10 11:49
谢谢
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |