[QE-users] 答复: Question on CPMD

Fri Jan 18 16:51:38 CET 2019

Dear Liang,

yes, CPMD can run in NVT and also in NPH, NPT. But I would solve your problem in NVE first, and then play with NVT.

In any case from your plots and inputs it is not hard to understand what's wrong in your calculations.
Please also refer to this paper

Physical Review A 44 (10), 6334, 1991

in order to learn how the CP method works and what you should take care of.
I try to explain you here quickly the meaning of the symbols in the file .evp you mentioned earlier, and the parameters to be tuned in your input.

First of all you should know that:

- Ekinc = kinetic energy of the coefficients of the electronic wavefunctions
- Etot = total energy of the electrons (i.e. potential on which the ions move)
- Econs = Etot + Ekin_ions = total energy of the whole real system (ions+electrons)
- Econt = Etot + Ekin_ions + Ekinc = total energy of the whole fake system (i.e. real system + electronic wavefunctions)
 (This holds in a NVE simulation. For NVT you would also have the contribution of the thermostat.)
  As Riccardo pointed out, Econt must be constant. Econs should oscillate around a constant.

In the CP method you are propagating not only the ions, but also the electronic wavefunctions, that have a mass (emass) and a kinetic energy (ekinc). The mass is arbitrary - you decide. But if the mass is too big, the kinetic energy of the electronic wavefunctions becomes too close to the one of the ions and the two systems will exchange energy, so that the kinetic energy of the electronic wavefunctions (ekinc) will grow again and again. This is exactly what is happening in your calculations: you have ekinc continuously increasing. In a NVE simulation this means that T_ions (i.e. E_kin_ions) decreases, whereas this cannot happen in a NVT simulation, because you have a thermostat that keeps T_ions fixed.
In any case, if E_kinc increases, E_cont increases, and your calculations are wrong as your constant of motion is not conserved.

I see that your E_kin_ion (= E_cons - E_tot) is ~0.1eV and your starting ekinc is ~0.2eV !!! You have definitely a too high ekinc, and you should decrease emass. Go to emass = 300 au or even less, and check that ekinc < Ekin_ions. In general
Ekinc < 1/10 Ekin_ions should be enough. Then you should take care that your timestep is not too high for this small mass, and possibly choose a shorter dt. A too large dt would cause again Econt to drift during your simulation.

In summary:
- read the paper above, or any good paper about the CP method.
- run NVE first. Keep in mind that Econs - Etot = Ekin_ions. In NVE you also have that Econt = Econs + ekinc.
- try several emass and dt, looking each time at:
 a) the comparison between ekinc and Ekin_ions (ekinc must be at least 1/10 of Ekin_ions)
 b) the conservation of Econt (and the oscillation of Econs around a constant): if Econt still grows, then choose a shorter timestep. Your plots should look like Riccardo's.

Finally, please use a stricter threshold for ekinc: ekinc_conv_thr = 1d-4 is too high, I would rather go to 1d-9.

Good luck!

  Giuliana

-- **************************************************** Giuliana Materzanini Theory and Simulation of Materials École Polytechnique Fédérale de Lausanne ****************************************************

On Jan 18, 2019 09:01, LEUNG Clarence <liangxy123 at hotmail.com> wrote:
Dear Riccardo Bertossa,

Can CPMD run in NVT? Now I perform a system at 300K.
[cid:699ebdcb-eca4-4293-bb11-3ec1fa1e40db]

Ekinc is go up all the time:
[cid:8e608b01-ff5b-495c-854f-8a5da861e9c3]
ETOT and Econs go up all the time:
[cid:4cc11732-f5e7-4524-965a-575db70e4f00]
Econt is not a constant:
[cid:f77d6356-1df6-473c-a91a-0c0da4b153a0]

I have run 0K simulation to go to the ground state and give 300K, the input file is

&CONTROL
  calculation='cp' ,
  etot_conv_thr = 3.5D-6 ,
  forc_conv_thr = 4.0D-4 ,
  ekin_conv_thr = 1d-4 ,
  pseudo_dir='/home/qeuser/SSSP_acc_PBE' ,
  nstep = 70000 ,
  tstress = .false. ,
  tprnfor = .false. ,
  dt=2.067,
  isave = 100,
  iprint = 10,
  ndr = 56,
  ndw = 57,
  restart_mode = 'restart',
  verbosity = 'low' ,
/

&SYSTEM
  ibrav=14,
  celldm(1)=36.1187133640d0, celldm(2)=0.6119121863d0, celldm(3)=1.3079965678d0,
  celldm(4)=0.0000000000d0, celldm(5)=0.0000000000d0, celldm(6)=0.0000000000d0,
  nat=63,
  ntyp=2,
  ecutwfc=50,
  ecutrho=400,
  input_dft='PBE',
  vdw_corr = 'DFT-D3' ,
  nr1b = 16 ,
  nr2b = 16,
  nr3b = 16 ,
  nosym = .true. ,
/
&ELECTRONS
  emass = 400
  emass_cutoff = 2.50,
  electron_dynamics = 'verlet' ,
/
&IONS
  ion_dynamics  = 'verlet',
  tempw = 300 ,
  ion_temperature = 'nose' ,
  fnosep=6.6666,
  ion_radius(1)=1.0,
  ion_radius(2)=1.0,
/

ATOMIC_SPECIES
  C 12.010700d0 C_pbe_v1.2.uspp.F.UPF
  P  30.973800d0 P.pbe-n-rrkjus_psl.1.0.0.UPF

ATOMIC_POSITIONS {crystal}
P        0.207062953   0.994563394   0.532135751
P        0.042823633   0.994796583   0.465770116
C        0.157559602   0.119473623   0.509211511
C        0.092142096   0.120398444   0.486938482
P        0.453113260   0.016237452   0.533698293
P        0.289980002   0.986689982   0.464529400
C        0.399646081   0.133053286   0.505994952
C        0.335858066   0.117232823   0.482213070
P        0.704320783   0.007706925   0.532982064
P        0.538824580   0.008408401   0.468045855
C        0.654037234   0.132822314   0.511474768
C        0.587848371   0.135094109   0.489114957
P        0.956731588  -0.000006745   0.532145072
P        0.791554157   0.002228032   0.466394516
C        0.907495292   0.125510515   0.509696521
C        0.841931621   0.126642797   0.487430517
P        0.208266817   0.243782922   0.529541853
P        0.043965854   0.246589308   0.464392818
C        0.158728245   0.369838186   0.509463992
C        0.093217169   0.370873210   0.486709312
P        0.424177849   0.286343852   0.516852430
P        0.290089532   0.244400471   0.459582516
C        0.349540011   0.348792825   0.487513040
P        0.705574516   0.255089554   0.532107492
P        0.542012221   0.260506971   0.464249905
C        0.656129286   0.380833823   0.512136142
C        0.590214180   0.381152001   0.490082742
P        0.957108151   0.250073183   0.531244245
P        0.792998637   0.251922321   0.464998669
C        0.907321789   0.375462368   0.510050868
C        0.841746979   0.376334296   0.487596817
P        0.208402877   0.492352243   0.532960869
P        0.043963546   0.496372026   0.465432581
C        0.160144417   0.618157610   0.509696507
C        0.094342113   0.619600026   0.487659242
P        0.471073678   0.513328655   0.548894023
P        0.309129724   0.486227608   0.474070256
C        0.414925399   0.629235003   0.523346545
C        0.350961159   0.616373438   0.499640968
P        0.705395802   0.505478539   0.534378721
P        0.537192856   0.506995560   0.474222102
C        0.656825419   0.631333460   0.511436510
C        0.590845794   0.632466501   0.490389008
P        0.956482750   0.499986748   0.532254831
P        0.791270143   0.501641183   0.466922839
C        0.906923580   0.625196225   0.510129362
C        0.841214140   0.626345388   0.488309739
P        0.207667541   0.743674910   0.533472131
P        0.043998636   0.744658734   0.465730512
C        0.158052412   0.868776604   0.510578825
C        0.092761413   0.869723388   0.488195850
P        0.461194613   0.761061943   0.539728957
P        0.299530254   0.737206731   0.471960596
C        0.407654987   0.881987587   0.515312903
C        0.343486860   0.870097791   0.492007540
P        0.706175574   0.757137135   0.532107627
P        0.539960743   0.756974729   0.469099032
C        0.655720285   0.882343139   0.510585276
C        0.589586421   0.883272199   0.489712695
P        0.956945276   0.749651294   0.531828394
P        0.791991567   0.752132650   0.465926241
C        0.907222485   0.875386768   0.509910655
C        0.841486983   0.876756958   0.488249976

K_POINTS {gamma}

Thanks very much.

LIANG Xiongyi
________________________________
发件人: users <users-bounces at lists.quantum-espresso.org> 代表 Riccardo Bertossa <rbertoss at sissa.it>
发送时间: 2019年1月18日 2:24
收件人: users at lists.quantum-espresso.org
主题: Re: [QE-users] Question on CPMD

dear LIANG,

when you run in the NVE you should have something like this:

[cid:part1.2BBB1DF7.418C4223 at sissa.it]

[cid:part2.309ED1CB.07CFEBCF at sissa.it]

if I am not wrong:

ekinc is the kinetic energy of the electrons

etot is the total potential energy (calculated with DFT)

econs is the energy that to make a good simulation you want to make as constant as possible, but it will oscillate because electrons have a mass and a kinetic energy, so it will never be constant.

econt is the constant of motion of the CP lagrangian, and this must always remain approximately constant, within the integration error of the verlet algorithm. If it is not constant your simulation is not good and you have to change some parameters.

best regards,

Riccardo Bertossa

SISSA

On 17/01/19 15:54, LEUNG Clarence wrote:
Dear QE users,

From CPMD, we can get    EKINC  ETOT   ENTHAL ECONS  and ECONT of the system. Which one is the total energy of the system?
ETOT? or Ekinc + Etot ? or Etot - Ekinc ?

Thanks.

LIANG Xiongyi
City University of Hong Kong

_______________________________________________
users mailing list
users at lists.quantum-espresso.org<mailto:users at lists.quantum-espresso.org>
https://lists.quantum-espresso.org/mailman/listinfo/users