计算化学公社
标题: 求助：如何利用内积的方法确定振动模式 [打印本页]
作者Author: wangy0822 时间: 2022-3-10 14:34
标题: 求助：如何利用内积的方法确定振动模式
本帖最后由 wangy0822 于 2022-3-10 14:38 编辑
不好意思。上次问题没有描述清楚，重新请教一下大家，举个例子说明应该更清楚点，下面的两个图是我建的两个不同的模型（其实差别不大）的Gaussian计算频率的结果，虽然比较明确的知道这两个模型的第一个振动模式是一样的，但是就是想通过振动模式的内积来确定一下，好像是可以通过振动模式对应矩阵的内积来做，然后看数值是不是接近1来判断振动模式是不是一样。
请问下大家有没有了解这个怎么去做，而且矩阵内积的话两个矩阵不应该是m*n，n*m这种形状吗，我在图中红线框出的不能做内积吧，如果不是这样做的话应该怎么去做？
作者Author: wzkchem5 时间: 2022-3-10 15:47
这个内积指的是把这两个3*N的矩阵分别展开成1*3N的向量，再做内积。等价于求这两个矩阵的element-wise product再求和。
作者Author: wangy0822 时间: 2022-3-10 16:01
谢谢您的解答，那么如何将这个3*N的矩阵展开成1*3N的向量
作者Author: wzkchem5 时间: 2022-3-10 17:22
把第i行第j列的矩阵元放到向量的第3(i-1)+j个位置
作者Author: wangy0822 时间: 2022-3-11 09:49
本帖最后由 wangy0822 于 2022-3-11 09:51 编辑
感谢您的解答，您的意思应该是下面那张图的意思吧，还有几个问题请教您一下
（1），按照您说的方法做内积的话对应元素相乘再相加之后数值是接近于多少才是振动模式一样呢？
（2），那如果模型原子数不相同的话那么元素的个数也不一样，那是不是就比较不了了呢？
（3），VASP里面的数值和Gaussian里面的数值差的挺多的（下图所示），是VASP里面需要除以某个数吗
作者Author: wzkchem5 时间: 2022-3-11 16:36
对
（1）只能说，如果第一个分子的第m个振动模和第二个分子的第n个振动模的内积的绝对值大于1/sqrt(2)，就可以说第一个分子的第m个振动模和第二个分子的第n个振动模的相似度，比和第二个分子的其他任一个振动模的相似度都高。至于如何判断一样，显然取决于研究目的和所需精度，不能一概而论，绝对值越接近1越一样，但是一样和不一样的分界线画在哪里没有普适标准。
（2）对
（3）你圈出来vasp里的是原子坐标吧？不是简正模吧？
作者Author: wangy0822 时间: 2022-3-11 19:38
本帖最后由 wangy0822 于 2022-3-11 19:58 编辑
感谢您的解答，上面的那个X,Y,Z下面的应该是坐标信息，每个振动模式都一样，dx，dy，dz下面是简正振动模吗？我看VASP里面的数值比Gaussian里面小不少，我算了下不同泛函算的同一模型在同一个振动模式（大致判断是同一个振动模式，可视化是一样的）乘了一下绝对值才0.1多，感觉不太对，还是说VASP的输出文件里没有简正振动模.。。我上传了用VASP两个泛函算振动的OUTCAR，您可以帮忙看下吗？感谢感谢，我算的是在237个振动模式的内积数值绝对值才0.1多。
作者Author: wangy0822 时间: 2022-3-12 10:25
还有个问题，Gaussian输出文件里面是简正模吗。。感谢解答
作者Author: wzkchem5 时间: 2022-3-12 16:04
高斯输出的那个是的。我没用VASP算过简正模，不知道是不是从这里读的
作者Author: wangy0822 时间: 2022-3-12 16:26
好的，谢谢您的解答，从数值上来看VASP里面的应该不是简正模，数值要比Gaussian里面小一个数量级；这个简正模是和相对分子质量有关吗？还是和什么有关？
作者Author: wzkchem5 时间: 2022-3-12 17:46
简正模是mass-weighted Hessian矩阵对角化得到的，mass-weighted Hessian矩阵是由各个原子的质量和Hessian矩阵计算得到的（也就是说不由分子总质量唯一决定，而是和每个原子各自的质量也有关），Hessian矩阵是由分子各个自由度的Hooke劲度系数以及自由度之间的耦合强度决定的
作者Author: wangy0822 时间: 2022-3-12 20:19
感谢楼上这位大佬的一直的解答，剩下最后一个问题问一下其他人。。问题：这个VASP算的振动里怎么去得到简正振动模式？有没有什么公式
作者Author: granvia 时间: 2022-3-12 20:35
先把各3*N维向量归一化，然后再求内积的平方
作者Author: wangy0822 时间: 2022-3-14 09:17
感谢您的解答，还有几个问题
（1）您的意思是把dx，dy，dz这三列分别归一化吗？
（2）归一化是将数据限制在[0，1]之内，最小的为0，最大的为1，是这个意思吗？
（3）您说的再求内积的平方是什么意思，求内积之后再平方？还是什么意思，可以再解释清除一下吗？
作者Author: wangy0822 时间: 2022-3-14 14:36
我试了一下，VASP结果与VASP结果对比，先把每一列归一化（数据限制在[0，1]之间），然后内积之后平方出来的数大于1了，不太对吧；然后VASP结果与Gaussian结果对比，就是把VASP结果归一化，然后平方之后与Gaussian结果做内积，这样对吗？还请大家不吝赐教
作者Author: granvia 时间: 2022-3-16 00:02
wangy0822 发表于 2022-3-14 14:36
我试了一下，VASP结果与VASP结果对比，先把每一列归一化（数据限制在[0，1]之间），然后内积之后平方出来 ...
向量的归一化就是向量除以它的模长（等于它自己与自己内积的平方根）。两个归一化向量的内积平方一定小于一
作者Author: wangy0822 时间: 2022-3-16 09:18
感谢您的解答，我昨天请教了一下老师。说VASP里面算的得到的结果是简正模，Gaussian里面计算频率时要加入freq=hpmodes这个参数可以得到简正模
作者Author: wangy0822 时间: 2022-3-16 11:09
到底VASP里面输出的是简正模还是Gaussian里面输出的是简正模
。。。哪个输出的数据需要处理，哪个数据可以直接使用
作者Author: granvia 时间: 2022-3-16 18:15
wangy0822 发表于 2022-3-16 11:09
到底VASP里面输出的是简正模还是Gaussian里面输出的是简正模。。。哪个输出的数据需要处理，哪个数据 ...
很简单，找到简单分子HF, H2O来测试一下就知道了
作者Author: wangy0822 时间: 2022-3-16 19:37
感谢您的解答，HF的测试结果，怎么看是不是简正模？
作者Author: wangy0822 时间: 2022-3-16 19:42
我现在是单纯的VASP结果可以做内积验证，单纯的Gaussian结果也可以做内积验证，就是两者的数据不是一种类型，结合不起来，不能交叉验证。
作者Author: granvia 时间: 2022-3-16 21:18
本帖最后由 granvia 于 2022-3-16 21:58 编辑
Gaussian默认输出的振动位移（对应于质量权重的Hessian矩阵的特征向量）精度是小数点后2位。使用FREQ=HPModes之后，可输入小数点后5位的振动位移向量，你可以验证它们就是归一化的。比如下面是H2O分子的结果（B3LYP/6-31G*）：
- Harmonic frequencies (cm**-1), IR intensities (KM/Mole), Raman scattering
- activities (A**4/AMU), depolarization ratios for plane and unpolarized
- incident light, reduced masses (AMU), force constants (mDyne/A),
- and normal coordinates:
- 1 2 3
- A1 A1 B2
- Frequencies --- 1713.2088 3728.7514 3850.4642
- Reduced masses --- 1.0825 1.0454 1.0810
- Force constants --- 1.8719 8.5633 9.4426
- IR Intensities --- 75.7520 1.6945 19.3415
- Coord Atom Element:
- 1 1 8 -0.00000 0.00000 -0.00000
- 2 1 8 0.00000 -0.00000 -0.06986
- 3 1 8 -0.07058 0.05004 -0.00000
- 1 2 1 0.00000 0.00000 -0.00000
- 2 2 1 0.42879 0.58401 0.55440
- 3 2 1 0.56005 -0.39709 -0.43612
- 1 3 1 0.00000 -0.00000 0.00000
- 2 3 1 -0.42879 -0.58401 0.55440
- 3 3 1 0.56005 -0.39709 0.43612
最后3列就是3个9x1的振动位移向量（9=3*N，N为原子数），然后就可以求出对应的3个质量权重的Hessian矩阵的特征列向量，即将上述各振动位移向量中的各元素乘以相应原子的质量的平方根，然后再归一化（因为Gaussian打印振动位移向量时做了归一化，所以你变换回Hessian特征列向量时就需要自己归一化一次）。你可以验证，上述求出的Hessian特征列向量之间是正交归一的，也就是简正模式的各特征向量。VASP我不太清楚，可以做类似尝试
作者Author: wangy0822 时间: 2022-3-16 22:27
本帖最后由 wangy0822 于 2022-3-16 22:29 编辑
感谢您的解答你的这个输出文件好像和我的不太一样，我的每个频率下面都有x，y，z
还有，我把VASP和Gaussian里面都化成正交归一就可以做内积比较了吧
作者Author: granvia 时间: 2022-3-16 23:57
wangy0822 发表于 2022-3-16 22:27
感谢您的解答你的这个输出文件好像和我的不太一样，我的每个频率下面都有x，y，z
还有，我把VASP和Gauss ...
输出文件中再往上找找，就能找到
作者Author: wangy0822 时间: 2022-3-17 14:19
感谢大家的解答
好像有个问题，就是我Gaussian算的团簇模型。VASP算的晶体结构模型，虽然都是为了模拟某种分子晶体的振动。但是毕竟结构不一样，所以是不是虽然振动频率接近，但是做出来的内积都非常小（大概0.01），就是从内积上看基本不是一个振动模式的那种。
我之前就一直认为频率接近它们的内积就要接近与1，所以一直认为数据处理的有问题。而且我也可视化看了振动模式，晶体结构的振动模式完全和团簇的不一样，所以这是不是就能说得通为什么内积会这么小了。
作者Author: granvia 时间: 2022-3-17 18:25
wangy0822 发表于 2022-3-17 14:19
感谢大家的解答
好像有个问题，就是我Gaussian算的团簇模型。VASP算的晶体结构模型，虽然都是为了模拟某种 ...
你得保证你原子的排列顺序和取向都尽量一致吧？ 因为简正模式的特征向量是3*N维的，你得保证各个维度是一致的
作者Author: wangy0822 时间: 2022-3-17 19:23
您看一下，这是我算的团簇模型和晶体结构模型在某一频率的振动模式，明显不一样吧，其实也不能说不一样，反正做内积肯定不是接近1的
作者Author: wangy0822 时间: 2022-3-17 19:23
本帖最后由 wangy0822 于 2022-3-17 19:25 编辑
您看一下，这是我算的团簇模型和晶体结构模型在某一频率的振动模式，明显不一样吧，其实也不能说不一样，反正做内积肯定不是接近1的
作者Author: granvia 时间: 2022-3-17 19:57
wangy0822 发表于 2022-3-17 19:23
您看一下，这是我算的团簇模型和晶体结构模型在某一频率的振动模式，明显不一样吧，其实也不能说不一样， ...
这谁看的出来哈 您究竟想解决什么问题？如果不能准确描述清楚，我也没法帮您出主意
作者Author: wangy0822 时间: 2022-3-17 21:01
本帖最后由 wangy0822 于 2022-3-17 21:10 编辑
不好意思，可能没解释清楚。
我之前呢是想采用振动模内积的方法对比Gaussian计算的团簇模型的振动模式与VASP计算的晶体结构的振动模式是否一样。
就像上面两张图，我就想用振动模内积看看它们是不是接近于1来判断是不是同一个振动模式，然后尝试几天，发现数据怎么做内积都没办法接近1
今天好像发现了错误：不是振动模式一不一样的问题，就是本身这两个结构就差别很大，就像您之前说的，原子的排列和取向都不同（上面两张图），本身就不可能得到内积接近于1. 然后振动模式就没有什么可比性。 是不是这样？
作者Author: granvia 时间: 2022-3-17 21:19
wangy0822 发表于 2022-3-17 21:01
不好意思，可能没解释清楚。
我之前呢是想采用振动模内积的方法对比Gaussian计算的团簇模型的振动模式与 ...
即使同一分子的同一振动模式，原子排序不一致以及x,y,z轴取向不一致都会使其内积远小于1。 但我不认为两个相似模型的振动模式没可比性，应该存在一种elegant的方法来描述这种相似性的，可能需要好好想想
作者Author: wangy0822 时间: 2022-3-17 21:26
感谢您的解答，问题就是我的团簇模型和晶体结构模型差距还是很大的，您可以看我前面的两个振动模式展示的模型，我认为不能达到相似这种程度。我感觉这种elegant方法还无从得知。
作者Author: smutao 时间: 2022-3-18 07:41
对离域性比较强的振动模式，直接对比团簇和晶体，很难做到
对一些比较定域的振动（比如C=O振动），直接对比可以做
假如想对比团簇体系和晶体体系中同一种化学键的强度（力常数或者振动频率），可以使用LModeA-nano做到，这个程序本身就是从振动的角度来计算化学键强度的，也许满足你需要
LModeA-nano 的长处就是可以直接把晶体体系和分子体系的化学键强度进行对比 （前提是两者在同一个级别下计算）
作者Author: sobereva 时间: 2022-3-18 07:51
对于算晶体体系的正则振动模式，团簇模型既没法代替晶体模型，也和周期性算的没有可比性
参看Phys. Chem. Chem. Phys., 2021, 23, 20038里的讨论。而且光是从你的截图里就已经体现这一点了，晶体计算给出的是集体振动模式。就算簇模型和周期性计算给你的频率表面上看可能有的相差不大，但振动模式却也可能相差甚巨。
作者Author: wangy0822 时间: 2022-3-18 08:48
感谢各位老师的解答
作者Author: smutao 时间: 2022-3-18 09:36
本帖最后由 smutao 于 2022-3-18 10:16 编辑
假如你的团簇和晶体包含同一种分子，想比较它在两个体系里的振动
有一种方案是将其视为片段 然后考察这个片段的振动模式
这个方案需要使用GSVA方法 (见10.1021/acs.jctc.7b01171 和 10.1007/s00214-021-02727-y) 相关帖文：http://bbs.keinsci.com/thread-23222-1-1.html
团簇体系很好解决 现有的程序就能做了
晶体方面 需要提取其3M-6个振动对应的Hessian （输入的信息为整个体系的Hessian，可以用UniMoVib做GSVA）
你想继续深入的话可以告诉我
作者Author: wangy0822 时间: 2022-3-18 09:51
好的好的，谢谢老师，我先了解一下您说的这个
作者Author: wangy0822 时间: 2022-3-23 20:56
今天请教了下老师，关于这个两种不同的模型如何去做内积比较振动模式一不一样，此帖之前已经讨论过一种团簇模型一种晶体结构模型，原子排列和取向不一样没办法去做内积，今天老师说可以旋转原子去让两种模型排列，取向一致然后再去做内积比较，请教一下大家，这种方法可不可行，团簇里的原子能旋转到晶体结构里面的一样或者近似？感觉有点悬呢。。
作者Author: zjxitcc 时间: 2022-3-23 23:16
你可以先找个只含几个原子的单胞试试看，保证孤立体系和单胞里的目标原子有最大重叠（甚至完全重合），然后做计算，到底差的大不大，一看就知道
作者Author: wangy0822 时间: 2022-3-28 21:02
感谢您的解答，大家有没有知道这种方法，采用一个3*3的旋转矩阵，然后让两个模型的rmsd最小，大家有没有了解过这个3*3的矩阵应该怎么去求？模型见图
作者Author: zjxitcc 时间: 2022-3-28 21:32
本帖最后由 zjxitcc 于 2022-3-28 21:36 编辑
可以用VMD，指定部分原子进行旋转和比较，看Sob老师这篇介绍《在VMD中计算RMSD衡量两个结构间的差异以及叠合两个结构》http://sobereva.com/290
也可以自己手写，我自己就写过，可以参考https://gitlab.com/jxzou/rmsd
注：要比较的那部分的朝向可以不一样，但要保证1号文件的每个原子与2号文件的每个原子一一对应，碳对碳，氢对氢，甲基对甲基，这边1号碳原子对应那边的1号碳原子，否则无论是VMD还是我写的程序都没法给出正确结果。
作者Author: wangy0822 时间: 2022-3-28 21:37
非常感谢，我先去试试
作者Author: zjxitcc 时间: 2022-3-28 21:44
你的体系太大了，建议先弄个简单情况测试一下，发现可行之后再考虑目标体系，不然折腾半天发现不可行就费时间了。
作者Author: wangy0822 时间: 2022-3-28 21:57
您好，我看了一下Sob老师的那篇博文，还有您自己写的那个，Sob老师博文里是讲了如何去计算两个模型的RMSD，但其实我的目标是找到一个3*3的旋转矩阵，使得两个模型的RMSD最小，本意不是去计算RMSD，您自己写的里面有个Note2: only the Gaussian gjf format is supported.但是我的是一个Gaussian的gjf，一个是VASP的POSCAR，所以还是无法计算。
作者Author: zjxitcc 时间: 2022-3-28 22:41
这还不简单，坐标复制出来自己写个gjf就得了。gjf也是文本格式，就是坐标加几行内容而已。旋转矩阵的话，你可以从我写的代码里加点打印语句，提取。
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |