计算同位素如何设置

asl1994

想用VASP做同位素交换，体系中两种同位素发生互换，但是不知道应该如何设置，看到有人直接修改POTCAR中的POMASS，可是这样体系中仍然只有一种同位素。请问怎么做到两种同位素共存呢，比如H和D？谢谢！

高阁

AutMaple 发表于 2025-4-26 10:03
这个应该是u=hv/kbT吧，你没有带上温度吗？

有温度，是图里的分母上的2。不好意思，太久了忘记了。请问(1)频率单位是hz吗？用hz应该是这样，结果不需要除以3，是359 N/m.之前不知道单位改来改去乱了。(1)是无量纲量，我不太确定应该用什么单位。

AutMaple

这个应该是u=hv/kbT吧，你没有带上温度吗？

高阁

AutMaple 发表于 2025-4-25 23:16
（7）北京高压的刘锦老师吗？你们组只有你一个做计算有点麻烦。。。
你去求求刘锦老师，说想去地化所的 ...

(7)是高科的刘老师；
(6)感谢提醒；
(5)我看网上说U值依赖软件和赝势(不知道是否依赖体系、温压)，找到的文献也都用的QE。我只找到一篇文献推荐了U，使用后电子步很难收敛，之前有人提醒我这说明结构可能不合理。请问应该怎么使用推荐的U呢？
(3)不是声学支，是声子谱的条数。原胞里有n个原子就有3n条声子谱，用j标记。222超胞有8个波矢k，其它波矢k是由周期边界条件引入，非物理。巨正则系综配分函数是对量子态i连乘，i应该既包括k又包括j，我觉得应该把8组连乘起来开8次方以归一化，但是不太清楚怎么在高对称路径上把物理k找到。您说的加权这个是怎么找到物理k并确定权重呢？
(2)这个脚本就是根据(1)算的，可以简述一下：同位素声子谱计算用的就是您讲的方法。phonopy给的band.dat是从左到右再从右到左……这样作图，脚本里先给它顺序调整一下再排成一列。轻重元素都这样就排成了两列，我默认两个声子谱排列方式是对应的，再严格照(1)计算，频率单位取了hz。Γ点声学支会导致分母上3个0，就乘了1。之前尝试用(6)算，phonopy对pdos归一化后为3，需要除掉。请问用(1)算是否也要除以3呢？除以3后似乎很接近正确值。
(1)用公式(1)得到的是β，但是实验测到的是F，大家作图也是F，所以我想算个F好做比较。目前感受是算β再得F准一点，因为用了温度，用(6)算F不准，因为没加温度，应改用准简谐近似声子谱。
此外用(1)还是有些问题，我之前说的是[size=1em]https://doi.org/10.1016/j.gca.2020.11.025，它说β几乎与1/T^2成正比。(1)中β是0K的，用它算2000K下F应该还是有偏差。

感谢指导！

AutMaple

高阁 发表于 2025-4-25 20:46
(1)我用的phonopy算的声子谱，之前算想算β再算错差的很多，目前是通过pdos对Fe积分算。偏差不算离谱，但 ...

（7）北京高压的刘锦老师吗？你们组只有你一个做计算有点麻烦。。。
你去求求刘锦老师，说想去地化所的刘耘老师那里问一问怎么做这些计算。刘耘老师那里，我们有专门适合VASP计算β的脚本，直接处理vasprun.xml文件，基于公式（1）。关于高压的问题，你可以顺便问问地化所的张飞武老师。这些东西在网上交流起来不太方便，但是线下就方便得多。
（6）一系列的实验晶格参数和实验频率，和这些实验晶格参数用软件计算出的压力，两个压力能得到一个依赖关系，这个方法比较糙。细致一点的方法可以问问张飞武老师。
（5）没有晶格参数，就最好用别人推荐的U值
（4）好吧哈哈哈
（3）声学支？不是光学支吗？使用公式（1）计算的时候，带入的都是光学支的频率。原胞的声学支是3，2*2*2就是8个原胞，也就是24个声学支，也没问题呀。
（2）我看不太懂matlab的代码，我只能读懂python的，而且我想看的是你基于公式（1）给出的计算脚本
（1）先算β，再用β算<F>?是不是本末倒置了？

高阁

AutMaple 发表于 2025-4-24 21:42
我经常看Dauphas用公式（6）去算，这算是他们组的传统吧。因为Dauphas之前是做实验的，看他12年的文章， ...

有篇文章说积分的那个公式可以变形，变成NRIXS测的数据积分，那样敏感度会更低。我猜他们这么做是为了和实验数据对比。

高阁

AutMaple 发表于 2025-4-25 16:03
（1）你用什么软件做的频率分析？
（2）你算RPFR的时候用的什么脚本，能给我看一下吗？
（3）你可以只 ...

(1)我用的phonopy算的声子谱，之前算想算β再算<F>错差的很多，目前是通过pdos对Fe积分算<F>。偏差不算离谱，但是我结构还有问题并且没用到T；
(2)我用的matlab，写得很丑，麻烦了；
(3)实验值在300K约为330 N/m，用pdos算出的值在390N/m左右，通过gamma点数据算β因子最终算出的结果为954 N/m，但是我不太清楚是否要除以3，pdos归一化积分要除以3。另有一个问题想请教您，β因子是靠巨正则系综的配分函数算的，为什么下标i不是k与j的合并指标呢？我如果原胞8个原子，扩胞为2*2*2就该有8个物理k，24条声子支，为什么是gamma点24个频率连乘而不是8*24个频率连乘呢？
(4)这个好像我记错了；
(5)我没有晶格参数，请问能不能通过<F>来调整U？
(6)我想请教一下这个怎么矫正呢？实验9对应计算10，是不是计算20就对应实验18呢？
(7)我们组就我一个人做计算。
感谢指点！

AutMaple

高阁 发表于 2025-4-24 23:52
用(1)我不太清楚怎么把轻重同位素的i对好了，我直接从下往上一个个比，结果似乎不对。还有晶体中指标除了 ...

（1）你用什么软件做的频率分析？
（2）你算RPFR的时候用的什么脚本，能给我看一下吗？
（3）你可以只用gamma点的频率，或者调一些高对称点，按照每个点的权重，都算一遍RPFR再加权平均
（4）哪个文献说的力场数应该和温度成正比？
（5）你用DFT+U的方法，要是没有数据，这个对RPFR的影响还是有一些的。我对针铁矿试过，不同的U值对应不同的分馏系数，不过我当时关注的是针铁矿的氧。所以，要是没有推荐的U值数据，你就调你的U值，拟合低压的晶格参数，要是连晶格参数都没有，那你就别调整U了。我之前给你的建议主要是针对计算的晶格参数对实验晶格参数的拟合，调节的是结构优化里的压力设置。
（6）看你提到POSCAR，你应该是用的VASP。要么你加压，拟合一套计算压力和实验压力的对比，比如说，你设置10bar的压力，给出的结构对应实验9bar的压力，那么你就知道计算压力和实验压力的区别，然后拟合一套校正线，用于设置高压下应当设置多少的压力。
（7）你是谁的学生？你们组没有人带你做计算么？你老师就让你自己摸索吗。。。心有点大。

高阁

AutMaple 发表于 2025-4-24 21:18
公式（1）的推导过程就表示了，这个公式是在简谐近似下获得的，所以带入的频率都是简谐频率，简谐频率是 ...

温度我看有文献说键力常数<F>应该和T成正比。我如果直接用0K下的声子谱带入(6)得到的<F>对应T=0，但是公式应该要求温度要足够高才对。我看到说phonopy算声子谱只与提供的unitcell结构有关，如果我提供2000K下的结构就可以算出准简谐近似的声子谱请问是这样吗？如果是简谐近似，热膨胀应该不改变结构，请问我能不能用third-order Birch-Murnaghan isothermal equation算状态方程，得到V然后等比例算结构再算声子谱呢？还有请问有限位移法算声子谱需要单点计算POSCAR-0**结构，这个时候需要加压吗？不好意思，才入门，问题比较多。

高阁

AutMaple 发表于 2025-4-24 21:42
我经常看Dauphas用公式（6）去算，这算是他们组的传统吧。因为Dauphas之前是做实验的，看他12年的文章， ...

用(1)我不太清楚怎么把轻重同位素的i对好了，我直接从下往上一个个比，结果似乎不对。还有晶体中指标除了声子分支数i还有波矢k，但是我不知道用哪个k。请问这两点怎么确定呢？如果我采用DFT+U计算调整U值能算准吗？我没有低压数据……

AutMaple

AutMaple 发表于 2025-4-24 21:18
公式（1）的推导过程就表示了，这个公式是在简谐近似下获得的，所以带入的频率都是简谐频率，简谐频率是 ...

我经常看Dauphas用公式（6）去算，这算是他们组的传统吧。因为Dauphas之前是做实验的，看他12年的文章，他们似乎用实验获得声子谱，如果你单纯从计算的角度出发，用公式（1）就可以了。

AutMaple

高阁 发表于 2025-4-17 21:11
我想从(7)算(6)，原始公式指的是(1)，我不知道怎么从声子谱里找到需要的ω并对应上。感谢指导！

公式（1）的推导过程就表示了，这个公式是在简谐近似下获得的，所以带入的频率都是简谐频率，简谐频率是不和温度相关的。因为它来自于Hessian矩阵，Hessian矩阵是能量对原子坐标的二阶偏导矩阵，能量是电子能，和温度无关。所以你不能用所谓的任意温度下的简谐声子谱。
公式（6）这个和Bigeleisen和Mayer的1947年的文章中的高温近似式是一样的，它是（1）的高温近似，这个你也知道。你为啥要从（7）去计算（6）呢？直接从公式（1）不能计算么？
另外，铁同位素分馏最关键的不是温度，而是你的结构对不对。核幔边界下的压力这么高，没有实验数据对比，单纯加压算出来的结构是有问题的，你最好从一系列的低压的实验结构来外推。

高阁

我想从(7)算(6)，原始公式指的是(1)，我不知道怎么从声子谱里找到需要的ω并对应上。感谢指导！

高阁

AutMaple 发表于 2025-4-15 00:16
你用Dynaphopy计算的非谐频率，是指的对应的高温下简谐振动项不能解释的热膨胀和热传导，而考虑的非谐性 ...

谢谢老师回答！
我算的是地球深部的铁同位素分馏，在2000K下确实可以忽略非简谐效应。我想用<F>=M/\hbar^2\int_0^{+∞}E^2g(E)dE和lnβ=(1/M_1-1/M_2)\hbar^2/(8k^2T^2)。在准简谐近似成立下它和Urey模型的原始公式应该是等价的，但是声子谱得到后原始公式中的ω我不太找的准。我不太清楚BM公式到底指哪个，应该和我现在用的等价。
除了分馏系数，<F>也是我想要的值，但是phonopy得到的声子谱是0K下的。请问怎么得到任意温度下的简谐声子谱呢？我原本想既然非简谐效应能忽略，就用dynaphopy得到2000K下的非简谐声子谱代替简谐声子谱。请问这样是否可行？

AutMaple

高阁 发表于 2025-3-17 00:15
请问能否给出一些参考教程或者是参考文献呢？如果只修改POMASS可否作为一个简单的替代方案呢？此外我用Ur ...

你用Dynaphopy计算的非谐频率，是指的对应的高温下简谐振动项不能解释的热膨胀和热传导，而考虑的非谐性项？
这种计算出来的非谐频率是不能直接带入Urey Model或者BM公式的。具体可以看地化所刘琪副研究员2010年的文章“On the proper use of the Bigeleisen–Mayer equation and corrections to it in the calculation of isotopic fractionation equilibrium constants”。
使用phonopy改变质量计算的计算FORCE_CONSTANTS是没有变化的，如果这个力场数不是经过一些质权处理后的力常数的话，具体可以搜一搜sob关于同位素替代、hessian矩阵、力场数矩阵和质权的介绍文章，他写的很详细。
高温下的2000k的非谐是不需要考虑的，这个在Dupuis的文章“Importance of a Fully Anharmonic Treatment of Equilibrium Isotope Fractionation Properties of Dissolved Ionic Species As Evidenced by Li+”可以看出来，高温下使用PIMD计算的分馏系数和简谐近似下的结果是没啥差别的。
我建议你直接使用BM公式计算2000K下的分馏系数，不需要考虑非谐效应。就算你计算的是含氢体系，你也不需要考虑非谐性。至于为什么不能用POMASS替换来计算β值或者RPFR值，这是因为两次计算会产生一些数值误差，最后获得的RPFR值很奇怪。之前厦大的一位师兄就遇到了这样的问题。你直接使用phonopy产生力场数，在这个力常数的基础上，使用两套质量，让MS给你找一套高对称的点，让phonopy输出两套频率，再带入BM公式就可以计算出β值了，至于温度的依赖性BM公式里是有的。
如果你要问从高温的MD，到获取构型，计算频率，最后计算RPFR值的这一套流程，你可以看现在在成都理工大学的高才洪研究员的文章，First-principles calculations of equilibrium bromine isotope fractionations，或者科大的黄方课题组出来的王文忠老师早年的文章“Wang,2019-Equilibrium Mg isotope fractionation among aqueous Mg2+, carbonates, brucite and lizardite Insights from first-principles molecular dynamics simulations”。
最后，你们组为啥要计算高温的分馏呀？2000k大部分分馏都没有啦。

高阁

AutMaple 发表于 2022-3-18 10:56
如果想计算平衡同位素分馏，我猜测你是想使用Urey模型或者BM方程来计算。因此需要用到同种物质下，不同质量 ...

请问能否给出一些参考教程或者是参考文献呢？如果只修改POMASS可否作为一个简单的替代方案呢？此外我用Urey模型计算平衡同位素分馏需要加温到2000K，还需要用Dynaphopy进一步计算。请问是否用phonopy改变质量计算FORCE_CONSTANTS，不用改变MD模拟就足够了呢？