[QE-users] NVE simulation in cp.x - large fluctuation in temperature

[Stefano de Gironcoli]

Dear Jie Peng,

  suppose you were running a model harmonic system in 1 dimension.

   M a = - K x

  at fixed energy  E.

  The kinetic energy would fluctuate harmonically between 0 (at 
maximum/minimum elongation) and E at the equilibrium distance .

  On average the Kinetic energy would be E/2 and its fluctuation some 
big fraction of E^2

  Something like sigma^2 = 1/T \int_0^T (E cos^2(2pi t/T) -E/2)^2 dt = 
E^2 1/T \int_0^T (cos(4pi t/T)/2)^2 dt = (E/2)^2 1/2pi \int_0^2pi 
cos^2(x) dx = (E/2)^2 / 2

   or

   sigma =  1/sqrt(2) * E/2 = 1/sqrt(2) avg EKin

   with 1 degree of freedom the mean square fluctuation of the kinetic 
energy is 70% of its average !

   you have 3 atoms in your cell hence 9 degrees of freedom. Assuming 
each contributes independently to the average this goes down by a factor 
1/sqrt(9)=1/3

   actually more likely just 1/sqrt(6) as the total momentum is 
conserved so only 6 modes at Gamma are actually excited...

   If you perform your simulation in a bigger supercell with more atoms 
(more degrees of freedom) the average will be more stable     ( 
proportionally to  1/sqrt(#deg.of.freedom-3 )  ... moreover the thermal 
excitations of vibrational modes will be sampled more faithfully.

   best

stefano


On 13/04/2018 21:39, Jie Peng wrote:
> Dear all
>
> I have been running MD simulations on HfS2 using cp.x code in Quantum 
> espresso. I start from initial configuration obtained from pwscf 
> vc-relax, and relax the system using cp.x by consecutive steps of: 
> electron relaxation->ionic relaxation->cell relaxation. Then, I just 
> directly start a NVE simulation starting from the equilibrium 
> configuration. I expect the system to almost stay stationary or the 
> temperature should be very small since I am allowing dynamics in a 
> system that is already in equilibrium. However, what I see is a huge 
> fluctuation in the /tmpp/ output of cp.x, as I attach a figure showing 
> variation of tmpp (Ionic temperature) with simulation time
>
> I did this because it is suggested in the user guide you should apply 
> an initial displacement to the atoms in your system after the 
> relaxation since otherwise there will not be any dynamics. But what I 
> see here is a large fluctuation of the system temperature.
>
> The thinking or questions here are
>
> 1.Does the tmpp represents the physical temperature of the system 
> here? I think it should be since it is the temperature corresponding 
> to kinetic energy of the ions.
>
> 2.It above point is true, why is the temperature varying so fiercely? 
> Am I setting incorrect parameters, for instance the timestep or the 
> fictitious mass? But I took those from previous simulation steps where 
> I did the relaxation, and they all worked well since they successfully 
> drived my system to equilibrium, satisfying the convergence threshold 
> on total energy, forces acting on atoms, and the fictitious electron 
> kinetic energy. I am confused at this point.
>
> The input file for NVE simulation is attached here:
>
> /&control/
> /    calculation='cp',/
> /    title='Halfnium disulfide'/
> /    restart_mode='restart',/
> /    ndr=53,/
> /    ndw=54,/
> /    nstep=50000,/
> /    iprint=10/
> /    isave=100,/
> /    tstress = .true./
> /    tprnfor = .true./
> /    dt=10,/
> /    wf_collect=.true./
> /    etot_conv_thr=1e-6/
> /    forc_conv_thr=1e-3/
> /    ekin_conv_thr=1e-5/
> /    prefix='HfS2',/
> /    pseudo_dir='/home/jpeng/HfS2/potential'/
> /    outdir='./tmp/',/
> / //
> / &system/
> /    ibrav= 4,/
> /    a=3.6529/
> /   c=5.6544/
> /    nat=  3, ntyp= 2,/
> /    ecutwfc =50/
> /    vdw_corr='DFT-D',/
> / !   lspinorb=.true./
> / !   noncolin=.true./
> / !   ecutrho=300/
> / !   nbnd=14/
> /!    occupations='smearing'/
> /!    smearing='gaussian'/
> /!   degauss=0.01/
> / !  nspin=2/
> / !   starting_magnetization(1)=0.1/
> /! Hf  95.94  Hf.pbe-mt_fhi.UPF/
> /! S  32.065  S.pbe-mt_fhi.UPF/
> ///
> / &electrons/
> /    electron_dynamics='verlet'/
> /    electron_velocities='zero'/
> /    emass=400/
> /    emass_cutoff=1/
> ///
> / &ions/
> /    ion_dynamics = 'verlet'/
> /    ion_damping=0.1/
> /!    ion_nstepe=10/
> / //
> / &cell/
> /    cell_dynamics = 'none'/
> //
> ///
> /ATOMIC_SPECIES/
> / Hf  95.94  Hf.pbe-mt_fhi.UPF/
> / S  32.065  S.pbe-mt_fhi.UPF/
> /ATOMIC_POSITIONS (crystal)/
> /Hf      -0.000000000  -0.000000000  -0.000000000/
> /S        0.666666667   0.333333333   0.257234636/
> /S        0.333333333   0.666666667  -0.257234636/
>
> Anyone could help me on it? Thank you very much.
>
> Best
> Jie
> -- 
> ------------------------------------------------------------------------------------------------------------------------
> Jie Peng
> PhD student
> 2134 Glenn Martin Hall, Mechanical Engineering, University of Maryland
> College Park, Maryland, USA
> Phone:(+1) 240-495-9445
> Email: jiepeng at umd.edu <mailto:jiepeng at umd.edu>
>
>
>
> _______________________________________________
> users mailing list
> users at lists.quantum-espresso.org
> https://lists.quantum-espresso.org/mailman/listinfo/users

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20180414/87526b84/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 56649 bytes
Desc: not available
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20180414/87526b84/attachment.png>