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

Stefano de Gironcoli degironc at sissa.it
Tue Apr 17 16:08:24 CEST 2018


Dear Jie Peng,

    one way to set the temperature of your system in a MD simulation is 
by the use of the

    ion_temperature and tempw variable pair in the &IONS namelist

    read the description of these variable in the input documentation ( 
http://www.quantum-espresso.org/Doc/INPUT_PW.htm )

    one thing that you can do is to start from, or close to, the 
equilibrium geometry and define ion_temperature = 'initial' setting 
tempw to twice the desired temperature.

    the system will start at the potential energy minumum. after a 
little while the initial kinetic energy will redistribute among all the 
degrees of freedom (kinetic and potential-related ones) and the kinetic 
energy should oscillate roughly around the desired temperature.

    If you want to further modify the system temperature you can then 
play with other options of the ion_temperature variable.

     best,

stefano


On 14/04/2018 20:16, Jie Peng wrote:
> Stefano:
>
> Thanks! This is so informative, and indeed I tried larger unit cells. 
> But then I have another issue which is: the temperature never reach 
> desired one. Below is what I did
>
> For a 2*2*2 supercell with input file:
>
> /&control/
> /    calculation='vc-cp',/
> /    title='Halfnium disulfide'/
> /    restart_mode='restart',/
> /    ndr=52,/
> /    ndw=53,/
> /    nstep=1200000,/
> /    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=7.3058/
> /   c=11.3088/
> /    nat=24, 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/
> /    tempw=300/
> /    fnosep=6/
> /!    ion_nstepe=10/
> / //
> / &cell/
> /    cell_dynamics = 'pr'/
> //
> ///
> /ATOMIC_SPECIES/
> / Hf  95.94  Hf.pbe-mt_fhi.UPF/
> / S  32.065  S.pbe-mt_fhi.UPF/
> /ATOMIC_POSITIONS angstrom/
> /   Hf    0.0000000000    0.0000000000    0.0000000000/
> /   S    1.8264500020    1.0545013980    1.4545075260/
> /   S   -0.0000000020    2.1090027990   -1.4545075260/
> /   Hf    0.0000000000    0.0000000000    5.6544000000/
> /   S    1.8264500020    1.0545013980    7.1089075260/
> /   S   -0.0000000020    2.1090027990    4.1998924740/
> /   Hf   -1.8264500000    3.1635041980    0.0000000000/
> /   S    0.0000000020    4.2180055960    1.4545075260/
> /   S   -1.8264500020    5.2725069970   -1.4545075260/
> /   Hf   -1.8264500000    3.1635041980    5.6544000000/
> /   S    0.0000000020    4.2180055960    7.1089075260/
> /   S   -1.8264500020    5.2725069970    4.1998924740/
> /   Hf    3.6529000000    0.0000000000    0.0000000000/
> /   S    5.4793500020    1.0545013980    1.4545075260/
> /   S    3.6528999980    2.1090027990   -1.4545075260/
> /   Hf    3.6529000000    0.0000000000    5.6544000000/
> /   S    5.4793500020    1.0545013980    7.1089075260/
> /   S    3.6528999980    2.1090027990    4.1998924740/
> /   Hf    1.8264500000    3.1635041980    0.0000000000/
> /   S    3.6529000020    4.2180055960    1.4545075260/
> /   S    1.8264499980    5.2725069970   -1.4545075260/
> /   Hf    1.8264500000    3.1635041980    5.6544000000/
> /   S    3.6529000020    4.2180055960    7.1089075260/
> /   S    1.8264499980    5.2725069970    4.1998924740 /
>
> Now the temperature is still oscillating around 80-90K after 50 ps 
> simulation as if that is the desired temperature
>
>
> And for a 3*3*3 supercell:
>
> /&control/
> /    calculation='vc-cp',/
> /    title='Halfnium disulfide'/
> /    restart_mode='restart',/
> /    ndr=52,/
> /    ndw=53,/
> /    nstep=1200000,/
> /    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=10.9587/
> /   c=16.9632/
> /    nat=81, 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/
> /    tempw=300/
> /    fnosep=6/
> /!    ion_nstepe=10/
> / //
> / &cell/
> /    cell_dynamics = 'pr'/
> //
> ///
> /ATOMIC_SPECIES/
> / Hf  95.94  Hf.pbe-mt_fhi.UPF/
> / S  32.065  S.pbe-mt_fhi.UPF/
> /ATOMIC_POSITIONS angstrom/
> /Hf 0.0000000000 0.0000000000 0.0000000000/
> /S 1.8264500020 1.0545013980 1.4545075260/
> /S -0.0000000020 2.1090027990 -1.4545075260/
> /Hf 0.0000000000 0.0000000000 5.6544000000/
> /S 1.8264500020 1.0545013980 7.1089075260/
> /S -0.0000000020 2.1090027990 4.1998924740/
> /Hf 0.0000000000 0.0000000000 11.3088000000/
> /S 1.8264500020 1.0545013980 12.7633075260/
> /S -0.0000000020 2.1090027990 9.8542924740/
> /Hf -1.8264500000 3.1635041980 0.0000000000/
> /S 0.0000000020 4.2180055960 1.4545075260/
> /S -1.8264500020 5.2725069970 -1.4545075260/
> /Hf -1.8264500000 3.1635041980 5.6544000000/
> /S 0.0000000020 4.2180055960 7.1089075260/
> /S -1.8264500020 5.2725069970 4.1998924740/
> /Hf -1.8264500000 3.1635041980 11.3088000000/
> /S 0.0000000020 4.2180055960 12.7633075260/
> /S -1.8264500020 5.2725069970 9.8542924740/
> /Hf -3.6529000000 6.3270083960 0.0000000000/
> /S -1.8264499980 7.3815097940 1.4545075260/
> /S -3.6529000020 8.4360111950 -1.4545075260/
> /Hf -3.6529000000 6.3270083960 5.6544000000/
> /S -1.8264499980 7.3815097940 7.1089075260/
> /S -3.6529000020 8.4360111950 4.1998924740/
> /Hf -3.6529000000 6.3270083960 11.3088000000/
> /S -1.8264499980 7.3815097940 12.7633075260/
> /S -3.6529000020 8.4360111950 9.8542924740/
> /Hf 3.6529000000 0.0000000000 0.0000000000/
> /S 5.4793500020 1.0545013980 1.4545075260/
> /S 3.6528999980 2.1090027990 -1.4545075260/
> /Hf 3.6529000000 0.0000000000 5.6544000000/
> /S 5.4793500020 1.0545013980 7.1089075260/
> /S 3.6528999980 2.1090027990 4.1998924740/
> /Hf 3.6529000000 0.0000000000 11.3088000000/
> /S 5.4793500020 1.0545013980 12.7633075260/
> /S 3.6528999980 2.1090027990 9.8542924740/
> /Hf 1.8264500000 3.1635041980 0.0000000000/
> /S 3.6529000020 4.2180055960 1.4545075260/
> /S 1.8264499980 5.2725069970 -1.4545075260/
> /Hf 1.8264500000 3.1635041980 5.6544000000/
> /S 3.6529000020 4.2180055960 7.1089075260/
> /S 1.8264499980 5.2725069970 4.1998924740/
> /Hf 1.8264500000 3.1635041980 11.3088000000/
> /S 3.6529000020 4.2180055960 12.7633075260/
> /S 1.8264499980 5.2725069970 9.8542924740/
> /Hf 0.0000000000 6.3270083960 0.0000000000/
> /S 1.8264500020 7.3815097940 1.4545075260/
> /S -0.0000000020 8.4360111950 -1.4545075260/
> /Hf 0.0000000000 6.3270083960 5.6544000000/
> /S 1.8264500020 7.3815097940 7.1089075260/
> /S -0.0000000020 8.4360111950 4.1998924740/
> /Hf 0.0000000000 6.3270083960 11.3088000000/
> /S 1.8264500020 7.3815097940 12.7633075260/
> /S -0.0000000020 8.4360111950 9.8542924740/
> /Hf 7.3058000000 0.0000000000 0.0000000000/
> /S 9.1322500020 1.0545013980 1.4545075260/
> /S 7.3057999980 2.1090027990 -1.4545075260/
> /Hf 7.3058000000 0.0000000000 5.6544000000/
> /S 9.1322500020 1.0545013980 7.1089075260/
> /S 7.3057999980 2.1090027990 4.1998924740/
> /Hf 7.3058000000 0.0000000000 11.3088000000/
> /S 9.1322500020 1.0545013980 12.7633075260/
> /S 7.3057999980 2.1090027990 9.8542924740/
> /Hf 5.4793500000 3.1635041980 0.0000000000/
> /S 7.3058000020 4.2180055960 1.4545075260/
> /S 5.4793499980 5.2725069970 -1.4545075260/
> /Hf 5.4793500000 3.1635041980 5.6544000000/
> /S 7.3058000020 4.2180055960 7.1089075260/
> /S 5.4793499980 5.2725069970 4.1998924740/
> /Hf 5.4793500000 3.1635041980 11.3088000000/
> /S 7.3058000020 4.2180055960 12.7633075260/
> /S 5.4793499980 5.2725069970 9.8542924740/
> /Hf 3.6529000000 6.3270083960 0.0000000000/
> /S 5.4793500020 7.3815097940 1.4545075260/
> /S 3.6528999980 8.4360111950 -1.4545075260/
> /Hf 3.6529000000 6.3270083960 5.6544000000/
> /S 5.4793500020 7.3815097940 7.1089075260/
> /S 3.6528999980 8.4360111950 4.1998924740/
> /Hf 3.6529000000 6.3270083960 11.3088000000/
> /S 5.4793500020 7.3815097940 12.7633075260/
> /S 3.6528999980 8.4360111950 9.8542924740/
>
> The temperature never goes up to even 5K!
>
>
>
> And FYI, in the post I sent couple of days ago, I have already checked 
> that cp relaxation calculation on 2*2*2 and 3*3*3 supercell gives me 
> an accurate description of the lattice structure.
>
> So below I write down some of my thought that lead to some questions:
>
> 1. What determines the time spent on heating up my system? I can tell 
> from the plots that for 2*2*2 it took maybe 2-3 ps to heat the system 
> to around 80-90K while for 3*3*3 it took 0.5ps to heat it to around 
> 1K. But I am not sure which parameter controls the rate at which the 
> temperature is being raised. I know it is difficult to determine, 
> since your system size, thermal properties and the thermostat all come 
> into play. But do you have a reference for me to check on this topic 
> or according to your experience, what parameters dominate the process 
> of raising temperature of a system.
>
> 2. Instead of the desired 300K, both my systems stabilize at a 
> temperature much lower. There seems to be a "barrier" stopping 
> increase in temperature. Right now I do not have a clue on what is 
> going on here, intuitive thinking would be the thermostat is not able 
> to transfer kinetic energy to my system. The thermostat has a 
> /fnosep=6 /which is about in the middle of the phonon spectrum of my 
> system (around 0 to 11Thz) as suggested by the users guide. A possible 
> explanation would be that the thermostat is injecting energy into the 
> system too slowly that I can hardly see a temperature increase, while 
> only on a nanosecond scale or even larger scale will it become prominent.
>
> Thank you very much! I am so grateful for your explanation provided!
>
> Best
> Jie
>
> On Sat, Apr 14, 2018 at 11:20 AM, Stefano de Gironcoli 
> <degironc at sissa.it <mailto:degironc at sissa.it>> wrote:
>
>     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
>>     <mailto:users at lists.quantum-espresso.org>
>>     https://lists.quantum-espresso.org/mailman/listinfo/users
>>     <https://lists.quantum-espresso.org/mailman/listinfo/users>
>
>
>
>
> -- 
> ------------------------------------------------------------------------------------------------------------------------
> 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>
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20180417/578d163d/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 38406 bytes
Desc: not available
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20180417/578d163d/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 43876 bytes
Desc: not available
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20180417/578d163d/attachment-0001.png>
-------------- 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/20180417/578d163d/attachment-0002.png>


More information about the users mailing list