Thursday, December 27, 2012

cpmd paper published in JACS journal


paper 1: 2012

Ab Initio Dynamics of Cellulose Pyrolysis: Nascent Decomposition Pathways at 327 and 600 °C

J. Am. Chem. Soc.2012134 (36), pp 14958–14972

We modeled nascent decomposition processes in cellulose pyrolysis at 327 and 600 °C using Car–Parrinello molecular dynamics (CPMD) simulations with rare events accelerated with the metadynamics method. We used a simulation cell comprised of two unit cells of cellulose Iβ periodically repeated in three dimensions to mimic the solid cellulose. To obtain initial conditions at reasonable densities, we extracted coordinates from larger classical NPT simulations at the target temperatures. CPMD-metadynamics implemented with various sets of collective variables, such as coordination numbers of the glycosidic oxygen, yielded a variety of chemical reactions such as depolymerization, fragmentation, ring opening, and ring contraction. These reactions yielded precursors to levoglucosan (LGA)—the major product of pyrolysis—and also to minor products such as 5-hydroxy-methylfurfural (HMF) and formic acid. At 327 °C, we found that depolymerization via ring contraction of the glucopyranose ring to the glucofuranose ring occurs with the lowest free-energy barrier (20 kcal/mol). We suggest that this process is key for formation of liquid intermediate cellulose, observed experimentally above 260 °C. At 600 °C, we found that a precursor to LGA (pre-LGA) forms with a free-energy barrier of 36 kcal/mol via an intermediate/transition state stabilized by anchimeric assistance and hydrogen bonding. Conformational freedom provided by expansion of the cellulose matrix at 600 °C was found to be crucial for formation of pre-LGA. We performed several comparison calculations to gauge the accuracy of CPMD-metadynamics barriers with respect to basis set and level of theory. We found that free-energy barriers at 600 °C are in the order pre-LGA < pre-HMF < formic acid, explaining why LGA is the kinetically favored product of fast cellulose pyrolysis.

Two temperatures are tested using CPMD+MTD. I am not sure about how they get the energy barriers. But it is interesting. Also, the CPMD details should be read carefully. I am curious how they set all the parameters too. For example, what about the couplings between NOF and EOF.

Friday, December 7, 2012

擅长的和喜欢的事情,转载加部分思考

不知道自己擅长什么吗?这没必要经过多么严密的考虑,只要看到身边的事能在心里有反应:“这件事我起来费不费劲?那件事呢?”。这样自问自答,答案就很明了了。

我擅长的事情就是我可以做的很好的东西。你平时对那类东西比较拿手 ,朋友们都喜欢让你帮他们出哪些主意,这都是不错的办法,还有你觉得什么对你来说特别简单,而别人做起来却没有你出色。

“比如我的中文阅读能力很强,英文就困难”。
“比如我喜欢数学公式,别人就未必。”


因此,比起“想做什么”,知道“自己擅长什么,适合什么”对你的人生更有帮助。

选择自己比较擅长的作为职业,自己也能相对地更加活跃。作为结果,自己能更出人头地,获得更多的薪水,这么多“更多”累计起来,也就能从这个工作中收获幸福。
最后,我想用创业天才,一手缔造了SGI和网景两大传奇公司的詹姆斯.克拉克的话作为结尾:
“我获得过最高评价之一,是皮克曼说的「Jim特别擅长于发现自己不擅长的东西」。如果谬赞属实的话,这无疑是我身上所有特质中对我帮助最大的。”

其他占了绝大部分的,平凡的,普通的人,还是不要把自己的兴趣作为职业比较好。
为什么呢?因为对平凡人来说,靠自己去发现自己喜欢什么,实际上是不太可能的。大部分的人作为普通职员,并不能靠自己创造出工作,只有极少数想乔布斯那样的人才能从事自己想做的工作。

像求职活动这样场合,对自己进行分析的意义也不大。所谓自我分析,大抵是为了找出自己想做什么。经过多番分析后得出的结果很可能是这样“我喜欢看到人们的笑容。所以去服务业吧~ 因为去和民可以看到人们的笑容,就定那了。”听起来似乎有点像黑色幽默,选择优衣库也是基于同样的逻辑。

Tuesday, December 4, 2012

questions about CPMD

1. When should I use NOSE MASSIVE?

For a simulation of a molecule in the gas phase the use of the MASSIVE thermostat is strongly recommended. In this context, massive does not refer to a more 'strict' or 'powerful' thermostat, but to a separate thermostat chain for each degree of freedom, i.e. massive referes to the total number of thermostats. This way a proper sampling of phase space is ensured, even if the various vibrational modes of the molecule(s) are only weakly coupled to the heat bath.

2. cpmd的几个问题。
   (1) symmetry 0 : ISOLATED system in a cubic/orthorhombic box with ISOLATED MOLECULE option activated.

cpmd _always_ uses periodic boundary conditions,
symmetry defines the shape of your supercell. SYMMETRY 0 is
a special case, where the effect of the periodic images ( In a box that is too small, a macromolecule may interact with its own image in a neighboring box, which is functionally equivalent to a molecule's "head" interacting with its own "tail". This produces highly unphysical dynamics in most macromolecules, ) is cancelled out, but for that to work correctly, there must not be any electron
density close to the borders of your supercell, so PBC or not
is no problem there.


           模拟结果和 
           symmetry 1: Simple CUBIC 
           or
           The ISOLATED MOLECULE keyword has only an effect on the calculation of the degrees of
freedom (3N-6 vs. 3N-3 for periodic systems).
           有何不同??
           center molecule
           The main purpose of this is to center the molecule (center
of mass) in the box. This is needed for the HOCKNEY Poisson solver. This solver gives wrong
results if the charge density is not centered in the computational box.
     (2)包括GAUSSIAN, 这些软件是如何放置分子的?在CPMD里,就是按Cartesian坐标放的吗?怎么保证他们在同一个盒子
            里?
     (3)使用不同的虚拟电子质量和时间间隔究竟有何不同?
      (4)虚拟电子和原子核(EOF AND NOF)的耦合,热原子到虚拟冷电子的耦合究竟和真实的电子原子核耦合是一回事情吗?后者的BOA被打破和虚拟冷电子大量吸收原子核能量是一回事情吗?? 

转载-模守恒赝势和超软赝势


 Norm-conserving 赝势是相当有名的而且是经彻底验证的。在这种方法中,赝波函数在定义的核心区域的截止半径以上是符合全电子波函数的。它要求改造后的波函数其在截止半径Rc之内的总电荷量仍要等于未改造前Rc之内总量的大小,这样赝势的精确度能够大幅的提升。因此,我们取距原子中心Rc 处为划分点,赝势产生示意图Rc 以上波函数完全一样保留,而 Rc以内则对波函数加以改造。主要是要把振荡剧烈的波函数改造成一个变化缓慢的波函数,而它需要是没有节点的。少了剧烈振荡不但允许只以相对很少的平面波来展开波函数,没有节点的(径向)波函数也意味着没有比它本征值更低的量子态来与它正交。求解内层电子的需要就自动消失了。我们以这样一个假的赝势能够在同样的本征值的情况下给出一价电子近似解,我们把它叫做是赝势 Vpseudo(Vp)。在 CASTEP 中引用的是最佳化的方法,然而描述第一列(碳,氮,氧)或过渡金属(镍,铜,钯)等局域化价电子轨域的所需之截止能量仍然经常还是太高。norm-conserving 赝势能够在实空间或是倒空间的波函数来使用;实空间的方法提供了对于系统而言比较好的可测量性。
超软赝势(ultrasoft pseudopotential)其特色是让波函数变得更平滑,也就是所需的平面波基底函数更少。Vanderbilt所提出来的超软赝势的想法是不用释放非收敛性条件,用这样的方法来产生更软的赝势。在这个方法里,虚波函数在核心范围是被允许作成尽可能越软(平滑),以致于截止能量可以被大大的减小。就技术上而言,这是靠着广义的正交条件来达成的。为了要重建整个总的电子密度,波函数平方所得到电荷密度必须在核心范围在加以附加额外的密度进去。这个电子云密度因此就被分成两各部分,第一部分是一个延伸是整各单位晶包平滑部分,第二部分是一个局域化在核心区域的自旋部分。前面所提的附加部分是只出现在密度,并不在波函数。这和像LAPW那样的方法不同,在那些方法中类似的方式是运用到波函数。
超软赝势(USP)产生算法保证了在预先选择的能量范围内会有良好的散射性质,这导致了赝势更好的转换性与精确性。超软赝势(USP)通常也借着把多套每个角动量通道当作价电子来处理浅的内层电子态。这也会使精确度跟转换性更加提升,虽然计算代价会比较高。目前,超软赝势(USP)只可以在倒置空间中使用。


计算平面波数目可取决於所含的元素

Norm-Conserving 型的 pseudopotential 从最早的 HSC 提出,经 BHS 建立一组涵盖整个週期表的参数,之后的KerkerTM,一直到 Optimised Pseudopotential,就在朝着

在兼顾准确性的情况下,尽可能使必须使用的平面波基底数目越少越好,它是直接影响(决定)所需计算量大小的量。一个贗势所需的基底数多少,可由 E_tot  E_cut 的收敛性本判定即平面波截止动能 E_cut 用到多大时则固态计算所求得的系统总能就不再改变。所需的 E_cut 越小也就是所谓的贗势越“软”。



大凡是贗势方法,碰到径向波函数上有节点的价电子态,由於贗化后波函数变得没节点,总是能得到相当软的贗势Heine  Cancelation Theorem),例如 n s-potential2s, 3s, ...)n p-potential (3p, 4p, ...)nd-potential (4d, 5d, ...);然而对於如 1s2p3d 这些首次出见的 nl 组合,其径向波函数都没有有节点,而 Norm-conservation 的条件又要求 Rc 内电荷量守恒,另外主量子数小的(或角动量子数大的)轨域其电子云都局域化得比较厉害。因此具有 2p  3d 价电子者,也就是所谓的第一横列(first row;硼、碳、氮、氧、氟)  3d 过渡金属元素,都是出了名的“硬”元素。

使用 Optimised  TM potential 虽然己经能够把的 pseudopotential 变得很 soft,但 Norm-conserving 条件对於原本就己经没有节点的价电子云分佈其改造及最佳化的程度,以现今日渐普遍的 ultrasoft 的贗势(它不必遵守 norm-conserving条件)来比,节省计算的程度仍是有限。总之,在使用 NCP 的情况下,计算量的大小是取决於原子(贗势)的种类这一点,仍是十分明显而普遍的观察。也就是说不同种类元素其 potential 的软硬的差异会令人明显感受到。

同值得一提的是,不设定 norm-convervation 而又要使贗势同等准确是需要付出一些程式发展的代价的,也因此,简易稳定的 norm-conserving pseudopotential 仍常被使用,不特别昂贵的计算、或为了与己经建立公信力的结果作比较测试,尤其是新发展的程式工具,往往是先实作在 NCP 架构上,之后才推广到 ultrasoft

在这节中要介绍的超软贗势,其特色是让波函数变得更平滑,也就是所需的平面波基底函数更少,在计算上是很大的好处。

超软贗势的 pseudo 波函数已经不必再遵守 Norm-Conserving 条件,它是靠定义Augmentation Charge(附加电荷)来达到所谓的 Generalised Norm-conserving 条件,基本上它是用来来做到把被砍掉之较局域化的电子云补回去,并使散射行为在参考态的能量附近仍保有不错的散射性质(变化率一致)。然而据(我们)实作的经验,其 log-derivative 图的涵盖度则不像由 norm-conserving K-B from 者那么好,也因此超软贗势需要引进了每个投影算符可以使用一个以上的参考能量的技术,就像是前面才介绍过的( Bloechl 所提出的)“广义完全可分离非局域贗势”,可以把用一个以上的参考本征态使用在同一个角动量的技术,如此可以弥补因为不遵守 Norm-Conserving 而造成的散射性质误差。此外,为了维持原贗核心内的价电子数而采用的附加电荷,会导致重叠算符 overlap operator S 出现在求期望值的公式上,这会使得 Kohn-Sham 方程式最后变成是一个广义本征值问题(Generalised Eigenvalue Problem),而不是一般的本征值问题。所幸,在数值演算法的运作上并没有太大的差异。

Vanderbilt 所提出的超软贗势,由於捨弃了 norm-conserving 条件,所以让原本即无节点的波函数也能被改造得非常平滑,这可以让我们使用很少的平面波基底展开,也因此计算代价比使用 norm-conserving 型的贗势更划算。

Sunday, November 18, 2012

can cpmd run for chemical reaction?

The answer is yes. The following information is from CPMD list by Prof. Hutter.


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

动力学模拟: 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,反复强调你的创新性(一定要“反复”,因为审稿人会忽 视),一般没有什么问题。另外,因为审稿人是带着寻找问题的模式去评判文章的,所以在正文中的每一句话不要过度发散,否则很容易招致不严谨或者补充数据的 评语。