[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