Hi > Dear Prof. Hutter: > We are currently testing chemical reaction by CPMD. I have some > puzzles with the following questions: > > 1. How to construct pseudopotentials and pseudowavefunction for cation > and anion? > We have tested to use neutral atomic pseudopotential, for example > Cl, and use the "ATOMIC CHARGE" with value -1.0 to run Cl- anion. > However we always get the same neutral result. Should we use anion Cl- > pseudopotentials? The basic idea is to have one PP for all possible states of an atom (transferablility) including anion and cations. Of course you can also construct PP using an ion as reference state. In a ab initio program you only specify the total number of electrons of the whole system. Default is to use the number of electrons to achieve charge neutrality. This can be changed using the CHARGE keyword. ATOMIC CHARGE is only used for the initial guess and can usually be ignored. > > 2. In your paper, Chem Phys. Lett., 321(2000)225, about proton transfer > how to generate H+ cation pseudopotential? Use the normal H PP (see above) > According to Fig. 2 in that paper, it seems to me all of the atoms of > (H2O)3H+ move. Does H+ transfer to water 3?? > There is no distiction of H and H+. We see proton transfer several times. > 3. If the answer for question 1 is that we have to use negative charged > pseudopotentials, then this question is puzzle to us. > We want to study Cl- anion(incoming) + ClCCH2CN, S_N2 reaction, the > incoming Cl- has to use negative charged pseudopotentials and the > outgoing Cl at beginning has to use neutral pseudopotentials. After the > reaction the incoming particle has to use neutral pseudopotential and > the outgoing Cl has to use negative charged pseudopotentials. This > question is similar to question 2. > see above > If the answer for the anion and cation pseudopotentials should all use > neutral pseudopotentials, then how to make "ATOMIC CHARGE" works? > see above > 4. We want to study another S_N2 reaction of O(1D) +CH4, where 1D is the > term symbol. The oxygen is in its excited state. Can we use tddft to > excite oxygen atom to its excited state and then write out its > wavefunction with fixed term symbol. How and which files story this > information? > Then the next question is to read in the restart file for oxygen only > and put CH4 as initial run. Can this work? and how? > This is not possible. regards Juerg Hutter
Sunday, November 18, 2012
can cpmd run for chemical reaction?
The answer is yes. The following information is from CPMD list by Prof. Hutter.
动力学模拟： nvt, npt and nve
转载文章， 抱歉，原作者不详。
第一原理动力学的系综选择，和各种系综的输入文件Example
做Car-Paranello MD，无论是用VASP, DMOL还是CPMD/cp2k， 都面临着如何选择系综的问题。那该怎么选择呢？我的看法是，应该尽可能使用NVE系综，而不是NVT/NPT。为什么呢？因为NVE系综是真正体系能量守恒的。
像 使用Nose-Hoover chain 的NVT系综固然也能保证能量守恒,但是它使用一个或多个虚拟的储热罐来缓存由温度波动引起的能量波动，从而模仿一个NVE系综。这种扩展拉格朗日体系方 法对体系性质的热力学平均和采样的影响不大，但是，对具体的动力学性质的细节会产生一定影响，因为虚拟的储热罐的动力学行为随着其本身耦合频率，幅度等参 数的调整，对体系的演进有着一定影响，具体就是对分子的扩散，转动，化学反应的时间特征影响。另外，任何的拉格朗日虚拟扩展量都会增加实际的运动自由度， 所以会略微慢点。
NVE则没有这些缺点，它是完美的隔离体系，完全遵循能量守恒定律的条件，也没有虚拟的动力学变量。
但 NVE有个缺点妨害了它的使用。什么缺点呢？就是的它的能量守恒于起始能量。即守恒能量等于起始构型的势能+起始设定的动能。很可能就是你前面预平衡阶段 最后一帧的速度和构型。NVE开始后，能量在动能、势能之间波动式分配。如果一开始势能高，能量就向动能转化，平衡后温度就高；反之，动能向势能转化，平 衡后温度就低。显然，这更接近于真实体系的动能-势能的转化模式，但是这就使得其最终的温度其平衡之前很难预料，你本来想模拟300K，结果最终温度却稳 定在了250K, 或者400K！你怎么办？只好使用NVT, NPT去控制温度。
从上面分析看，使用NVT/NPT是不得以而强为之。还是NVE好。
NVE的温度问题也不是不能解决。一般用Velocity scaling积分算法跑一个短时间的NVT模拟，比如100fs，取最后一帧Restart一个NVE，一般都能把平衡后的温度控制在一定的范围内。
附， CPMD各种系综的设定。
1. NVE呢，不指定任何控温关键词的，又没有设定任何控压关键词的，即为NVE。For example:
&CPMD
MOLECULAR DYNAMICS CP
RESTART WAVEFUNCTION COORDINATES VELOCITIES CELL LATEST
EMASS
500
TIMESTEP // 0.096fs/步
4
MAXSTEP
50000
&END
如果不想使用前面模拟最后一帧的速度，可以加关键词：
TEMPERATURE
300
它的意思是在NVE的开始前，重新为每个原子/动力学变量设定动能，使得总的温度为300K。NVE开始后，并不控制体系守在这个温度。
2. NVT 有多种多种控制温度的关键词。
(1) Velocity Scaling控温可以帮助体系快速预平衡到设定温度/势能附近，但是不满足能量守恒，一般不用于采样阶段的模拟。关键词： TEMPCONTROL
&CPMD
MOLECULAR DYNAMICS CP
RESTART WAVEFUNCTION COORDINATES VELOCITIES CELL LATEST
TEMPCONTROL IONS, ELECTRONS
300 0.006
EMASS
500
TIMESTEP
4
MAXSTEP
50000
&END
其中，IONS 300设定原子的温度为300K; ELECTRON 0.006(au), 是设定电子，即波函数作为运动量(或叫扩展拉格朗日变量)的温度，需预先跑个100步的NVE，取最后几步的EKINC即电子动能，作为这里的控温参数。可比之稍大一点。
(2) NOSE-HOOVER Chain控温是遵循能量守恒的NVT系综，或者说它通过增加变量的方法让体系模拟一个隔离的NVE体系。第一原理模拟中广泛使用。
&CPMD
MOLECULAR DYNAMICS CP
RESTART WAVEFUNCTION COORDINATES VELOCITIES CELL ACCUMULATORS NOSEE NOSEP LATEST // 重接上次模拟继续
TIMESTEP
4
MAXSTEP
10000
TRAJECTORY XYZ SAMPLE //每隔10步采个样存到TRAJECTORY中
10
DIPOLE DYNAMICS WANNIER SAMPLE //每隔10步存储Dipole moment到DIPOLE文件中
10
STORE //每隔500步写一次RESTART.1文件
500
EMASS
500
NOSE IONS // NOSE控制原子温度和频率(单位波数，即cm-1)
300 1600
NOSE ELECTRONS // NOSE控制电子动能和频率(单位波数，即cm-1)
0.020 25000
&END
(3) NPT的设定：NPT的目的是最终使得盒子在每个方向上的受力与外部的压力相等。压力控制要在&CPMD部分使用PARRINELLO- RAHMAN NPT关键词，并在&SYSTEM中设定PRESSURE的大小。此外，还有一些其它的关键词控制压力的各向异性，请查阅手册。
&CPMD
MOLECULAR DYNAMICS CP
RESTART WAVEFUNCTION COORDINATES VELOCITIES CELL ACCUMULATORS NOSEE NOSEP LATEST // 重接上次模拟继续
TIMESTEP
4
MAXSTEP
10000
TRAJECTORY XYZ SAMPLE //每隔10步采个样存到TRAJECTORY中
10
DIPOLE DYNAMICS WANNIER SAMPLE //每隔10步存储Dipole moment到DIPOLE文件中
10
STORE //每隔500步写一次RESTART.1文件
500
EMASS
500
NOSE IONS // NOSE控制原子温度和频率(单位波数，即cm-1)
300 1600
NOSE ELECTRONS // NOSE控制电子动能和频率(单位波数，即cm-1)
0.020 25000
PARRINELLO-RAHMAN NPT //设定NPT系综和控压方法
STRESS TENSOR //每隔100步输出一次压力值
100
&END
&SYSTEM
ANGSTROM
PRESSURE // 设定外部压力为100Kbar
100
CELL VECTORS // NPT模拟需要使用CELL VECTOR来设置盒子
13.080 0.000 0.000
0.000 11.050 0.000
-9.805 0.000 14.374
CUTOFF
100.0&END
Thursday, November 8, 2012
分享发高水平文章的经验-转载
杨磊的文章。
一觉醒来，发现有这么多关注，有些惶恐，为避免误导，把后记改为前言。
有不对之处，敬请指正，更多分享请参见我在新浪的博客
当 你做完大部分实验或计算之后，就要开始着手写论文了。对于Natured子刊、JACS和Advanced Materials这类杂志来说，论文撰写的重要性我觉得至少占40%。也就是说如果你能够切入一个非常有吸引力的角度，你可以让你的实验结果发到更好的 杂志。对于NS来说，我觉得实验的设计更为重要。如何能够写好一篇文章，我认为首先应该抛弃两个错误的看法。第一：不要鄙视烂的结果都能够发在好杂志上。 你需要思考如果你拿这些数据能够把文章写成怎样。你要学习你没有想到的“点”。比如说，性能可能并没有非常突出，但是他/她提出了一个非常有启发性的假 设。第二：不要认为审稿人误会你的评语愚蠢。我知道审稿人在审阅时（包括我在审Advanced Materials时）速度是非常快的。如果一个领域的评审人在短时间内都没有看出你的创新点，说明你没有表达清楚。我经常听到有人抱怨“我这篇文章其实 和以前不一样，审稿人却认为没有新东西”或者“我的性能明显要比别人的文章好，不知道为什么审稿人没有注意到”等等。出现这种情况后，要重新审视自己的文 章。思考怎样写别人不会忽视我的重点，怎样写不会让人误解。一个小窍门是让你的同学（大方向一致但不是一个小领域的）快速浏览一下你的文章，让他指出不确 定的东西，然后加以改正。
我觉得写文章最重要也最难写的就是Introduction。这是审稿人看得比较认真而且容易理解的部分。 而且我发现一个规律，越好的杂志，审稿人越喜欢攻击introduction。可能是因为你的实验设计已经很好，不太容易有问题。但是对于 introduction，审稿人却非常容易下手。比如这篇文章没有新意，或者你在introduction提到的问题，在正文中没有解决等等。在读好文 章时一定要学习他们在组织introduction时的思路。其次，一定要有一个吸引人的标题。不要过于中立。我以前投一篇文章的时候，刚开始拟定为 Sulfur Poisoning Behavior of ....。后来偶然看到Berkeley物理系的一片不相干的文章，用了New Insights into ..。我就把这个模式套用到我的文章上，我导师认为这个标题立马让文章档次提高。我的一个经验，经常收集那些好文章的title （不需要局限你的领域），以备将来时灵活运用。至于正文，只要围绕你的Introduction，反复强调你的创新性（一定要“反复”，因为审稿人会忽 视），一般没有什么问题。另外，因为审稿人是带着寻找问题的模式去评判文章的，所以在正文中的每一句话不要过度发散，否则很容易招致不严谨或者补充数据的 评语。
一觉醒来，发现有这么多关注，有些惶恐，为避免误导，把后记改为前言。
有不对之处，敬请指正，更多分享请参见我在新浪的博客
前
言：这三个部分分享了很多关于提升自己成果的经验，但是大家不要进入一个误区：为了发文章而做实验。
发牛文更多是因为你的研究热情和辛勤付出，因为科研成果的内核还是你能否真正解决前人未能解决的问题。当然，从营销学角度，我们去探寻并运用这些规律也是
无可厚非的。在当今商业高度发达的年代，你不仅要有过硬的产品（真正的实验突破），也要学会去包装和做增值服务（系统性地设计实验），更要懂得去推广自己
的产品（写一篇好的文章）。
一 选择自己领域中热门和能够突破的方向
要 想在Nature 或者Science （以下简称NS）上发表文章，首先要对自己领域最近10年有哪些文章发表在这些刊物上，并进行分类。以氧化物燃料电池领域为例，在2002-2012年区 间总共有8篇文章发表在这两个杂志上。如果你研究的小领域没有文章在NS或者Nature的子刊上，那说明杂志编辑认为你的领域不具有很广的关注度。同 时，要分析是些什么样的突破发表在NS上。比如在这8篇文章，有6篇文章直接与燃料电池的阳极材料有关。这说明如果你能在阳极的研究中有所突破，存在在 NS上发表的可能性。再进一步分析其细节，你会发现更多的规律。比如，燃料电池阳极的最主要的问题是碳氢燃料在高温下的裂解导致碳沉积和硫在镍表面吸附导 致阳极硫中毒。早前的SN上的文章主要关注怎样防止在阳极上的碳沉积，在2006年首先出现了一个新的阳极材料同时有抗碳沉积和抗硫中毒。这篇文章给了我 一个启发，说明现有的阳极材料必须能够同时解决这两个问题，才有可能在NS上出现。当然这也是合理的，因为碳氢燃料包含碳和硫。
当 然，并不是说你知道这些趋势，你一定能够在这样上面有所突破，但是能够给你一个非常具有指引性的思路。比如说，当时我的研究课题是做电解质的，因为师兄毕 业需要移交阳极的课题，我学习了一段时间。我把我所研究的新电解质去做阳极的抗硫测试，发现具有不可思议的抗硫性能。在我多次重复加以确认之后，我意识到 了其重要性。其实当时有人建议说可以用这个结果在Advanced Materials上投一篇文章，但是在我分析这些年在SN上发表的氧化物燃料电池文章，我决定继续研究该阳极的抗碳沉积特性，然后进一步优化。这个做法 非常重要，使我后来冲击Science奠定了重要的基础。
要 想在Nature 或者Science （以下简称NS）上发表文章，首先要对自己领域最近10年有哪些文章发表在这些刊物上，并进行分类。以氧化物燃料电池领域为例，在2002-2012年区 间总共有8篇文章发表在这两个杂志上。如果你研究的小领域没有文章在NS或者Nature的子刊上，那说明杂志编辑认为你的领域不具有很广的关注度。同 时，要分析是些什么样的突破发表在NS上。比如在这8篇文章，有6篇文章直接与燃料电池的阳极材料有关。这说明如果你能在阳极的研究中有所突破，存在在 NS上发表的可能性。再进一步分析其细节，你会发现更多的规律。比如，燃料电池阳极的最主要的问题是碳氢燃料在高温下的裂解导致碳沉积和硫在镍表面吸附导 致阳极硫中毒。早前的SN上的文章主要关注怎样防止在阳极上的碳沉积，在2006年首先出现了一个新的阳极材料同时有抗碳沉积和抗硫中毒。这篇文章给了我 一个启发，说明现有的阳极材料必须能够同时解决这两个问题，才有可能在NS上出现。当然这也是合理的，因为碳氢燃料包含碳和硫。
当 然，并不是说你知道这些趋势，你一定能够在这样上面有所突破，但是能够给你一个非常具有指引性的思路。比如说，当时我的研究课题是做电解质的，因为师兄毕 业需要移交阳极的课题，我学习了一段时间。我把我所研究的新电解质去做阳极的抗硫测试，发现具有不可思议的抗硫性能。在我多次重复加以确认之后，我意识到 了其重要性。其实当时有人建议说可以用这个结果在Advanced Materials上投一篇文章，但是在我分析这些年在SN上发表的氧化物燃料电池文章，我决定继续研究该阳极的抗碳沉积特性，然后进一步优化。这个做法 非常重要，使我后来冲击Science奠定了重要的基础。
二、系统性的设计实验
据 我了解，很多最为关键或者突破的实验数据都是意外得到的，或者超过自己预期的 （当然也存在像Goodenough教授这种牛人能够从理论上设计材料）。当你获得比以前文献中更好的性能时，就要开始考虑怎么设计一系列系统的试验，以 能够将来写出一篇有完整故事情节的文章。因为现在已经不是“酒香不怕巷子深”的年代了，除非你的结果能够改变人类的认知，否则都需要思考围绕该突破的实验 设计。其工作量大约是一般长文的2~3倍。除了最为关键的4个图放在正文，其余的将放到补充材料里面。
实验该怎么设计才会对主编和审 稿人的口味？当然不同领域有不同的文章结构。一个简单的方法就是你尽可能把自己领域中不同小方向在Journal of the American Chemical Society, Angewandte Chemie International Edition 和Advanced Materials 上面的文章综合起来。比如，这些杂志上有专注于合成的、有专注于表征的或者专注机理理解的文章。你如果能够把这些文章的最有特色的东西有机的糅合在一起， 你的文章就已经高于这些杂志的档次了。以催化和表面化学为例，SN上的实验设计思路一般来说就是一个比较新颖的纳米结构，比较高档的表征（如STM或同步 辐射）、优异的性能和分子动力学的理论计算。如果你去详细比较SN上某一篇文章每一小部分和JACS上类似的的全文，你会发现其实JACS上的水平更专。 根据这个思路，你就可以设计完整的实验，寻找合作对象，相互促进，最终达到一个完美的实验结论。我的那篇Nature Communication 就是以这种思路设计的。当时需要对我们现有性能的理论解释，我们寻求了与布鲁克海文国家实验室的合作。他们给我们提供了很好的思路，继续优化实验，与他们 的理论达到了较好的融合。虽然在投Nature主刊40多天后被拒，但是审稿人对实验设计非常肯定：This paper has really nice science；The science is top notch等等。这篇文章本身的实验结果没有我Science上那篇文章的突破大，但是好的实验设计让这篇文章被子刊接受。
三、撰写完整且吸引人的文章据 我了解，很多最为关键或者突破的实验数据都是意外得到的，或者超过自己预期的 （当然也存在像Goodenough教授这种牛人能够从理论上设计材料）。当你获得比以前文献中更好的性能时，就要开始考虑怎么设计一系列系统的试验，以 能够将来写出一篇有完整故事情节的文章。因为现在已经不是“酒香不怕巷子深”的年代了，除非你的结果能够改变人类的认知，否则都需要思考围绕该突破的实验 设计。其工作量大约是一般长文的2~3倍。除了最为关键的4个图放在正文，其余的将放到补充材料里面。
实验该怎么设计才会对主编和审 稿人的口味？当然不同领域有不同的文章结构。一个简单的方法就是你尽可能把自己领域中不同小方向在Journal of the American Chemical Society, Angewandte Chemie International Edition 和Advanced Materials 上面的文章综合起来。比如，这些杂志上有专注于合成的、有专注于表征的或者专注机理理解的文章。你如果能够把这些文章的最有特色的东西有机的糅合在一起， 你的文章就已经高于这些杂志的档次了。以催化和表面化学为例，SN上的实验设计思路一般来说就是一个比较新颖的纳米结构，比较高档的表征（如STM或同步 辐射）、优异的性能和分子动力学的理论计算。如果你去详细比较SN上某一篇文章每一小部分和JACS上类似的的全文，你会发现其实JACS上的水平更专。 根据这个思路，你就可以设计完整的实验，寻找合作对象，相互促进，最终达到一个完美的实验结论。我的那篇Nature Communication 就是以这种思路设计的。当时需要对我们现有性能的理论解释，我们寻求了与布鲁克海文国家实验室的合作。他们给我们提供了很好的思路，继续优化实验，与他们 的理论达到了较好的融合。虽然在投Nature主刊40多天后被拒，但是审稿人对实验设计非常肯定：This paper has really nice science；The science is top notch等等。这篇文章本身的实验结果没有我Science上那篇文章的突破大，但是好的实验设计让这篇文章被子刊接受。
当 你做完大部分实验或计算之后，就要开始着手写论文了。对于Natured子刊、JACS和Advanced Materials这类杂志来说，论文撰写的重要性我觉得至少占40%。也就是说如果你能够切入一个非常有吸引力的角度，你可以让你的实验结果发到更好的 杂志。对于NS来说，我觉得实验的设计更为重要。如何能够写好一篇文章，我认为首先应该抛弃两个错误的看法。第一：不要鄙视烂的结果都能够发在好杂志上。 你需要思考如果你拿这些数据能够把文章写成怎样。你要学习你没有想到的“点”。比如说，性能可能并没有非常突出，但是他/她提出了一个非常有启发性的假 设。第二：不要认为审稿人误会你的评语愚蠢。我知道审稿人在审阅时（包括我在审Advanced Materials时）速度是非常快的。如果一个领域的评审人在短时间内都没有看出你的创新点，说明你没有表达清楚。我经常听到有人抱怨“我这篇文章其实 和以前不一样，审稿人却认为没有新东西”或者“我的性能明显要比别人的文章好，不知道为什么审稿人没有注意到”等等。出现这种情况后，要重新审视自己的文 章。思考怎样写别人不会忽视我的重点，怎样写不会让人误解。一个小窍门是让你的同学（大方向一致但不是一个小领域的）快速浏览一下你的文章，让他指出不确 定的东西，然后加以改正。
我觉得写文章最重要也最难写的就是Introduction。这是审稿人看得比较认真而且容易理解的部分。 而且我发现一个规律，越好的杂志，审稿人越喜欢攻击introduction。可能是因为你的实验设计已经很好，不太容易有问题。但是对于 introduction，审稿人却非常容易下手。比如这篇文章没有新意，或者你在introduction提到的问题，在正文中没有解决等等。在读好文 章时一定要学习他们在组织introduction时的思路。其次，一定要有一个吸引人的标题。不要过于中立。我以前投一篇文章的时候，刚开始拟定为 Sulfur Poisoning Behavior of ....。后来偶然看到Berkeley物理系的一片不相干的文章，用了New Insights into ..。我就把这个模式套用到我的文章上，我导师认为这个标题立马让文章档次提高。我的一个经验，经常收集那些好文章的title （不需要局限你的领域），以备将来时灵活运用。至于正文，只要围绕你的Introduction，反复强调你的创新性（一定要“反复”，因为审稿人会忽 视），一般没有什么问题。另外，因为审稿人是带着寻找问题的模式去评判文章的，所以在正文中的每一句话不要过度发散，否则很容易招致不严谨或者补充数据的 评语。
Tuesday, November 6, 2012
cpmd 系统能量收敛问题
【原创】CPMD NOSE-HOOVER NaN问题
作者: ChemiAndy
能量NaN是说数值太大，写不下了。问题三个原因：EMASS选择不当，TIMESTEP选择不当，和NOSE ELECTRON的EKIN选择不当。
所有三个原因，都归结于一个原因，体系物理总能和虚拟电子动能的之间的绝热出了问题，这个是CPMD的几乎所有问题的根源。
本来呢，体系的哈密顿应该等于体系物理总能(ECLASSICAL)=离子动能+势能(Kohn-Sham)，介个是从拉格朗日原理推导出来的，Shordinger方程就是照这个写的。
但是在CPMD中，
Hamitonian(EHAM)=体系物理总能(ECLASSICAL)+虚拟电子动能(EKIN)
多出来一项。多一项不要紧，拉格朗日方程L(q,p)并没有限定说什么不能作为体系变量。所以嘛，我们给体系增加一个自由度，相当于增加一个约束而已。设计系综的人经常这么搞。NOSE-HOOVER啊啥的，全都是增加虚拟自由度，以实现对体系的动能约束。这叫扩展拉格朗日方法。问题是，CPMD为神马要加这一项呢？虚拟电子动能， DFT中的Kohn-Sham方程中已经包括电子动能项了啊（电子密度的函数），为啥又要搞个动能项出来？
这正是Car & Prrinello 的Niubility的地方。MD主要问题是啥？算力啊。算出来每个原子受力，好知道下一步跑到哪。用DFT算原子受力，得先算波函数，做SCF自洽场计算，波函数对坐标一阶导，就是力啦。可SCF这玩意儿老费时间，算不起啊，怎么办呢？C&P俩老哥琢磨着，原子核跑，电子也在跟着跑啊。电子是啥？就是波函数。电子往哪边跑一下，对平面波函数而言，无非就是单电子波函数（轨道）乘以某个系数而已。既然这样，还算啥波函数啊，干脆把电子作为一个运动变量，让它自己跑得了。每一步根据它的受力情况，直接预测下一时间步的位置（轨道系数）。那电子受力怎么算呢？就是总势能泛函(EKS)对单个轨道(psi-i)的一阶导（叫“场力”）。问题是，一个东西没有质量你是不可能计算它受到的力的，F=ma嘛，幼儿园都学过了。电子当然有质量，可是太小了，这样的东西，即使只受到一个极其微弱的力，一个MD时间步之后，都跑到天上去了。不行，得加大这个电子质量，就是虚拟电子质量。只要找到一个合适的虚拟质量，那么电子就会老老实实跟着原子核跑，而不至于跑到天上去。
那你老弟问了，让波函数这样跑，你怎么能保证它下一步出现在它实际应该出现的地方呢？能跟实际优化出来的波函数对的上吗？实践验证，跟BOMD(每步都优化波函数)比较，居然出奇地一致。为什么捏。。。回帖可见。.....(In this terminology, “low electronic temperature” or “cold electrons” means that the electronic subsystem is close to its instantaneous minimum energy min{Φi } E KS i.e. close to the exact Born–Oppenheimer surface. Thus, a ground–state wavefunction optimized for the initial configuration of the nuclei will stay close to its ground state also during time evolution if it is kept at a sufficiently low temperature.)即，只要电子足够冷，它就会贴紧原子核运动，就像你的衣服一样，给它一点动能，它也不会飞走。虚拟电子质量只要在一个合理的范围内，大一点，小一点，对电子动力学行为影响很小。
这样说不是说电子质量的设置不重要，恰恰相反，电子质量是整个CPMD成败的关键。为什么这么说呢？因为尽管我们希望电子很冷，但是伙计你看，它贴着一个高温的原子呢。控制不好，原子就会不断把能量泄漏给电子。波函数乱跑，体系能量就急剧升高。波函数都错了，电子电子叠到了一起，能量还不飙到天上去？ 你就得到一堆NaN
(这就是不每步优化波函数的下场，哼哼)这就是传说中的动能耦合问题，或曰“绝热”问题， adiabaticity. 保持电子与原子核的绝热性，是CPMD的首要问题。
避免不绝热，就要保持电子的质量足够小。道理很简单，原子核和电子之间的能量传递方式是振动。如果电子质量大，它的振动频率就可能低到3000-4000cm-1，而这恰好是H原子的振动频率范围。频率接近，那么两个人就会玩车震。。错了，是共振，能量就在共振中传递。反过来，如果电子质量就一点点，其振动频率就会很高。原子体系对此不感冒。从这儿也可以看出，只要体系都是重原子，电子虚拟质量高一点没关系。在CPMD中，电子虚拟质量EMASS默认设置是400au.但是文献表明，对于水这个体系，EMASS最佳值是270.
电子虚拟质量高一点有啥好处？低一点有啥坏处？问得好。电子虚拟质量高，就可以允许较大的时间步长TIMESTEP。MD的时间步长取决于最高振动频率。比如H原子，每10fs振动一次，那么你的时间步长就应该在1fs左右。否则可能存在一步走过的风险，成为失足青年。一般的体系，如果EMASS选400, TIMESTEP可以是4au(~0.1fs)。但是，如果你的体系有H，而且在我的忽悠下你的EAMSS选了270,那老兄你的TIMESTEP就得选2.5(0.06fs)，否则还是一堆NaN。杯具啊青年，老板给你1个月，你丫机器占了2个月，大伙儿恨不得把你吊起来往你身上滴蜡烛。
下面来说说为什么NOSE ELECTRON中的参数EKIN也会导致NaN杯具。EKIN是电子虚拟动能的热浴值。CPMD中的NOSE HOOVER热浴是扩展拉格朗日体系方法，是最严格的能量守恒的热浴方法。普通的热浴，要么直接Scale原子速度(Andersen)，要么加一个活塞(Beredsen)，原理过于简单，都是属于温度高了移走点能量，温度低了输入进去点能量，不可避免会造成能量恶意不守恒。看起来没问题，系综等概率假定已经失陷，统计热力学贞操已经失守。NOSE-HOOVER怎么搞呢？原子核的温度还是高了砍，低了补，只是增加或减少的能量都来自哈密顿里的一个储热罐(一个带质量，有动能，并乘以系数Q。原子核和这个虚拟物质通过Q传递热量)。这样就能保持总动能一定，无非那边多点这边少点的问题。具体到CPMD，则分别为原子核和电子搞了两个一定容量的储热罐在哈密顿里。电子的虚拟动能热浴值的设置很重要。原理和上面一样，设置过大或过小，都会改变绝热性。想避免不绝热问题，必须选择合适的EKIN。其实也简单，假如要跑一个100K的体系，那么让这个体系跑个初始温度100K不带热浴的NVE，看看EKIN的平均能量是多少。把这个值拿过来做NOSE ELECTRON的EKIN参数即可。至于频率值，大概设在20000-30000就行，不很严肃。
所有三个原因，都归结于一个原因，体系物理总能和虚拟电子动能的之间的绝热出了问题，这个是CPMD的几乎所有问题的根源。
本来呢，体系的哈密顿应该等于体系物理总能(ECLASSICAL)=离子动能+势能(Kohn-Sham)，介个是从拉格朗日原理推导出来的，Shordinger方程就是照这个写的。
但是在CPMD中，
Hamitonian(EHAM)=体系物理总能(ECLASSICAL)+虚拟电子动能(EKIN)
多出来一项。多一项不要紧，拉格朗日方程L(q,p)并没有限定说什么不能作为体系变量。所以嘛，我们给体系增加一个自由度，相当于增加一个约束而已。设计系综的人经常这么搞。NOSE-HOOVER啊啥的，全都是增加虚拟自由度，以实现对体系的动能约束。这叫扩展拉格朗日方法。问题是，CPMD为神马要加这一项呢？虚拟电子动能， DFT中的Kohn-Sham方程中已经包括电子动能项了啊（电子密度的函数），为啥又要搞个动能项出来？
这正是Car & Prrinello 的Niubility的地方。MD主要问题是啥？算力啊。算出来每个原子受力，好知道下一步跑到哪。用DFT算原子受力，得先算波函数，做SCF自洽场计算，波函数对坐标一阶导，就是力啦。可SCF这玩意儿老费时间，算不起啊，怎么办呢？C&P俩老哥琢磨着，原子核跑，电子也在跟着跑啊。电子是啥？就是波函数。电子往哪边跑一下，对平面波函数而言，无非就是单电子波函数（轨道）乘以某个系数而已。既然这样，还算啥波函数啊，干脆把电子作为一个运动变量，让它自己跑得了。每一步根据它的受力情况，直接预测下一时间步的位置（轨道系数）。那电子受力怎么算呢？就是总势能泛函(EKS)对单个轨道(psi-i)的一阶导（叫“场力”）。问题是，一个东西没有质量你是不可能计算它受到的力的，F=ma嘛，幼儿园都学过了。电子当然有质量，可是太小了，这样的东西，即使只受到一个极其微弱的力，一个MD时间步之后，都跑到天上去了。不行，得加大这个电子质量，就是虚拟电子质量。只要找到一个合适的虚拟质量，那么电子就会老老实实跟着原子核跑，而不至于跑到天上去。
那你老弟问了，让波函数这样跑，你怎么能保证它下一步出现在它实际应该出现的地方呢？能跟实际优化出来的波函数对的上吗？实践验证，跟BOMD(每步都优化波函数)比较，居然出奇地一致。为什么捏。。。回帖可见。.....(In this terminology, “low electronic temperature” or “cold electrons” means that the electronic subsystem is close to its instantaneous minimum energy min{Φi } E KS i.e. close to the exact Born–Oppenheimer surface. Thus, a ground–state wavefunction optimized for the initial configuration of the nuclei will stay close to its ground state also during time evolution if it is kept at a sufficiently low temperature.)即，只要电子足够冷，它就会贴紧原子核运动，就像你的衣服一样，给它一点动能，它也不会飞走。虚拟电子质量只要在一个合理的范围内，大一点，小一点，对电子动力学行为影响很小。
这样说不是说电子质量的设置不重要，恰恰相反，电子质量是整个CPMD成败的关键。为什么这么说呢？因为尽管我们希望电子很冷，但是伙计你看，它贴着一个高温的原子呢。控制不好，原子就会不断把能量泄漏给电子。波函数乱跑，体系能量就急剧升高。波函数都错了，电子电子叠到了一起，能量还不飙到天上去？ 你就得到一堆NaN
(这就是不每步优化波函数的下场，哼哼)这就是传说中的动能耦合问题，或曰“绝热”问题， adiabaticity. 保持电子与原子核的绝热性，是CPMD的首要问题。
避免不绝热，就要保持电子的质量足够小。道理很简单，原子核和电子之间的能量传递方式是振动。如果电子质量大，它的振动频率就可能低到3000-4000cm-1，而这恰好是H原子的振动频率范围。频率接近，那么两个人就会玩车震。。错了，是共振，能量就在共振中传递。反过来，如果电子质量就一点点，其振动频率就会很高。原子体系对此不感冒。从这儿也可以看出，只要体系都是重原子，电子虚拟质量高一点没关系。在CPMD中，电子虚拟质量EMASS默认设置是400au.但是文献表明，对于水这个体系，EMASS最佳值是270.
电子虚拟质量高一点有啥好处？低一点有啥坏处？问得好。电子虚拟质量高，就可以允许较大的时间步长TIMESTEP。MD的时间步长取决于最高振动频率。比如H原子，每10fs振动一次，那么你的时间步长就应该在1fs左右。否则可能存在一步走过的风险，成为失足青年。一般的体系，如果EMASS选400, TIMESTEP可以是4au(~0.1fs)。但是，如果你的体系有H，而且在我的忽悠下你的EAMSS选了270,那老兄你的TIMESTEP就得选2.5(0.06fs)，否则还是一堆NaN。杯具啊青年，老板给你1个月，你丫机器占了2个月，大伙儿恨不得把你吊起来往你身上滴蜡烛。
下面来说说为什么NOSE ELECTRON中的参数EKIN也会导致NaN杯具。EKIN是电子虚拟动能的热浴值。CPMD中的NOSE HOOVER热浴是扩展拉格朗日体系方法，是最严格的能量守恒的热浴方法。普通的热浴，要么直接Scale原子速度(Andersen)，要么加一个活塞(Beredsen)，原理过于简单，都是属于温度高了移走点能量，温度低了输入进去点能量，不可避免会造成能量恶意不守恒。看起来没问题，系综等概率假定已经失陷，统计热力学贞操已经失守。NOSE-HOOVER怎么搞呢？原子核的温度还是高了砍，低了补，只是增加或减少的能量都来自哈密顿里的一个储热罐(一个带质量，有动能，并乘以系数Q。原子核和这个虚拟物质通过Q传递热量)。这样就能保持总动能一定，无非那边多点这边少点的问题。具体到CPMD，则分别为原子核和电子搞了两个一定容量的储热罐在哈密顿里。电子的虚拟动能热浴值的设置很重要。原理和上面一样，设置过大或过小，都会改变绝热性。想避免不绝热问题，必须选择合适的EKIN。其实也简单，假如要跑一个100K的体系，那么让这个体系跑个初始温度100K不带热浴的NVE，看看EKIN的平均能量是多少。把这个值拿过来做NOSE ELECTRON的EKIN参数即可。至于频率值，大概设在20000-30000就行，不很严肃。
控温吗，无非设置了一个储热罐，然后需要设置目标温度，以及储热罐与体系之间交换能量的频率和幅度。
下面是一个输入的例子：
TEMPERATURE
100
NOSE IONS
100 1600
NOSE ELECTRONS
0.0003 10000
NOSE IONS第一个参数100等于The target temperature in Kelvin
第二个参数1600是耦合频率，thermostat frequency in cm-1, 即恒温器与体系多久交换一次能量，用波数作为单位，一般取体系特征红外振动频率，在水体系中那就是O-H键的振动频率，水有1600, 3600-3800 cm-1等几个峰。注意这个取值并不严格，Having a precise value for this frequency is not important, as one only wishes to insure that the thermostat will couple to the mode of interest.。
NOSE ELECTRONS第一个参数0.0003eV是电子动能，这个并非是预先可知的。
因此对于每一个体系都需要跑一个NVE测试一下，一般可跑一个500步的NVE，取后面200步的电子动能EKIN的平均值。NVE就是去掉所有的控温控压参数的CPMD。虽然精确值并不重要，但是这个值太小的话，意味着电子过冷，原子就会传递能量给电子，造成原子动能衰减，即 cause spurious damping of the ions；如果这个值太大，意味着电子的温度很高，动能太大，容易飞走，脱离B-O面(departures from the Born-Oppenheimer surface)。
NOSE ELECTRONS第二个参数10000是电子的耦合频率。它的取值要参考我前面帖子中的分析，即它最好是你整个体系中最高振动频率的2-3倍，保持CP绝热.注意这个最高振动频率不是O-H键的最高振动频率3700， 而是声子频率，即凝聚态体系由于晶格间振动耦合，会产生更高的频率。好在水这种体系的晶格间的振动耦合不明显。因此这里取10000-20000都差不多。
关于控温的幅度，即chain length，要通过NOSE PARAMETERS 关键词修改。由于默认值4对绝大多数体系是足够的，所以，除非有特别的非平衡模拟比如加温，否则不需要修改。
你的第二个问题超出我所理解的范围。水有磁性的吗？由于所有原子配对了，没有净的自旋，何来自旋密度？如果你设置了自旋激发态计算，是否设置正确？由于我在spin-polarization方面经验很少，所以欢迎虫友补充/。
下面是一个输入的例子：
TEMPERATURE
100
NOSE IONS
100 1600
NOSE ELECTRONS
0.0003 10000
NOSE IONS第一个参数100等于The target temperature in Kelvin
第二个参数1600是耦合频率，thermostat frequency in cm-1, 即恒温器与体系多久交换一次能量，用波数作为单位，一般取体系特征红外振动频率，在水体系中那就是O-H键的振动频率，水有1600, 3600-3800 cm-1等几个峰。注意这个取值并不严格，Having a precise value for this frequency is not important, as one only wishes to insure that the thermostat will couple to the mode of interest.。
NOSE ELECTRONS第一个参数0.0003eV是电子动能，这个并非是预先可知的。
因此对于每一个体系都需要跑一个NVE测试一下，一般可跑一个500步的NVE，取后面200步的电子动能EKIN的平均值。NVE就是去掉所有的控温控压参数的CPMD。虽然精确值并不重要，但是这个值太小的话，意味着电子过冷，原子就会传递能量给电子，造成原子动能衰减，即 cause spurious damping of the ions；如果这个值太大，意味着电子的温度很高，动能太大，容易飞走，脱离B-O面(departures from the Born-Oppenheimer surface)。
NOSE ELECTRONS第二个参数10000是电子的耦合频率。它的取值要参考我前面帖子中的分析，即它最好是你整个体系中最高振动频率的2-3倍，保持CP绝热.注意这个最高振动频率不是O-H键的最高振动频率3700， 而是声子频率，即凝聚态体系由于晶格间振动耦合，会产生更高的频率。好在水这种体系的晶格间的振动耦合不明显。因此这里取10000-20000都差不多。
关于控温的幅度，即chain length，要通过NOSE PARAMETERS 关键词修改。由于默认值4对绝大多数体系是足够的，所以，除非有特别的非平衡模拟比如加温，否则不需要修改。
你的第二个问题超出我所理解的范围。水有磁性的吗？由于所有原子配对了，没有净的自旋，何来自旋密度？如果你设置了自旋激发态计算，是否设置正确？由于我在spin-polarization方面经验很少，所以欢迎虫友补充/。
Monday, November 5, 2012
记者萨苏揭露日本侵华内幕，被日本政府控罪抄家
萨苏，《环球时报》签约记者，旅居日本多年，是一名普通的IT工程师。他利用业余时间收集抗日战争的日文史料，著有《国破山河在--从日本史料揭秘中国抗战》、《尊严不是无代价的--从日本史料揭秘中国抗战》、《与鬼为邻》等书，用翔实的资料和流畅的文笔深刻揭露了大量鲜为人知的日本侵华战争内幕。近日被日本政府控告以间谍罪和涉嫌骗取15万日圆失业救济，并以此为理由进入萨苏家搜查。
Subscribe to:
Posts (Atom)