[BUG FIXED] Structural Optimization of FCC Ag using FHI PP

JEJohns · Post by **JEJohns** » Sun May 02, 2010 5:48 pm

I spent yesterday afternoon trying to setup to do a set of calculations of a 100 super cell of Ag(100). Before doing them, I wanted to test the psuedopotential I was using (the FHI GGA using PBE) so I set up a set of convergence studies for the silver FCC Structure. I have two questions:

1) While using the optcell 1 argument, I had to specify a very large dilatmx (~>1.3) because the routine kept trying to expand the unit cell way out on the first few broyden steps. Is this normal? I know from the tutorial this doesn't happen with LDA & Aluminum.

2) After running convergences in Ecut (8Ha - 28 Ha), the size of the kpoint grid (6x6x6 - 18x18x18), tsmear (0.2, 0.1, 0.05, 0.02, 0.01, 0.005), the unit cell converged to a lattice constant of 7.95 Bohr, and the actual value is well established at 7.729 Bohr. Is a 2.8% expansion typical of convergence for the unit cell of a noble metal? It seems a lot worse than the tutorial which got the convergence of aluminium to 0.2% with much softer parameters.

mverstra · Post by **mverstra** » Tue May 04, 2010 1:09 am

JEJohns wrote:1) While using the optcell 1 argument, I had to specify a very large dilatmx (~>1.3) because the routine kept trying to expand the unit cell way out on the first few broyden steps. Is this normal? I know from the tutorial this doesn't happen with LDA & Aluminum.

dilatmx doesn't keep it from expanding, it only tolerates much larger variations in cells. Your pseudopotential is probably bad if it is expanding the lattice way beyond the experimental value. What does it relax to? What happens with optcell 2?

2) After running convergences in Ecut (8Ha - 28 Ha), the size of the kpoint grid (6x6x6 - 18x18x18), tsmear (0.2, 0.1, 0.05, 0.02, 0.01, 0.005), the unit cell converged to a lattice constant of 7.95 Bohr, and the actual value is well established at 7.729 Bohr. Is a 2.8% expansion typical of convergence for the unit cell of a noble metal? It seems a lot worse than the tutorial which got the convergence of aluminium to 0.2% with much softer parameters.

d-metals will always be harder potentials. The precision is pretty bad indeed, but sometimes tolerated. Again, probably your potential is to blame. Where did you get it (which table?). I have used ftp://ftp.abinit.org/pub/abinitio/Psps/ ... Ag.GGA.fhi with decent results over the years (although not brilliant either, maybe 2% errors) - it requires a much higher ecut, though, as I recall.

Matthieu

JEJohns · Post by **JEJohns** » Tue May 04, 2010 6:47 pm

Thanks for the Response Matthieu

mverstra wrote:dilatmx doesn't keep it from expanding, it only tolerates much larger variations in cells. Your pseudopotential is probably bad if it is expanding the lattice way beyond the experimental value. What does it relax to? What happens with optcell 2?

It converges to 7.954 Bohr

d-metals will always be harder potentials. The precision is pretty bad indeed, but sometimes tolerated. Again, probably your potential is to blame. Where did you get it (which table?). I have used ftp://ftp.abinit.org/pub/abinitio/Psps/ ... Ag.GGA.fhi with decent results over the years (although not brilliant either, maybe 2% errors) - it requires a much higher ecut, though, as I recall.
Matthieu

I came across this article http://iopscience.iop.org/1367-2630/10/ ... 063020.pdf comparing the lattice parameters of d metals & some work functions using different GGA psuedopotentials within the PAW method. So it seems like 2-3% is a little high, but not much worse.

Also I was trying to calculate the work function for my slab, by taking the difference between the Fermi energy and the Potential at some point in the vacuum, so I put "prtpot 1" in the input file. It is now outputting a potential file for each broyden step and dataset (outo_geom_DS#_TIM#_POT), but when I open them with cut3d, the potential it reads out is zero everywhere. Is it a problem with my understanding of how prtpot should work or something else?

#Making the intial slab of 100 Ag using PBE IXC
#And Converged FCC Values from the FCC Structure
#By 7 slabs of vacuum, efermi had only coverged to the 20meV mark.
#but calc was taking much longer. Using the 7 vacuum layers & adding atoms
#looks like I used 6.5 slabs of vacuum. oops
kptopt 1
#Define a Unit Cell
#Definition of the unit cell
acell 7.954 7.954 7.954 #This is a 2.8% increase in unit cell size from
#experimental
rprim: 0.5 -0.5 0.0 # Unit Cell with 2 atoms in FCC geometry
0.5 0.5 0.0 #With 6.5 unit of vacuum
0.0 0.0 8.0
rprim+ 0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.5

ionmov 3
ecutsm 0.5
ndtset 6
iprcel 45
ntime 75
nstep 75
#Define the atoms involved
ntypat 1
znucl 47

#Define and position the atoms
natom: 4 natom+ 1 #Starting with 4 atoms for this Ag(100) slab
typat 1 1 1 1 1 1 1 1 1 1
nband: 30 nband+ 7

xcart 0.0 0.0 0.0
0.0 3.977 3.977
0.0 0.0 7.954
0.0 3.977 11.931
0.0 0.0 15.908
0.0 3.977 19.885
0.0 0.0 23.862
0.0 3.977 27.839
0.0 0.0 31.816
#Define the basis set

ecut 26

ngkpt 12 12 12 # This is a 2x2x2 FCC grid, based on the primitive vectors

#Exchange-correlation functional
ixc 11 #11=> GGA, Perdew-Burke-Ernzerhof GGA functional
tsmear 0.01
occopt 4
prtpot 1
toldff 1.0d-6 #difference in forces for SCF Convergence
tolmxf 1.0d-5 #Total difference in gradiants of forces for broyden
#convergence
#toldfe 1.0d-7
timopt 2

mverstra · Post by **mverstra** » Thu May 06, 2010 10:16 am

That is strange - try with a single dataset. Your trick of filling typat and the rest for the maximal number of atoms is clever but may be confusing abinit. If you examine the file produced by cut3d is it full of 0.0000D00 or some very small values /= 0?

Matthieu

JEJohns · Post by **JEJohns** » Thu May 06, 2010 9:16 pm

mverstra wrote:That is strange - try with a single dataset. Your trick of filling typat and the rest for the maximal number of atoms is clever but may be confusing abinit. If you examine the file produced by cut3d is it full of 0.0000D00 or some very small values /= 0?

Matthieu

I'm doing that calculation right now, but when I tell cut3d to print out the formatted 3d grid of potentials, I get numerical values. For example, in one run, my ngfft was 27 27 320, and the output of cut3d for the indexed 3D data (menu item #6), i get a block of numbers that are 233280rows by 4 columns, where the 4th column has reasonable values for the potential (27x27x320 is 233280). I assume that this corresponds to a fourier transform of the potential, so i don't know why it's not giving me the real space equivalent. Is there a non trivial way to back that information out from the data I have?

mverstra · Post by **mverstra** » Fri May 07, 2010 9:42 am

JEJohns wrote:
I'm doing that calculation right now, but when I tell cut3d to print out the formatted 3d grid of potentials, I get numerical values. For example, in one run, my ngfft was 27 27 320, and the output of cut3d for the indexed 3D data (menu item #6), i get a block of numbers that are 233280rows by 4 columns, where the 4th column has reasonable values for the potential (27x27x320 is 233280). I assume that this corresponds to a fourier transform of the potential, so i don't know why it's not giving me the real space equivalent. Is there a non trivial way to back that information out from the data I have?

Ok - this means your pot file is probably fine - the 3d indexed data is in real space I believe. So the problem is in the output method you chose. Perhaps you compiled with too much optimization. Give us your cut3d input values (which output file format were you using?), and output (unless it is too long). I mean the screen output of cut3d, not the file contents, which you said were 0.

Matthieu

JEJohns · Post by **JEJohns** » Sat May 08, 2010 9:10 pm

mverstra wrote:
JEJohns wrote:
I'm doing that calculation right now, but when I tell cut3d to print out the formatted 3d grid of potentials, I get numerical values. For example, in one run, my ngfft was 27 27 320, and the output of cut3d for the indexed 3D data (menu item #6), i get a block of numbers that are 233280rows by 4 columns, where the 4th column has reasonable values for the potential (27x27x320 is 233280). I assume that this corresponds to a fourier transform of the potential, so i don't know why it's not giving me the real space equivalent. Is there a non trivial way to back that information out from the data I have?

Ok - this means your pot file is probably fine - the 3d indexed data is in real space I believe. So the problem is in the output method you chose. Perhaps you compiled with too much optimization. Give us your cut3d input values (which output file format were you using?), and output (unless it is too long). I mean the screen output of cut3d, not the file contents, which you said were 0.

Matthieu

Ok, so here's the output from the single dataset calculation

Cut3d']
.Version 6.0.3 of CUT3D
.(MPI version, prepared for a x86_64_linux_gnu4.4 computer)

.Copyright (C) 1998-2010 ABINIT group .
CUT3D comes with ABSOLUTELY NO WARRANTY.
It is free software, and you are welcome to redistribute it
under certain conditions (GNU General Public License,
see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt).

ABINIT is a project of the Universite Catholique de Louvain,
Corning Inc. and other collaborators, see ~abinit/doc/developers/contributors.txt .
Please read ~abinit/doc/users/acknowledgments.html for suggested
acknowledgments of the ABINIT effort.
For more information, see http://www.abinit.org .

.Starting date : Thu 6 May 2010.

What is the name of the 3D function (density, potential or wavef) file ?
[/quote]
Ag100_7atoms12121o_geom_TIM12_POT
[quote="Cut3d']
=> Your 3D function file is : Ag100_7atoms12121o_geom_TIM12_POT
Does this file contain formatted 3D ASCII data (=0)
or unformatted binary header + 3D data (=1) ?
or ETSF binary (=2) ?
[/quote]
1
[quote="Cut3d wrote:===============================================================================
ECHO of the ABINIT file header

First record :
.codvsn,headform,fform = 6.0.3 57 102

Second record :
bantot,intxc,ixc,natom = 1218 0 11 7
ngfft(1:3),nkpt = 27 27 375 21
nspden,nspinor = 1 1
nsppol,nsym,npsp,ntypat = 1 16 1 1
occopt,pertcase,usepaw = 4 0 0
ecut,ecutdg,ecutsm = 2.6000000000E+01 2.6000000000E+01 5.0000000000E-01
ecut_eff = 2.6000000000E+01
qptn(1:3) = 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00
rprimd(1:3,1) = 3.9770000000E+00 -3.9770000000E+00 0.0000000000E+00
rprimd(1:3,2) = 3.9770000000E+00 3.9770000000E+00 0.0000000000E+00
rprimd(1:3,3) = 0.0000000000E+00 0.0000000000E+00 7.9540000000E+01
stmbias,tphysel,tsmear = 0.0000000000E+00 0.0000000000E+00 1.0000000000E-02

Third record :
istwfk= 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
nband = 58 58 58 58 58 58 58 58 58 58 58 58
58 58 58 58 58 58 58 58 58
npwarr=16048159841597615960159521593415992159541597015954
15946159641596615934159241592615900158381584015814
15718
so_psp= 1
symafm=
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
symrel=
1 0 0 0 1 0 0 0 1 -1 0 0 0 -1 0 0 0 -1
-1 0 0 0 1 0 0 0 -1 1 0 0 0 -1 0 0 0 1
-1 0 0 0 -1 0 0 0 1 1 0 0 0 1 0 0 0 -1
1 0 0 0 -1 0 0 0 -1 -1 0 0 0 1 0 0 0 1
0 1 0 1 0 0 0 0 1 0 -1 0 -1 0 0 0 0 -1
0 -1 0 1 0 0 0 0 -1 0 1 0 -1 0 0 0 0 1
0 -1 0 -1 0 0 0 0 1 0 1 0 1 0 0 0 0 -1
0 1 0 -1 0 0 0 0 -1 0 -1 0 1 0 0 0 0 1
type = 1 1 1 1 1 1 1
kptns = (max 50 k-points will be written)
4.166667E-02 4.166667E-02 5.000000E-01
1.250000E-01 4.166667E-02 5.000000E-01
2.083333E-01 4.166667E-02 5.000000E-01
2.916667E-01 4.166667E-02 5.000000E-01
3.750000E-01 4.166667E-02 5.000000E-01
4.583333E-01 4.166667E-02 5.000000E-01
1.250000E-01 1.250000E-01 5.000000E-01
2.083333E-01 1.250000E-01 5.000000E-01
2.916667E-01 1.250000E-01 5.000000E-01
3.750000E-01 1.250000E-01 5.000000E-01
4.583333E-01 1.250000E-01 5.000000E-01
2.083333E-01 2.083333E-01 5.000000E-01
2.916667E-01 2.083333E-01 5.000000E-01
3.750000E-01 2.083333E-01 5.000000E-01
4.583333E-01 2.083333E-01 5.000000E-01
2.916667E-01 2.916667E-01 5.000000E-01
3.750000E-01 2.916667E-01 5.000000E-01
4.583333E-01 2.916667E-01 5.000000E-01
3.750000E-01 3.750000E-01 5.000000E-01
4.583333E-01 3.750000E-01 5.000000E-01
4.583333E-01 4.583333E-01 5.000000E-01
wtk =
0.03 0.06 0.06 0.06 0.06 0.06 0.03 0.06 0.06 0.06
0.06 0.03 0.06 0.06 0.06 0.03 0.06 0.06 0.03 0.06
0.03
Then the occ = which shows ~41 occupied bands

tnons =
0.000000 0.000000 0.000000 0.000000 0.000000 0.300000
0.000000 0.000000 0.300000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 0.300000
0.000000 0.000000 0.300000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 0.300000
0.000000 0.000000 0.300000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 0.300000
0.000000 0.000000 0.300000 0.000000 0.000000 0.000000
znucl= 47.00

Pseudopotential info :
title=silver, fhi98PP : Trouiller-Martins-type, GGA Perdew/Burke/Ernzerhof (1996), l= 1 local
znuclpsp= 47.00, zionpsp= 11.00, pspso= 0, pspdat= 11001, pspcod= 6, pspxc= 11
lmnmax = 4

Last record :
residm,etot,fermie= 1.618471E-06 -2.532602731854E+02 -1.700730E-01
xred =
0.000000E+00 0.000000E+00 4.478977E-04
-5.000000E-01 5.000000E-01 4.953898E-02
7.006492E-46 3.503246E-46 9.974912E-02
-5.000000E-01 5.000000E-01 1.500000E-01
0.000000E+00 0.000000E+00 2.002509E-01
-5.000000E-01 5.000000E-01 2.504610E-01
0.000000E+00 0.000000E+00 2.995521E-01
End the ECHO of the ABINIT file header
===============================================================================

===========================================================

ECHO important input variables ...

Dimensional primitive vectors (ABINIT equivalent : rprimd):
3.977000E+00 -3.977000E+00 0.000000E+00
3.977000E+00 3.977000E+00 0.000000E+00
0.000000E+00 0.000000E+00 7.954000E+01
Grid density (ABINIT equivalent : ngfft): 27 27 375
Number of atoms : 7
Number of atomic types: 1

# Atomic positions (cartesian coordinates - Bohr)
1 0.000000E+00 0.000000E+00 3.562578E-02
2 0.000000E+00 3.977000E+00 3.940330E+00
3 4.179723E-45 -1.393241E-45 7.934045E+00
4 0.000000E+00 3.977000E+00 1.193100E+01
5 0.000000E+00 0.000000E+00 1.592796E+01
6 0.000000E+00 3.977000E+00 1.992167E+01
7 0.000000E+00 0.000000E+00 2.382637E+01

This file is a Density or Potential file

3D function was read. Ready for further treatment.
===========================================================

What is your choice ? Type:
0 => exit
1 => point (interpolation of data for a single point)
2 => line (interpolation of data along a line)
3 => plane (interpolation of data in a plane)
4 => volume (interpolation of data in a volume)
5 => 3D formatted data (output the bare 3D data - one column)
6 => 3D indexed data (bare 3D data, preceeded by 3D index)
7 => 3D Molekel formatted data
8 => 3D data with coordinates (tecplot ASCII format)
9 => output .xsf file for XCrysDen
10 => output .dx file for OpenDx
11 => compute atomic charge using the Hirshfeld method
12 => NetCDF file
14 => Gaussian/cube wavefunction module

1

Cut3d wrote: Your choice is 1

Select the coordinate system:
Type 1) for cartesian coordinates
or 2) for crystallographic coordinates

1

Cut3d wrote: Input point Cartesian Coord: X Y Z

0 0 0

Cut3d wrote:Crystallographic coordinates: 0.000000E+00 0.000000E+00 0.000000E+00
X coordinate, r1 is: 0.000000E+00
Y coordinate, r2 is: 0.000000E+00
Z coordinate, r3 is: 0.000000E+00

---------------------------------------------
Non-spin-polarized value= 0.000000E+00
---------------------------------------------
Task 1 has been done !

More analysis of the 3D file ? (1=default=yes,0=no)
[/quote

Similarly, if I use any Cartesian point, or crystallographic point, i still get a 0.0000000E+00 answer

[qoute="Cut3d']
Input point Cartesian Coord: X Y Z
1.2 1.3 14
Crystallographic coordinates: -1.257229E-02 3.143073E-01 1.760121E-01
X coordinate, r1 is: 9.874277E-01
Y coordinate, r2 is: 3.143073E-01
Z coordinate, r3 is: 1.760121E-01

---------------------------------------------
Non-spin-polarized value= 0.000000E+00
---------------------------------------------
Task 1 has been done !

[/quote]

But If I do this with the Density file from the same run, I get very sensible answers, for example, the Density at 0,0,0 is 4.3685E-04. Thanks again for all the help

mverstra · Post by **mverstra** » Sun May 09, 2010 12:58 pm

Ok, this is much more specific. Maybe the "point" interpolation routine has been broken, and it was not noticed. This is not something we do very routinely, as you usually want the whole cell potential on a grid for visualization or something.

Could you try the other options: line, volume interpolation, and perhaps the xsf/dx ouput files? If all of these are non 0.00 then we know it's the point interpolation and I can have a look at it.

Matthieu

JEJohns · Post by **JEJohns** » Sun May 09, 2010 5:27 pm

The point, Line, Plane, and Volume interpolations all give the potential as being pure ZERO. Outputting the potential to an opendx file, or NetCDF file, or indexed 3D or Raw 3D give reasonable potentials, once you visualize them properly. The *.cube file for gaussian loads the structure of the atoms correctly, but it is missing some header or some other information, and guassview can't load any surfaces, and I don't have a working install of XCrysDens to check that option.
--James

mverstra · Post by **mverstra** » Sun May 16, 2010 12:40 am

You are entirely correct! This bug only shows up for the point line plane interpolations, and for potentials only, densities are fine. I will add this patch to the 6.0.5 version, but in the mean time you can add the following line inside src/98_main/cut3d.F90, at line 515 (obviously without the +).

Great that this has come out!
Cheers,

Matthieu

Code: Select all

   else if(fform0==101 .or. fform0==102)then    ! Potential case
     if(ispden==0)then
       grid(:,:,:)=grid_full(:,:,:,1)
     else if(ispden==1 .or. ispden==2)then
       grid(:,:,:)=grid_full(:,:,:,ispden)
     else
       write(*,*)' Error: bad ispden value'
       stop
     end if
+     gridtt = grid
   end if

ABINIT Discussion Forums

[BUG FIXED] Structural Optimization of FCC Ag using FHI PP

[BUG FIXED] Structural Optimization of FCC Ag using FHI PP

Re: Structural Optimization of FCC Ag using FHI PP

Re: Structural Optimization of FCC Ag using FHI PP

Re: Structural Optimization of FCC Ag using FHI PP

Re: Structural Optimization of FCC Ag using FHI PP

Re: Structural Optimization of FCC Ag using FHI PP

Re: Structural Optimization of FCC Ag using FHI PP

Re: Structural Optimization of FCC Ag using FHI PP

Re: Structural Optimization of FCC Ag using FHI PP

Re: Structural Optimization of FCC Ag using FHI PP