I am trying to simulate Fe at high pressure. so I have to modify muffin-tin radius in Fe.in. I would like to know whether this is the only parameter that I need to modify or there are other related parameters. A related question is whether the muffin-tin radius should be kept constant for all pressures or changed according to cell size? For the moment, I have some difficulty for obtaining the correct lattice c/a ratio for hcp iron. the c/a versus energy curve is not smooth enough and c/a at the lowest energy does not agree well with experimental data. I have used rgkmax=8.0, gmaxvr=16.0, epspot=1.0e-6, and ngrid=8 8 6.
My third question concerns the local orbitals, does inclusion of the local orbital is invoked automatically in calculation?
Best regards!
zhangyg
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
(1) run the program "species" to change muffin-tin radius.
(2) increase cutoff "rgkmax" to make energy-volume curve smooth.
(3) local orbital invoked automatically.
with best wishes!
zhangyg
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
You should only need to change the muffin-tin radius in Fe.in without running species. The other parameters (determined by the species code) should be almost unchanged with respect to radius.
In principle, the default species file should give good behaviour. Were you using autormt=.true. when you obtained a jagged curve?
Best wishes
Kay.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Yes, when I used autormt=.true. the energy~c/a curve is not very smooth (I meaned the slopes at higher and lower c/a are very different, and that the c/a at the lowest energy did not correspond well to experiments). The problem can be remedied by using fixed muffin-tin radius. A better cutoff and k-points gives even better curve and a c/a in good agreement with experiments.
When I run "species" with Rmt=1.8, nothing changed except sprmax (from 42.4036 to 41.5249). the c/a obtained is about the same for Rmt=2.0 and 1.8. But the energy is different in these two cases.
Because of thess observations, I plan to use a fixed Rmt (at 1.8) for all simulations from low to high pressures. I will report to you the results once they come out.
Best regards!
zhangyg
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Here I want to tell you some of my results on hcp iron.
As a first step, I want to study Pressure-Volume relation of hcp iron at high pressures and T=0. To do that, first at each volume I made simulations at different lattice ratio c/a to find the c/a that corresponds to the lowest energy. Then the energies found and the corresponding volumes are fitted with an EOS to get bulk modulus and its pressure derivative, which permit me to draw a P-V curve. Unfortunately, the P-V curve show a large difference with experimental data (e.g., Mao et al, 1990, J. Geophys. Res.) and with previous theoretical studies using similar techniques (e.g.,Stixrude et al, 1994, PRB 50:6442~), especially when P is increased the difference becomes larger (pressure predicted by EXCITING is much too low). The c/a ratio is about right though, which is ~1.58.
I have used the same procedure for the code PWscf. The results I obtained are quite similar with previous experimental and theoretical studies that I mentioned above. So I really do not know the reason for the problem I encountered here and seek help from you.
Following is a typical "exciting.in" that I used:
tasks
0
!max length of G+k vectors
rgkmax
9.0
!max length of G for exapnding the interstitial density and potential
gmaxvr
18.0
I think you have made a mistake with your second atom: it should be at (0.5 0.5 0.5), to make the structure hcp.
If I use this (along with the default Fe.in but with rmt=1.9; a c/a ratio of 1.6; GGA; and a shifted 4x4x4 k-point set) then I get the following energy-volume set:
I think I am still confused here, because textbooks taught me that for a hcp primitive cell, the a-axis is at (1,0,0),the b-axis is at (-0.5,0.866,0), and the c-axis is (0,0,1); The a- and b-axis have 120 degree apart; the first atom is at (000), the second atom is at (1/3,2/3,1/2). You said that the second atom is at (1/2,1/2,1/2), I am still thinking why this is so.
Curiously, when I put (1/2,1/2,1/2) and (1/3,2/3,1/2) separatly into Exciting, and calculate the forces between atoms, they are always zero! The total energy is similar, with the (1/3,2/3,1/2) case slightly at lower energy. The number of crystal symmetry is 12 with no inversion for (1/3,2/3,1/2) and is 8 with inversion.
I have to think more to tell you.
With best regards!
zhangyg
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
I got correct results now. It was due to an error of mine in writing the EOS equation (3rd-order Birch-Murnaghan) for fitting V~E. The results are now in very good agreement with previous experimental and theoretical studies.
Thanks a lot for having the chance to discuss with you!
With best regards!
zhangyg
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Dear exciting users:
I am trying to simulate Fe at high pressure. so I have to modify muffin-tin radius in Fe.in. I would like to know whether this is the only parameter that I need to modify or there are other related parameters. A related question is whether the muffin-tin radius should be kept constant for all pressures or changed according to cell size? For the moment, I have some difficulty for obtaining the correct lattice c/a ratio for hcp iron. the c/a versus energy curve is not smooth enough and c/a at the lowest energy does not agree well with experimental data. I have used rgkmax=8.0, gmaxvr=16.0, epspot=1.0e-6, and ngrid=8 8 6.
My third question concerns the local orbitals, does inclusion of the local orbital is invoked automatically in calculation?
Best regards!
zhangyg
I think I solved my problem.
(1) run the program "species" to change muffin-tin radius.
(2) increase cutoff "rgkmax" to make energy-volume curve smooth.
(3) local orbital invoked automatically.
with best wishes!
zhangyg
Hi Zhangyg,
You should only need to change the muffin-tin radius in Fe.in without running species. The other parameters (determined by the species code) should be almost unchanged with respect to radius.
In principle, the default species file should give good behaviour. Were you using autormt=.true. when you obtained a jagged curve?
Best wishes
Kay.
Dear Kay:
Thanks for your kind response.
Yes, when I used autormt=.true. the energy~c/a curve is not very smooth (I meaned the slopes at higher and lower c/a are very different, and that the c/a at the lowest energy did not correspond well to experiments). The problem can be remedied by using fixed muffin-tin radius. A better cutoff and k-points gives even better curve and a c/a in good agreement with experiments.
When I run "species" with Rmt=1.8, nothing changed except sprmax (from 42.4036 to 41.5249). the c/a obtained is about the same for Rmt=2.0 and 1.8. But the energy is different in these two cases.
Because of thess observations, I plan to use a fixed Rmt (at 1.8) for all simulations from low to high pressures. I will report to you the results once they come out.
Best regards!
zhangyg
Dear Kay and EXCITING users:
Here I want to tell you some of my results on hcp iron.
As a first step, I want to study Pressure-Volume relation of hcp iron at high pressures and T=0. To do that, first at each volume I made simulations at different lattice ratio c/a to find the c/a that corresponds to the lowest energy. Then the energies found and the corresponding volumes are fitted with an EOS to get bulk modulus and its pressure derivative, which permit me to draw a P-V curve. Unfortunately, the P-V curve show a large difference with experimental data (e.g., Mao et al, 1990, J. Geophys. Res.) and with previous theoretical studies using similar techniques (e.g.,Stixrude et al, 1994, PRB 50:6442~), especially when P is increased the difference becomes larger (pressure predicted by EXCITING is much too low). The c/a ratio is about right though, which is ~1.58.
I have used the same procedure for the code PWscf. The results I obtained are quite similar with previous experimental and theoretical studies that I mentioned above. So I really do not know the reason for the problem I encountered here and seek help from you.
Following is a typical "exciting.in" that I used:
tasks
0
!max length of G+k vectors
rgkmax
9.0
!max length of G for exapnding the interstitial density and potential
gmaxvr
18.0
beta0
0.05
nempty
10
epspot
1.e-7
xctype
20
stype
1
swidth
0.01
!hcp lattice vectors
avec
1.0 0.0 0.0
-0.5 0.86602540378444 0.0
0.0 0.0 1.0
scale1
3.97070781412664025010
scale2
3.97070781412664025010
scale3
6.59137497145022281516
sppath
'/home/zhang/program/exciting/species/'
atoms
1 : nspecies
'Fe.in' : spfname
2 : natoms
0.000000000000 0.000000000000 0.0 0.0 0.0 0.0 : atposl, bfcmt
0.333333333333 0.666666666667 0.5 0.0 0.0 0.0 : atposl, bfcmt
ngridk
16 16 10
vkloff
0.0 0.0 0.5
I am looking forward to your kind help!
zhangyg
Dear Zhangyg,
I think you have made a mistake with your second atom: it should be at (0.5 0.5 0.5), to make the structure hcp.
If I use this (along with the default Fe.in but with rmt=1.9; a c/a ratio of 1.6; GGA; and a shifted 4x4x4 k-point set) then I get the following energy-volume set:
183.030484874810 -2545.52372768911
172.266838178647 -2545.53723003408
162.136179425880 -2545.54779202249
152.601283896320 -2545.55576959869
143.627115979077 -2545.55939598936
135.180700435540 -2545.55817631546
127.231001233119 -2545.55063111897
119.748807504522 -2545.53763120660
112.706626213535 -2545.51327710843
Using the Vinet equation of state, this gives:
V0=70.50462 a.u.^3
B0=297.8336 GPa
B0'=5.43962
This compares well with Gerd Steinle-Neumann et al. (Phys. Rev. B 60, 791 - 799 (1999)):
V0=69.0 a.u.^3
B0=292 GPa
B0'=4.4
although the B0' is somewhat different (experiment is 5.3). This may be due to my sparse k-point set or the different equation of state.
Let me know what you think.
Cheers,
Kay.
Dear Kay:
Thanks a lot for your calculations.
I think I am still confused here, because textbooks taught me that for a hcp primitive cell, the a-axis is at (1,0,0),the b-axis is at (-0.5,0.866,0), and the c-axis is (0,0,1); The a- and b-axis have 120 degree apart; the first atom is at (000), the second atom is at (1/3,2/3,1/2). You said that the second atom is at (1/2,1/2,1/2), I am still thinking why this is so.
Curiously, when I put (1/2,1/2,1/2) and (1/3,2/3,1/2) separatly into Exciting, and calculate the forces between atoms, they are always zero! The total energy is similar, with the (1/3,2/3,1/2) case slightly at lower energy. The number of crystal symmetry is 12 with no inversion for (1/3,2/3,1/2) and is 8 with inversion.
I have to think more to tell you.
With best regards!
zhangyg
You're absolutely right, and I'm absolutely wrong!
The coordinates are (1/3,2/3,1/2).
I'll let you know what the new results are.
Cheers,
Kay.
Dear Kay:
I got correct results now. It was due to an error of mine in writing the EOS equation (3rd-order Birch-Murnaghan) for fitting V~E. The results are now in very good agreement with previous experimental and theoretical studies.
Thanks a lot for having the chance to discuss with you!
With best regards!
zhangyg
Dear zhangyg,
Glad to hear it works OK!
Best wishes,
Kay.