计算同位素如何设置

asl1994

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

jiangning198511

把同一类元素在POSCAR里面分开，在POTCAR分别用修改质量的和没有修改的
例如你要计算H2，POSCAR 中写为
1 1
....
....

POTCAR
H
....
D
...
这样两种都有了

asl1994

jiangning198511 发表于 2017-8-13 10:53
把同一类元素在POSCAR里面分开，在POTCAR分别用修改质量的和没有修改的
例如你要计算H2，POSCAR 中写为
1 ...

明白了，谢谢！

AutMaple

如果想计算平衡同位素分馏，我猜测你是想使用Urey模型或者BM方程来计算。因此需要用到同种物质下，不同质量替换后的振动频率。我不建议直接修改POTCAR。亦可以使用PHONOPY来替换同位素质量，其原理是使用vasprun.xml中输出的hessian在解质权之后，加上你所设置的质权，进而输出不同同位素质量组合下的震动频率。其实按照Urey模型的假设，同位素质量的变化并不会导致其结构发生明显的变化。因此，你只需要用一个质量优化出结构，用有限位移法或者DFPT计算出hessian矩阵就行了。不过注意，如果你计算的结构固定了一部分原子，PHONOPY会出错，因为这时你计算的hessian矩阵维度和你的结构中的原子数目对不上，这时候你就需要手动写一个code来解开质权，加上同位素质量的质权，解开矩阵获得特征值计算质量了。

高阁

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

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

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 发表于 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下的非简谐声子谱代替简谐声子谱。请问这样是否可行？

高阁

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

AutMaple

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

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

AutMaple

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

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

高阁

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

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

高阁

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 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-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 21:42
我经常看Dauphas用公式（6）去算，这算是他们组的传统吧。因为Dauphas之前是做实验的，看他12年的文章， ...

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