[Pw_forum] finite displacement phonon calculation

Merlin Meheut merlin.meheut at gmail.com
Tue Apr 26 17:54:53 CEST 2016


Dear QE users and developpers,

I would like to compute phonon properties using the finite displacement
method, as proposed in the example of quantum espresso version 5.1
(PHonon/FD/example). For frequencies at q=0, the answer is reasonably close
to the results of DFPT, but for q-vectors different from zero, I find very
negative frequencies, and very different from DFPT.

Follows a description of my calculations. Any help would be thoroughly
appreciated.

My system is forsterite Mg2SiO4 (orthorhombic Pbnm):

&system
    ibrav =0, celldm(1)=8.985459,
    nat =28, ntyp = 3, ecutwfc =${a}.0, ecutrho = ${b}.0
/&end

CELL_PARAMETERS (alat=  8.98545900)
   1.011800931   0.000000000   0.000000000
   0.000000000   2.175636534   0.000000000
   0.000000000   0.000000000   1.273688789

ATOMIC_POSITIONS (crystal)
Mg      -0.000000000  -0.000000000  -0.000000000
Mg       0.500000000   0.500000000  -0.000000000
Mg      -0.000000000   0.000000000   0.500000000
Mg       0.500000000   0.500000000   0.500000000
Mg       0.991704573   0.277641327   0.250000000
Mg       0.491704573   0.222358673   0.750000000
Mg       0.008295427   0.722358673   0.750000000
Mg       0.508295427   0.777641327   0.250000000
Si       0.425939676   0.093656753   0.250000000
Si       0.925939676   0.406343247   0.750000000
Si       0.574060324   0.906343247   0.750000000
Si       0.074060324   0.593656753   0.250000000
O        0.765422396   0.091530172   0.250000000
O        0.265422396   0.408469828   0.750000000
O        0.234577604   0.908469828   0.750000000
O        0.734577604   0.591530172   0.250000000
O        0.222541031   0.447081956   0.250000000
O        0.722541031   0.052918044   0.750000000
O        0.777458969   0.552918044   0.750000000
O        0.277458969   0.947081956   0.250000000
O        0.277386450   0.162817013   0.033263088
O        0.777386450   0.337182987   0.966736912
O        0.722613550   0.837182987   0.533263088
O        0.222613550   0.662817013   0.466736912
O        0.722613550   0.837182987   0.966736912
O        0.222613550   0.662817013   0.033263088
O        0.277386450   0.162817013   0.466736912
O        0.777386450   0.337182987   0.533263088

On one had I have made a phonon (perturbation theory) calculation on a
2x2x2 grid of the Brillouin zone, deliberately without Born charges and
dielectric tensor (since I do not do it with finite displacements):

&inputph
   amass(1)=23.9850,
   amass(2)=27.9769,
   amass(3)=15.9949,
   alpha_mix(1) = 0.7,
   ldisp=.true., nq1=2, nq2=2, nq3=2,
   tr2_ph =  1.0D-16,
   prefix='FOR',
   fildyn='mat.$PREFIX',
   lraman=.false.,
   epsil =.false.,
   trans =.true.,
   zue = .false.,
   outdir='*/tmpdir/*$LOGNAME/'
/&end
0.0 0.0 0.0

Then I obtained an ifc matrix with q2r and re-calculated the frequencies
with matdyn on the 8 q-points (finding the same frequencies as directly
from the individual dynamical matrices).
For example, for q= (0.5,0,0) I find the following frequencies:

  152.7355  152.7355  159.2373  159.2373  159.6323  159.6323
  171.3805  171.3805  177.2753  177.2753  191.6439  191.6439
  202.0364  202.0364  244.3864  244.3864  244.9694  244.9694
  264.4823  264.4823  277.6309  277.6309  279.3910  279.3910
  287.4779  287.4779  288.9345  288.9345  295.9524  295.9524
  307.4176  307.4176  309.7795  309.7795  329.0229  329.0229
  332.7833  332.7833  339.4032  339.4032  346.7481  346.7481
  371.5853  371.5853  403.7708  403.7708  423.2761  423.2761
  436.9743  436.9743  458.7147  458.7147  459.0536  459.0536
  482.4561  482.4561  493.6059  493.6059  501.4041  501.4041
  530.9486  530.9486  550.9678  550.9678  560.3782  560.3782
  603.1617  603.1617  786.0294  786.0294  800.2429  800.2429
  820.2436  820.2436  842.1800  842.1800  846.5773  846.5773
  854.6477  854.6477  947.1591  947.1591  992.4286  992.4286

On the other hand, I realized a finite displacement calculation with fd.in
below:

&inputfd
 fd_prefix      = 'FOR'
 fd_outdir      = '*/tmpdir/mmeheut/*'
 fd_outfile     = 'displaced'
 fd_outfile_dir = './fd_files'
 nrx1           = 2
 nrx2           = 2
 nrx3           = 2
 innx           = 1
 de             = 0.01
/

and fd_ifc.in

&input
  prefix='FOR'
  nrx1=2,
  nrx2=2,
  nrx3=2,
  de=0.01,
  file_force='./fd_files/force',
  file_out='./For-rel-pbeUS.80Ryx4-222111.fd222',
  innx=1
  nodispsym=.false.
  noatsym=.false.
  verbose=.true.
  hex=.false.
/

the ifc matrix I obtain is a bit different from DFPT but it is difficult to
make clearer statement. When I use matdyn to obtain frequencies, at q=0 it
is correct, but at q= (0.5,0,0) for example, I have:

 -549.4798 -549.4798 -542.0251 -542.0251 -270.8103 -270.8103
 -270.3831 -270.3831 -268.2961 -268.2961 -255.8122 -255.8122
 -229.8735 -229.8735 -224.8183 -224.8183 -219.7591 -219.7591
 -205.2850 -205.2850 -200.4448 -200.4448 -197.4146 -197.4146
  131.2067  131.2067  143.8600  143.8600  155.3084  155.3084
  163.2342  163.2342  181.7740  181.7740  200.2658  200.2658
  213.1614  213.1614  243.0599  243.0599  254.4463  254.4463
  257.8945  257.8945  270.8913  270.8913  271.1355  271.1355
  277.1344  277.1344  280.2588  280.2588  296.9960  296.9960
  339.1382  339.1382  344.2783  344.2783  346.3596  346.3596
  360.5026  360.5026  369.8420  369.8420  387.3814  387.3814
  457.0229  457.0229  491.8080  491.8080  512.6835  512.6835
  576.2465  576.2465  588.0433  588.0433  596.5263  596.5263
  608.6300  608.6300  623.7900  623.7900  628.1992  628.1992

The symetry does seem correct, but the frequencies are completely off..

Thanks again for any help,

Sincerely,

 Merlin Méheut

-- 
Merlin Méheut, Géosciences et Environnement Toulouse,
OMP, 14 avenue Edouard Belin, 31400 Toulouse, France

phone +33 (0)5 61 33 26 17, fax +33 (0)5 61 33 25 60
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20160426/a604edfa/attachment.html>


More information about the users mailing list