Relaxation run and resp fn runs: different force results?

Phonons, DFPT, electron-phonon, electric-field response, mechanical response…

Moderators: mverstra, joaocarloscabreu

Locked
brehmj
Posts: 35
Joined: Thu Jan 20, 2011 3:18 pm

Relaxation run and resp fn runs: different force results?

Post by brehmj » 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

ilukacevic
Posts: 271
Joined: Sat Jan 16, 2010 12:05 pm
Location: Dept. of Physics, University J. J. Strossmayer, Osijek, Croatia
Contact:

Re: Relaxation run and resp fn runs: different force result

Post by ilukacevic » Tue Jul 17, 2012 9:13 am

Hi!

I think that it might not be just the difference in ngkpt, but also in the tolerance criterion. Tolvrs 10^-22 can give much more converged quantities than toldfe 10^-12. Because of that you get different forces between relax and rf. You should actually use the same input parameters (ecut, ngkpt, ecutsm, dilatmx, tolerances, grid shifts, scf options, etc.) for both types of calculation, unless you're absolutely sure that your quantities are completely converged.
Besides, I don't know if it makes any differences, but your diemac is too low for metals. The value of 12 is optimal for Si semiconductor. you should leave it to default.

Hope this helps.

Igor L.

User avatar
jzwanzig
Posts: 504
Joined: Mon Aug 17, 2009 9:25 am

Re: Relaxation run and resp fn runs: different force result

Post by jzwanzig » Tue Jul 17, 2012 6:26 pm

Igor is completely right--it is particularly important to keep ecutsm the same between the relaxation and the new runs with the relaxed structure. Toldfe is probably the weakest criterion, unless all the atoms are on special positions I use toldff for relaxtion and tolvrs for potentials. Also, concerning diemac, setting it wrong doesn't wreck the calculations, just makes the initial preconditioning step less efficient. For a metal leave it at its default, for a semiconductor like Si use 10 or so, for a wide gap insulator use 3.
Josef W. Zwanziger
Professor, Department of Chemistry
Canada Research Chair in NMR Studies of Materials
Dalhousie University
Halifax, NS B3H 4J3 Canada
jzwanzig@gmail.com

maxim
Posts: 78
Joined: Wed May 19, 2010 1:17 pm
Location: Institute of Silicate Chemistry of Russian Academy of Sciences, Saint-Petersburg, Russia

Re: Relaxation run and resp fn runs: different force result

Post by maxim » Mon Sep 10, 2012 11:05 am

Dear abinit users and developers,
I am trying to solve the same problem. I try to explore some of the properties of KNO3, for which I first decided to relax this structure.

brehmj wrote:actually use the same input parameters (ecut, ngkpt, ecutsm, dilatmx, tolerances, grid shifts, scf options, etc.) for both types of calculation

2jzwanzig: toldff, tolrff are not allowed for RF calculations. Maybe in order to keep the same tolerance input parameter for both types of calculation is it better to use tolvrs instead? And why no one uses tolwfr?

jzwanzig wrote:unless all the atoms are on special positions I use toldff for relaxtion

A newbie question: what is a special position - is it a position where ALL THREE coordinates unique (NOT 0, 1, 1/2, 1/3, 2/3, etc.)? And in my code these are atoms No7-15:

Code: Select all

xred 0.0000000000 0.0000000000 0.0000000000
0.6666666667 0.3333333333 0.3333333333
0.3333333333 0.6666666667 0.6666666667
0.0000000000 0.0000000000 0.4050000000
0.3333333333 0.6666666667 0.0716700000
0.6666666667 0.3333333333 0.7383300000
0.0713500000 0.5356400000 0.1006700000
0.4643600000 0.5356400000 0.1006700000
0.4643600000 0.9286500000 0.1006700000
0.1309800000 0.2620200000 0.4340000000
0.1309800000 0.8690200000 0.4340000000
0.7379800000 0.8690200000 0.4340000000
0.4046600000 0.2023300000 0.7673300000
0.7976700000 0.2023300000 0.7673300000
0.7976700000 0.5953400000 0.7673300000
?
Is it true that forces on all other atoms with general positions will disappear (when using toldff, tolrff)? Or not? Or these atoms with special positions will already break the balance of zero forces?

The last question about occupations: I have 47 bands only in my output file, and much more k-points, and then what does the phrase means:
brehmj wrote:occ fillings I see in the ouput file (at least one is still 0 for each kpt

if e.g. i have 2048 k-points and 47 bands only:

Code: Select all

              occ5     2.000000  2.000000  2.000000  2.000000  2.000000  2.000000
                       2.000000  2.000000  2.000000  2.000000  2.000000  2.000000
                       2.000000  2.000000  2.000000  2.000000  2.000000  2.000000
                       2.000000  2.000000  2.000000  2.000000  2.000000  2.000000
                       2.000000  2.000000  2.000000  2.000000  2.000000  2.000000
                       2.000000  2.000000  2.000000  2.000000  2.000000  2.000000
                       2.000000  2.000000  2.000000  2.000000  2.000000  2.000000
                       2.000000  2.000000  2.000000  0.000000  0.000000
?
Does that phrase have meaning only for another value of kptopt variable and does my values of occ ​​are normal?

Thank you for your time!
Attachments
t41.in
(2.09 KiB) Downloaded 342 times
M.Yu. Arsent'ev
Institute of Silicate Chemistry of RAS
tikhonov_p-a@mail.ru

maxim
Posts: 78
Joined: Wed May 19, 2010 1:17 pm
Location: Institute of Silicate Chemistry of Russian Academy of Sciences, Saint-Petersburg, Russia

Re: Relaxation run and resp fn runs: different force result

Post by maxim » Mon Jan 07, 2013 1:17 pm

jzwanzig wrote:Also, concerning diemac, setting it wrong doesn't wreck the calculations, just makes the initial preconditioning step less efficient. For a metal leave it at its default, for a semiconductor like Si use 10 or so, for a wide gap insulator use 3.

And what if we are dealing with ferroelectrics - their dielectric constant is much higher than 10.
Which value of diemac should I use: very high original, or value corresponding to their electronic conductivity (3-10)?
M.Yu. Arsent'ev
Institute of Silicate Chemistry of RAS
tikhonov_p-a@mail.ru

Locked