【原创】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方面经验很少，所以欢迎虫友补充/。
No comments:
Post a Comment