Relaxation run and resp fn runs: different force results?
Posted: Mon Jul 16, 2012 4:04 pm
Hi,
I have read this forum extensively, but am still striking out on some of my phonon calculations.
I have a structure which relaxes to a metal. My relaxation conditions are toldfe 1x10-12 and tolmxf 1x10-6.
I use occopt 7 and tsmear of 0.005.
I then attempt a resp. fn run using tolvrs 1.0-22. I copy/paste the acell and xred directly from my relaxation run script -- even using the "more" decimal places found in the *.log file as opposed to the *.out file. I use the typical 3 dataset script I have seen here in this forum and in the tutorials: first, do the tolvrs part, second do the ddk part, and third do the phonon part. [Admitedly there should be no need to do the second run as the structure is a metal, but it should not affect the phonon run adversely.] I have put the script I use at the end of this post. I use enough bands, as confirmed by the occ fillings I see in the ouput file (at least one is still 0 for each kpt. My ecut is locked in at 25 as I use OPIUM PSPS "optimized" to a certain kinetic energy.
What I see from the output from the first part of the first part of the resp fn run are forces that exceed that forces from my initial relaxation run by 1 or 2 orders of magnitude.
I do do one thing different from my relaxation run to my resp fn run: for the relax. run I use kpt grids that are dense, like 16x16x16. But for the resp fn, I lighten this up to 8x8x8. [sometimes, I lighten it up for speed, but other times, especially for hexagonal structures, I need to change the grid or the program basically does not find an answer to the ddk part (regardless if the structure is a metal or an insulator). For example, for a hexagonal type system, 6x6x6 works but a 16x16x16 (or 8x8x8) seems to just freeze. I beleive I use the correct shiftk: .5 .5 .5 for all 90/90/90 structures and 0 0 .5 for all 90/90/120 structs.]
Could this grid change be responsible for the problem? If so, what is the best action to be taken? Re-relax all structures with the lighter kpt grid and hope to achieve the 1x10-6 tolmxf with it?
Thanks,
JB
SCRIPT:
ndtset 3
#self-consistent run (do a relaxation prior to this to get the lattice params and xred you need)
iscf1 3
kptopt1 1
tolvrs1 10d-18
#Set 2 : Calculate the ddk s - needed for piezoelectric tensor and Born effective charges in dataset 3wf
getwfk2 -1
iscf2 -3
kptopt2 2
nqpt2 1
qpt2 0 0 0
rfelfd2 2
rfdir2 1 1 1
tolwfr2 1.0d-22
#Set 3 : response-function calculations for all needed perturbations
getddk3 -1
getwfk3 -2
iscf3 3
kptopt3 2
nqpt3 1
qpt3 0 0 0
rfphon3 1
rfatpol3 1 10
rfstrs3 3
rfdir3 1 1 1
tolvrs3 1.0d-10
#common input data
occopt 7
tsmear 0.005
angdeg 90 90 120
#Definition of the atom types and atoms
ntypat 3
znucl 82 22 16
natom 10
typat 2*1 2*2 6*3
nband 56
#nbdbuf 4
acell 1.21306650660996E+01 1.21306650660996E+01 1.05994861868081E+01
xred
3.33333333333333E-01 6.66666666666666E-01 7.50000000000000E-01
6.66666666666666E-01 3.33333333333333E-01 2.50000000000000E-01
0.00000000000000E+00 0.00000000000000E+00 2.50772128175421E-37
0.00000000000000E+00 6.12236641053275E-41 5.00000000000000E-01
3.48911517188473E-01 1.74455758594236E-01 7.50000000000000E-01
1.74455758594236E-01 3.48911517188473E-01 2.50000000000000E-01
8.25544241405764E-01 1.74455758594236E-01 7.50000000000000E-01
6.51088482811527E-01 8.25544241405764E-01 2.50000000000000E-01
8.25544241405764E-01 6.51088482811527E-01 7.50000000000000E-01
1.74455758594236E-01 8.25544241405764E-01 2.50000000000000E-01
diemac 12.0d0
nstep 120
ecut 25.0
ecutsm 0.5
#Exchange-correlation functional
ixc 2
#kptopt 1
ngkpt 6 6 6
nshiftk 1
shiftk 0 0 0.5
I have read this forum extensively, but am still striking out on some of my phonon calculations.
I have a structure which relaxes to a metal. My relaxation conditions are toldfe 1x10-12 and tolmxf 1x10-6.
I use occopt 7 and tsmear of 0.005.
I then attempt a resp. fn run using tolvrs 1.0-22. I copy/paste the acell and xred directly from my relaxation run script -- even using the "more" decimal places found in the *.log file as opposed to the *.out file. I use the typical 3 dataset script I have seen here in this forum and in the tutorials: first, do the tolvrs part, second do the ddk part, and third do the phonon part. [Admitedly there should be no need to do the second run as the structure is a metal, but it should not affect the phonon run adversely.] I have put the script I use at the end of this post. I use enough bands, as confirmed by the occ fillings I see in the ouput file (at least one is still 0 for each kpt. My ecut is locked in at 25 as I use OPIUM PSPS "optimized" to a certain kinetic energy.
What I see from the output from the first part of the first part of the resp fn run are forces that exceed that forces from my initial relaxation run by 1 or 2 orders of magnitude.
I do do one thing different from my relaxation run to my resp fn run: for the relax. run I use kpt grids that are dense, like 16x16x16. But for the resp fn, I lighten this up to 8x8x8. [sometimes, I lighten it up for speed, but other times, especially for hexagonal structures, I need to change the grid or the program basically does not find an answer to the ddk part (regardless if the structure is a metal or an insulator). For example, for a hexagonal type system, 6x6x6 works but a 16x16x16 (or 8x8x8) seems to just freeze. I beleive I use the correct shiftk: .5 .5 .5 for all 90/90/90 structures and 0 0 .5 for all 90/90/120 structs.]
Could this grid change be responsible for the problem? If so, what is the best action to be taken? Re-relax all structures with the lighter kpt grid and hope to achieve the 1x10-6 tolmxf with it?
Thanks,
JB
SCRIPT:
ndtset 3
#self-consistent run (do a relaxation prior to this to get the lattice params and xred you need)
iscf1 3
kptopt1 1
tolvrs1 10d-18
#Set 2 : Calculate the ddk s - needed for piezoelectric tensor and Born effective charges in dataset 3wf
getwfk2 -1
iscf2 -3
kptopt2 2
nqpt2 1
qpt2 0 0 0
rfelfd2 2
rfdir2 1 1 1
tolwfr2 1.0d-22
#Set 3 : response-function calculations for all needed perturbations
getddk3 -1
getwfk3 -2
iscf3 3
kptopt3 2
nqpt3 1
qpt3 0 0 0
rfphon3 1
rfatpol3 1 10
rfstrs3 3
rfdir3 1 1 1
tolvrs3 1.0d-10
#common input data
occopt 7
tsmear 0.005
angdeg 90 90 120
#Definition of the atom types and atoms
ntypat 3
znucl 82 22 16
natom 10
typat 2*1 2*2 6*3
nband 56
#nbdbuf 4
acell 1.21306650660996E+01 1.21306650660996E+01 1.05994861868081E+01
xred
3.33333333333333E-01 6.66666666666666E-01 7.50000000000000E-01
6.66666666666666E-01 3.33333333333333E-01 2.50000000000000E-01
0.00000000000000E+00 0.00000000000000E+00 2.50772128175421E-37
0.00000000000000E+00 6.12236641053275E-41 5.00000000000000E-01
3.48911517188473E-01 1.74455758594236E-01 7.50000000000000E-01
1.74455758594236E-01 3.48911517188473E-01 2.50000000000000E-01
8.25544241405764E-01 1.74455758594236E-01 7.50000000000000E-01
6.51088482811527E-01 8.25544241405764E-01 2.50000000000000E-01
8.25544241405764E-01 6.51088482811527E-01 7.50000000000000E-01
1.74455758594236E-01 8.25544241405764E-01 2.50000000000000E-01
diemac 12.0d0
nstep 120
ecut 25.0
ecutsm 0.5
#Exchange-correlation functional
ixc 2
#kptopt 1
ngkpt 6 6 6
nshiftk 1
shiftk 0 0 0.5