Fermi surface

Queries about input and output files, running specific calculations, etc.


Moderators: Global Moderator, Moderator

Post Reply
Message
Author
User avatar
Kitaura
Newbie
Newbie
Posts: 10
Joined: Mon Oct 08, 2007 11:07 pm
License Nr.: 756

Fermi surface

#1 Post by Kitaura » Tue Sep 23, 2008 7:48 am

Dear VASP users,

I would like to compute and draw fermi surface of various materials using VASP code. I think information needed to draw fermi surface is in the wavecar file. Dose anyone tell me how to compute fermi surface from wavecar file?
Any comments and advice will be greatly appriciated.

Thanks in advance.
Last edited by Kitaura on Tue Sep 23, 2008 7:48 am, edited 1 time in total.
Top
admin
Administrator
Administrator
Posts: 2921
Joined: Tue Aug 03, 2004 8:18 am
License Nr.: 458

Fermi surface

#2 Post by admin » Tue Sep 23, 2008 9:54 am

you need a dense k-grid to interpolate the Fermi surface from the EIGENVAL file. WAVECAR is not required for that putpose.
Last edited by admin on Tue Sep 23, 2008 9:54 am, edited 1 time in total.
Top
User avatar
Kitaura
Newbie
Newbie
Posts: 10
Joined: Mon Oct 08, 2007 11:07 pm
License Nr.: 756

Fermi surface

#3 Post by Kitaura » Tue Sep 23, 2008 2:01 pm

Dear admin,

Thank you very much for your quick reply.

I usually use EIGENVAL file to draw band structure, and I didn't know how to visualize fermi surface from EIGENVAL file. So, please tell me how to visualize fermi surface from the file. (how to visualize a sphere in the case of alkali metal, for example.)

Thanks in advance.
Last edited by Kitaura on Tue Sep 23, 2008 2:01 pm, edited 1 time in total.
Top
lcyin
Newbie
Newbie
Posts: 31
Joined: Fri Mar 18, 2005 7:07 am
License Nr.: 375

Fermi surface

#4 Post by lcyin » Fri Sep 18, 2009 12:44 pm

Hi Kitaura,
I think what you want to do is to plot the isosurface of the electron charge density, is that right?
If so, you can check the Band decompose charge section of the VASP Manual.
Last edited by lcyin on Fri Sep 18, 2009 12:44 pm, edited 1 time in total.
Top
weberju

Fermi surface

#5 Post by weberju » Thu Jun 06, 2013 10:32 pm

Here is a short C++ code that will convert a VASP OUTCAR file into a *.bxsf Fermi surface file that can be used by xcrysden (http://www.xcrysden.org/doc/fermi.html). Keep in mind that you need to use really high k-point densities to get a smooth Fermi surface.

C++ code:
------------

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <string>
#include <strings.h>
#include <cstring>
#include <cstdlib>
using namespace std;

#define NKPTS_MAX 100000
#define NBANDS_MAX 2000

int main ( int argc, char *argv[])
{
if( argc < 4)
{
printf("INPUT FORMAT: \n");
printf("./v2xsf.exe nkx nky nkz \n");
exit(0);
}
int nkx, nky,nkz;
nkx=atoi(argv[1]);
nky=atoi(argv[2]);
nkz=atoi(argv[3]);

fstream fid;
fid.open("OUTCAR",ios::in);
if ( !fid )
{
printf("cant open OUTCAR\n");
exit(0);
}

char str[255], stmp[255];

//keyword1[]="NKPTS";
int NKPTS,NBANDS;
do {
fid.getline(str,255);
if(strstr(str,"NKPTS") != NULL )
{
printf("str=%s\n", str);
sscanf(str,"%s%s%s%d",stmp,stmp,stmp,&NKPTS);
sscanf(str,"%s%s%s%s%s%s%s%s%s%s%s%s%s%s%d",stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,&NBANDS);
printf("NKPTS: %d, NBANDS: %d\n",NKPTS,NBANDS);
break;
}
}while( !fid.eof()) ;



float A[3][3], B[3][3];
do {
fid.getline(str,255);
if(strstr(str,"reciprocal lattice vectors") != NULL )
{
for(int i=0;i<3;i++)
{
fid.getline(str,255);
sscanf(str,"%f%f%f%f%f%f",&A[0],&A[1],&A[2],&B[0],&B[1],&B[2]);
}
printf("A1: %f %f %f \n",A[0][0],A[0][1],A[0][2]);
printf("A2: %f %f %f \n",A[1][0],A[1][1],A[1][2]);
printf("A3: %f %f %f \n",A[2][0],A[2][1],A[2][2]);
printf("B1: %f %f %f \n",B[0][0],B[0][1],B[0][2]);
printf("B2: %f %f %f \n",B[1][0],B[1][1],B[1][2]);
printf("B3: %f %f %f \n",B[2][0],B[2][1],B[2][2]);

break;
}
}while( !fid.eof());



float Ef;
do {
fid.getline(str,255);
if(strstr(str,"E-fermi") != NULL )
{
printf("str=%s\n", str);
sscanf(str,"%s%s%f",stmp,stmp,&Ef);
printf("Ef: %f\n",Ef);
break;
}
}while( !fid.eof()) ;



fid.getline(str,255);
fid.getline(str,255);

float kpt[NKPTS_MAX][3];
float Egv[NKPTS_MAX][NBANDS_MAX];
for(int ik =0; ik < NKPTS; ik++)
{
// k-point 1: kx ky kz
fid.getline(str,255);
sscanf(str,"%s%s%s%f%f%f",stmp,stmp,stmp,&kpt[ik][0],&kpt[ik][1],&kpt[ik][2]);
fid.getline(str,255);

for(int ib=0;ib<NBANDS;ib++)
{
int bandno,occ;
fid.getline(str,255);
sscanf(str,"%d%f%f",&bandno,&Egv[ik][ib],&occ);
// printf("%d %d : %f\n", ik,ib,Egv[ik][ib]);
}

fid.getline(str,255);
}

fid.close();


// find those cross Ef

int Nf=0;
static int nstore[NBANDS_MAX];
for( int ib =0 ; ib < NBANDS; ib ++)
{
float vmin=10000, vmax=-10000;
for( int ik =0 ; ik < NKPTS ; ik++ )
{
if ( Egv[ik][ib] < vmin)
{ vmin = Egv[ik][ib];}
if( Egv[ik][ib] > vmax )
{ vmax= Egv[ik][ib];}
}
if( vmax > Ef && vmin < Ef)
{ nstore[Nf++]=ib; }

}
int ib;
float Egvf[NKPTS][Nf];
for(int i=0; i<Nf; i++)
{
ib=nstore;
for( int ik=0 ; ik < NKPTS ; ik ++ )
{
Egvf[ik]=Egv[ik][ib];
}
}

//save to bxsf

FILE *fid2;
fid2=fopen("FS.bxsf","w");
if( fid2 == NULL)
{
printf( "can't create fout: FS.bxsf\n");
exit(0);
}

fprintf(fid2,"BEGIN_INFO\n");
fprintf(fid2," Fermi Energy: %7.5f\n",Ef);
fprintf(fid2,"END_INFO\n");


fprintf(fid2,"\n");

fprintf(fid2,"BEGIN_BLOCK_BANDGRID_3D\n");
fprintf(fid2," one_word_comment \n");
fprintf(fid2," BEGIN_BANDGRID_3D\n");
fprintf(fid2," %d \n",Nf);
fprintf(fid2," %d %d %d \n",nkx,nky,nkz);
fprintf(fid2," 0 0 0\n"); // origin Gamma
fprintf(fid2," %f %f %f \n", B[0][0],B[0][1],B[0][2]);//Vec(1,:));
fprintf(fid2," %f %f %f \n", B[1][0],B[1][1],B[1][2]);
fprintf(fid2," %f %f %f \n", B[2][0],B[2][1],B[2][2]);

int ik = 0;
for(int i=0; i < Nf; i++)
{
fprintf(fid2, " BAND: %d \n",nstore+1);
for(int ix=0 ; ix < nkx; ix ++)
{
for(int iy=0 ; iy < nky; iy ++)
{
for(int iz=0 ; iz <nkz; iz++)
{
fprintf(fid2," %f ",Egvf[ik++]);
}
fprintf(fid2,"\n");
}
fprintf(fid2,"\n");
}

}

fprintf(fid2," END_BANDGRID_3D\n");
fprintf(fid2,"END_BLOCK_BANDGRID_3D\n");


fclose(fid2);

}

<span class='smallblacktext'>[ Edited Thu Jun 06 2013, 10:33PM ]</span>
Last edited by weberju on Thu Jun 06, 2013 10:32 pm, edited 1 time in total.
Top
beck
Newbie
Newbie
Posts: 19
Joined: Tue Jun 14, 2005 3:22 pm

Fermi surface

#6 Post by beck » Wed Jun 12, 2013 6:31 am

Dear weberju, all

your nice code is great, but was a little buggy and could be improved a bit, so here it is:

improved C++ code from weberju:
(sorry for the "?" - I dont know where they come from! Just replace all of them with spaces)
------------------------

Code: Select all

#include?<stdio.h>
#include?<stdlib.h>
#include?<math.h>
#include?<iostream>
#include?<fstream>
#include?<string>
#include?<strings.h>
#include?<cstring>
#include?<cstdlib>
using?namespace?std;

int?main?(?int?argc,?char?**argv)
{
?int?nkx,?nky,nkz;

??char?str[255],?stmp[255];

??fstream?fid2a;
??fid2a.open("KPOINTS",ios::in);
??if?(?!fid2a?)
????{
??????printf("cant?open?KPOINTS\n");
??????exit(0);
????}
????printf("open?KPOINTS?successfully!\n");
????fid2a.getline(str,255);
????fid2a.getline(str,255);
????fid2a.getline(str,255);
????fid2a.getline(str,255);
????sscanf(str,"%d%d%d",&nkx,&nky,&nkz);
????printf("KPOINTS?found:?%d?%d?%d\n",nkx,nky,nkz);

??fstream?fid;
??fid.open("OUTCAR",ios::in);
??if?(?!fid?)
????{
??????printf("cant?open?OUTCAR\n");
??????exit(0);
????}
????printf("open?OUTCAR?successfully!\n");
??
??unsigned?long?NKPTS,NBANDS;
??do?{
????fid.getline(str,255);
????if(strstr(str,"NKPTS")?!=?NULL?)
??????{
	sscanf(str,"%s%s%s%ld",stmp,stmp,stmp,&NKPTS);
	sscanf(str,"%s%s%s%s%s%s%s%s%s%s%s%s%s%s%ld",stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,&NBANDS);
	printf("K-points?found:?%ld\n",NKPTS);
	if((nkx*nky*nkz)?==?NKPTS)
	??{
	????printf("nkx/y/z?and?total?number?of?K-points?consistent!?GOOD!\n");
	??}else{
	??printf("nkx/y/z?and?total?number?of?K-points?inconsistent!\nStopping?here!\n");
	??exit(-1);
	}
	printf("Total?number?of?BANDS?found:?%ld\n",NBANDS);
	break;
??????}
??}while(?!fid.eof())?;
??
??float?A[3][3],?B[3][3];
??do?{
????fid.getline(str,255);
????if(strstr(str,"reciprocal?lattice?vectors")?!=?NULL?)
??????{
	for(int?i=0;i<3;i++)
	??{
	????fid.getline(str,255);
	????sscanf(str,"%f%f%f%f%f%f",&A[i][0],&A[i][1],&A[i][2],&B[i][0],&B[i][1],&B[i][2]);
	??}
	printf("\nRealspace?Cell:\n%f?%f?%f?\n",A[0][0],A[0][1],A[0][2]);
	printf("%f?%f?%f?\n",A[1][0],A[1][1],A[1][2]);
	printf("%f?%f?%f?\n",A[2][0],A[2][1],A[2][2]);
	printf("Brillouinzone:\n%f?%f?%f?\n",B[0][0],B[0][1],B[0][2]);
	printf("%f?%f?%f?\n",B[1][0],B[1][1],B[1][2]);
	printf("%f?%f?%f?\n\n",B[2][0],B[2][1],B[2][2]);
	
	break;
??????}
??}while(?!fid.eof());
????
??float?Ef;
??do?{
????fid.getline(str,255);
????if(strstr(str,"E-fermi")?!=?NULL?)
??????{
	//	printf("str=%s\n",?str);
	sscanf(str,"%s%s%f",stmp,stmp,&Ef);
	printf("Fermienergy:?%f\n",Ef);
	break;
??????}
??}while(?!fid.eof())?;

??printf("Space?Requirements:?%ld?%ld?->?total?(Byte):?%ld\n", NKPTS*3*sizeof(double), NKPTS*NBANDS*sizeof(double),?NKPTS*3*sizeof(double)+ NKPTS*NBANDS*sizeof(double));
??
??float?**kpt;
??float?**Egv;
??float?**Egvf;
??kpt?=?new?float*[NKPTS];
??Egv?=?new?float*[NKPTS];
??for(unsigned?int?i1=0;i1<NKPTS;i1++)
????{
??????kpt[i1]?=?new?float[3];
??????Egv[i1]?=?new?float[NBANDS];
????}
??for(int?ik?=0;?ik?<?NKPTS;?ik++)
????{
??????fid.getline(str,255);
??????sscanf(str,"%s%s%s%f%f%f",stmp,stmp,stmp,&kpt[ik][0],&kpt[ik][1],&kpt[ik][2]);
??????fid.getline(str,255);
??????
??????for(int?ib=0;ib<NBANDS;ib++)
	{
	??int?bandno,occ;
	??fid.getline(str,255);
	??sscanf(str,"%d%f%d",&bandno,&Egv[ik][ib],&occ);
	}
??????
??????fid.getline(str,255);
????}
??
??fid.close();

??//?find?those?cross?Ef
??
??int?Nf=0;
??int?*nstore;
??nstore?=?new?int[NBANDS];
??for(?int?ib?=0?;?ib?<?NBANDS;?ib?++)
????{
??????float?vmin=10000,?vmax=-10000;
??????for(?int?ik?=0?;?ik?<?NKPTS?;?ik++?)
	{
	??if?(?Egv[ik][ib]?<?vmin)
	????{?vmin?=?Egv[ik][ib];}
	??if(?Egv[ik][ib]?>?vmax?)
	????{?vmax=?Egv[ik][ib];}
	}
??????if(?vmax?>?Ef?&&?vmin?<?Ef)
	{?nstore[Nf++]=ib;?}
??????
????}
??printf("Number?of?Bands?crossing?E-Fermi:?%d\n",Nf);
??Egvf?=?new?float*[NKPTS];
??for(unsigned?int?i1=0;i1<NKPTS;i1++)
????{
??????Egvf[i1]?=?new?float[Nf];
????}
??int?ib;
??for(int?i=0;?i<Nf;?i++)
????{
??????ib=nstore[i];
??????for(?int?ik=0?;?ik?<?NKPTS?;?ik?++?)
	{
	??Egvf[ik][i]=Egv[ik][ib];
	}
????}
??
??//save?to?bxsf
??
??FILE?*fid2;
??fid2=fopen("FermiSurface.bxsf","w");
??if(?fid2?==?NULL)
????{
??????printf(?"can't?create?fout:?FS.bxsf\n");
??????exit(0);
????}
??
??fprintf(fid2,"BEGIN_INFO\n");
??fprintf(fid2,"?Fermi?Energy:?%7.5f\n",Ef);
??fprintf(fid2,"END_INFO\n");


??fprintf(fid2,"\n");

??fprintf(fid2,"BEGIN_BLOCK_BANDGRID_3D\n");
??fprintf(fid2,"?Bugfixed_by_mbB_11.06.2013\n");
??fprintf(fid2,"?BEGIN_BANDGRID_3D\n");
??fprintf(fid2,"?%d?\n",Nf);
??fprintf(fid2,"?%d?%d?%d?\n",nkx,nky,nkz);
??fprintf(fid2,"?0?0?0\n");?//?origin?Gamma
??fprintf(fid2,"?%f?%f?%f?\n",?B[0][0],B[0][1],B[0][2]);
??fprintf(fid2,"?%f?%f?%f?\n",?B[1][0],B[1][1],B[1][2]);
??fprintf(fid2,"?%f?%f?%f?\n",?B[2][0],B[2][1],B[2][2]);

??unsigned?long?ik?=?0;
??for(int?i=0;?i?<?Nf;?i++)
????{
??????ik?=?0;
??????fprintf(fid2,?"?BAND:?%d?\n",nstore[i]+1);
??????for(int?ix=0?;?ix?<?nkx;?ix?++)
	{
	??for(int?iy=0?;?iy?<?nky;?iy?++)
	????{
	??????for(int?iz=0?;?iz?<nkz;?iz++)
		{
		??fprintf(fid2,"?%f?",Egvf[ik++][i]);
		}
	??????fprintf(fid2,"\n");
	????}
	??fprintf(fid2,"\n");
	}
??????
????}
??
??fprintf(fid2,"?END_BANDGRID_3D\n");
??fprintf(fid2,"END_BLOCK_BANDGRID_3D\n");
??
??printf("\nGenerated?File?FermiSurface.bxsf.?Now?start\n");
??printf("xcrysden?--bxsf?FermiSurface.bxsf\n");
??fclose(fid2);
}
----------------------------
end C++ code. Save it (eg. to fermisurface.cpp) and
compile it with the command:
g++ -O3 -lm -o fermisurface fermisurface.cpp

Now, you want to make a VASP calculation using the non-collinear version with ISYM=0 and LSORBIT=.TRUE. in the INCAR using - as already pointed out - more k-points than normal. After successful finishing it just run "fermisurface" in the folder of the calculation. Then you can visualize the fermi surface using xcrysden (the last line of output from the code repeats the correct command line for it).

Thanks again weberju! Finally somebody showed how easy the fermi surface can be generated!

Regards,
beck

Code: Select all

Code: Select all

Code: Select all







<span class='smallblacktext'>[ Edited Wed Jun 12 2013, 06:41AM ]</span>
Last edited by beck on Wed Jun 12, 2013 6:31 am, edited 1 time in total.
Top
hatdau
Newbie
Newbie
Posts: 39
Joined: Thu Jul 03, 2008 12:04 am
Location: Michigan State University
Contact:
Contact hatdau

Fermi surface

#7 Post by hatdau » Wed Jun 12, 2013 6:42 pm

That's great job you guys have done! However above code just works with one iteration calculation. You have to add one loop to take care of multiple iteration one.
Best,
Last edited by hatdau on Wed Jun 12, 2013 6:42 pm, edited 1 time in total.
Go green, Go white
Top
beck
Newbie
Newbie
Posts: 19
Joined: Tue Jun 14, 2005 3:22 pm

Fermi surface

#8 Post by beck » Fri Jun 14, 2013 1:00 pm

Dear hatdau,

> However above code just works with one iteration calculation.
> You have to add one loop
could have done by you also, but here we are. And since I am unable to get rid of the "?" whatever I try please contact me via email to get the code.

Regards.

<span class='smallblacktext'>[ Edited Fri Jun 14 2013, 01:07PM ]</span>
Last edited by beck on Fri Jun 14, 2013 1:00 pm, edited 1 time in total.
Top
hatdau
Newbie
Newbie
Posts: 39
Joined: Thu Jul 03, 2008 12:04 am
Location: Michigan State University
Contact:
Contact hatdau

Fermi surface

#9 Post by hatdau » Fri Jun 14, 2013 10:01 pm

Hi Beck,
I have been able to modify the code to take care of it as well.
Thank you so much.
Best,
Last edited by hatdau on Fri Jun 14, 2013 10:01 pm, edited 1 time in total.
Go green, Go white
Top
yhzhaosicnu

Fermi surface

#10 Post by yhzhaosicnu » Fri Aug 02, 2013 10:25 am

Dear weberju,

Thanks for giving us the code. There is still some questions which I need your help.

1. Should I use the KPOINTS in self-consistent calculation or band structure calculation?

2. This code give me the following bxsf file, which is not a complete one. Can you help me? Thanks a lot.

----------------------------------------------------------------------------------------------
BEGIN_INFO
Fermi Energy: 5.44160
END_INFO

BEGIN_BLOCK_BANDGRID_3D
one_word_comment
BEGIN_BANDGRID_3D
0
40 40 6
0 0 0
0.307542 0.177560 0.000000
0.000000 0.355119 0.000000
0.000000 0.000000 0.039179
END_BANDGRID_3D
END_BLOCK_BANDGRID_3D
-------------------------------------------------------------------------------------------------------------------

Best,

Yonghong

[quote="atoi(argv[1"]);
nky=atoi(argv[2]);
nkz=atoi(argv[3]);

fstream fid;
fid.open("OUTCAR",ios::in);
if ( !fid )
{
printf("cant open OUTCAR\n");
exit(0);
}

char str[255], stmp[255];

//keyword1[]="NKPTS";
int NKPTS,NBANDS;
do {
fid.getline(str,255);
if(strstr(str,"NKPTS") != NULL )
{
printf("str=%s\n", str);
sscanf(str,"%s%s%s%d",stmp,stmp,stmp,&NKPTS);
sscanf(str,"%s%s%s%s%s%s%s%s%s%s%s%s%s%s%d",stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,&NBANDS);
printf("NKPTS: %d, NBANDS: %d\n",NKPTS,NBANDS);
break;
}
}while( !fid.eof()) ;



float A[3][3], B[3][3];
do {
fid.getline(str,255);
if(strstr(str,"reciprocal lattice vectors") != NULL )
{
for(int i=0;i<3;i++)
{
fid.getline(str,255);
sscanf(str,"%f%f%f%f%f%f",&A[0],&A[1],&A[2],&B[0],&B[1],&B[2]);
}
printf("A1: %f %f %f \n",A[0][0],A[0][1],A[0][2]);
printf("A2: %f %f %f \n",A[1][0],A[1][1],A[1][2]);
printf("A3: %f %f %f \n",A[2][0],A[2][1],A[2][2]);
printf("B1: %f %f %f \n",B[0][0],B[0][1],B[0][2]);
printf("B2: %f %f %f \n",B[1][0],B[1][1],B[1][2]);
printf("B3: %f %f %f \n",B[2][0],B[2][1],B[2][2]);

break;
}
}while( !fid.eof());



float Ef;
do {
fid.getline(str,255);
if(strstr(str,"E-fermi") != NULL )
{
printf("str=%s\n", str);
sscanf(str,"%s%s%f",stmp,stmp,&Ef);
printf("Ef: %f\n",Ef);
break;
}
}while( !fid.eof()) ;



fid.getline(str,255);
fid.getline(str,255);

float kpt[NKPTS_MAX][3];
float Egv[NKPTS_MAX][NBANDS_MAX];
for(int ik =0; ik < NKPTS; ik++)
{
// k-point 1: kx ky kz
fid.getline(str,255);
sscanf(str,"%s%s%s%f%f%f",stmp,stmp,stmp,&kpt[ik][0],&kpt[ik][1],&kpt[ik][2]);
fid.getline(str,255);

for(int ib=0;ib<NBANDS;ib++)
{
int bandno,occ;
fid.getline(str,255);
sscanf(str,"%d%f%f",&bandno,&Egv[ik][ib],&occ);
// printf("%d %d : %f\n", ik,ib,Egv[ik][ib]);
}

fid.getline(str,255);
}

fid.close();


// find those cross Ef

int Nf=0;
static int nstore[NBANDS_MAX];
for( int ib =0 ; ib < NBANDS; ib ++)
{
float vmin=10000, vmax=-10000;
for( int ik =0 ; ik < NKPTS ; ik++ )
{
if ( Egv[ik][ib] < vmin)
{ vmin = Egv[ik][ib];}
if( Egv[ik][ib] > vmax )
{ vmax= Egv[ik][ib];}
}
if( vmax > Ef && vmin < Ef)
{ nstore[Nf++]=ib; }

}
int ib;
float Egvf[NKPTS][Nf];
for(int i=0; i<Nf; i++)
{
ib=nstore;
for( int ik=0 ; ik < NKPTS ; ik ++ )
{
Egvf[ik]=Egv[ik][ib];
}
}

//save to bxsf

FILE *fid2;
fid2=fopen("FS.bxsf","w");
if( fid2 == NULL)
{
printf( "can't create fout: FS.bxsf\n");
exit(0);
}

fprintf(fid2,"BEGIN_INFO\n");
fprintf(fid2," Fermi Energy: %7.5f\n",Ef);
fprintf(fid2,"END_INFO\n");


fprintf(fid2,"\n");

fprintf(fid2,"BEGIN_BLOCK_BANDGRID_3D\n");
fprintf(fid2," one_word_comment \n");
fprintf(fid2," BEGIN_BANDGRID_3D\n");
fprintf(fid2," %d \n",Nf);
fprintf(fid2," %d %d %d \n",nkx,nky,nkz);
fprintf(fid2," 0 0 0\n"); // origin Gamma
fprintf(fid2," %f %f %f \n", B[0][0],B[0][1],B[0][2]);//Vec(1,:));
fprintf(fid2," %f %f %f \n", B[1][0],B[1][1],B[1][2]);
fprintf(fid2," %f %f %f \n", B[2][0],B[2][1],B[2][2]);

int ik = 0;
for(int i=0; i < Nf; i++)
{
fprintf(fid2, " BAND: %d \n",nstore+1);
for(int ix=0 ; ix < nkx; ix ++)
{
for(int iy=0 ; iy < nky; iy ++)
{
for(int iz=0 ; iz <nkz; iz++)
{
fprintf(fid2," %f ",Egvf[ik++]);
}
fprintf(fid2,"\n");
}
fprintf(fid2,"\n");
}

}

fprintf(fid2," END_BANDGRID_3D\n");
fprintf(fid2,"END_BLOCK_BANDGRID_3D\n");


fclose(fid2);

}

[/quote]
Last edited by yhzhaosicnu on Fri Aug 02, 2013 10:25 am, edited 1 time in total.
Top
weberju

Fermi surface

#11 Post by weberju » Mon Sep 23, 2013 10:52 pm

yhzhaosicnu,

Beck sent me the updated code that he modified (Thanks a lot Beck!). Please try this version as it fixed a lot of the issues with the old version. As pointed out, make sure that the symmetry is turned off. It shouldn't matter if it is a self-consistent calculation or not. I hope this helps.


/////////////////////////////////////////////
//
// FermiSurface.cpp : Code to extract FermiSurface
// from a VASP OUTCAR (calculated with ISYM=0 and
// LSORBIT=.TRUE. using the noncollinear version of VASP
//
// Original from weberju, 11.06.2013
// Bugfixed and extended: beck, 14.06.2013
// Bugfixed and extended to 2 spins: beck, 28.06.2013
// added usage/--help option: beck, 06.09.2013
//
// LICENSE:
// GPLv3 License-http://www.gnu.org/licenses/gpl.html
//
/////////////////////////////////////////////
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <string>
#include <strings.h>
#include <cstring>
#include <cstdlib>
using namespace std;

int main ( int argc, char **argv)
{
int nkx, nky,nkz;
int run;
char str[255], stmp[255];
fstream fid2a;
char fname1[255],fname2[255];
bool run_kpoints;

sprintf(fname1,"OUTCAR");
sprintf(fname2,"KPOINTS");
run = 0;
run_kpoints = true;

if(argc == 2 && (strstr(argv[1],"--help") || strstr(argv[1],"-h")))
{
printf("Usage:\nFermiSurface (for all standard filenames)\n");
printf("FermiSurface <OUTCAR-NAME> (for all standard filenames expect the OUTCAR-file, which correct name is given)\n");
printf("FermiSurface <KPOINTS-NAME> <OUTCAR-NAME> (for the correct OUTCAR and KPOINTS filenames)\n");
printf("FermiSurface nkx nky nkz (for all standard filenames with the k-point mesh given here)\n");
printf("FermiSurface nkx nky nkz <OUTCAR-NAME> (k-point mesh and a different OUTCAR filename)\n");
exit(0);
}
switch(argc)
{
default:
break;
case 2: // just OUTCAR name is given
sprintf(fname1,"%s",argv[1]);
break;
case 3: // KPOINTS and OUTCAR name/path are given
sprintf(fname2,"%s",argv[1]);
sprintf(fname1,"%s",argv[2]);
break;
case 4: // only nkx nky nkz are given
nkx = atoi(argv[1]);
nky = atoi(argv[2]);
nkz = atoi(argv[3]);
run_kpoints = false;
break;
case 5: // nkx nky nkz and OUTCAR-name given
nkx = atoi(argv[1]);
nky = atoi(argv[2]);
nkz = atoi(argv[3]);
run_kpoints = false;
sprintf(fname1,"%s",argv[4]);
break;
}
if(run_kpoints)
{
fid2a.open(fname2,ios::in);
if ( !fid2a )
{
printf("cant open KPOINTS\n");
exit(0);
}
printf("open KPOINTS successfully!\n");
fid2a.getline(str,255);
fid2a.getline(str,255);
fid2a.getline(str,255);
fid2a.getline(str,255);
sscanf(str,"%d%d%d",&nkx,&nky,&nkz);
fid2a.close();
}
printf("KPOINTS found: %d %d %d\n",nkx,nky,nkz);

fstream fid;
fid.open(fname1,ios::in);
if ( !fid )
{
printf("cant open OUTCAR\n");
exit(0);
}
printf("open OUTCAR successfully!\n");

unsigned long NKPTS,NBANDS;
do {
fid.getline(str,255);
if(strstr(str,"NKPTS") != NULL )
{
sscanf(str,"%s%s%s%ld",stmp,stmp,stmp,&NKPTS);
sscanf(str,"%s%s%s%s%s%s%s%s%s%s%s%s%s%s%ld",stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,stmp,&NBANDS);
if(run == 0)
printf("K-points found: %ld\n",NKPTS);
if((nkx*nky*nkz) == NKPTS)
{
if(run == 0)
printf("nkx/y/z and total number of K-points consistent! GOOD!\n");
}else{
printf("nkx/y/z and total number of K-points inconsistent!\nStopping here!\n");
exit(-1);
}
if(run == 0)
printf("Total number of BANDS found: %ld\n",NBANDS);
break;
}
}while( !fid.eof()) ;

do{
float A[3][3], B[3][3];
do {
fid.getline(str,255);
if(strstr(str,"reciprocal lattice vectors") != NULL )
{
for(int i=0;i<3;i++)
{
fid.getline(str,255);
sscanf(str,"%f%f%f%f%f%f",&A[0],&A[1],&A[2],&B[0],&B[1],&B[2]);
}
if(run == 0)
{
printf("\nRealspace Cell:\n%f %f %f \n",A[0][0],A[0][1],A[0][2]);
printf("%f %f %f \n",A[1][0],A[1][1],A[1][2]);
printf("%f %f %f \n",A[2][0],A[2][1],A[2][2]);
printf("Brillouinzone:\n%f %f %f \n",B[0][0],B[0][1],B[0][2]);
printf("%f %f %f \n",B[1][0],B[1][1],B[1][2]);
printf("%f %f %f \n\n",B[2][0],B[2][1],B[2][2]);
}
break;
}
}while( !fid.eof());

float Ef;
do {
fid.getline(str,255);
if(strstr(str,"E-fermi") != NULL )
{
sscanf(str,"%s%s%f",stmp,stmp,&Ef);
if(run == 0)
printf("Fermienergy: %f\n",Ef);
break;
}
}while( !fid.eof()) ;

if(run == 0)
printf("Space Requirements: %ld %ld -> total (Byte): %ld\n",NKPTS*3*sizeof(double),NKPTS*NBANDS*sizeof(double),NKPTS*3*sizeof(double)+NKPTS*NBANDS*sizeof(double));

if(!fid.eof()){
int do_spin;
do_spin = 1;
fid.getline(str,255);
while(strlen(str) < 1 && !fid.eof())
fid.getline(str,255);

if(strstr(str,"spin") != NULL)
{
printf("2 Spins found!\n");
do_spin = 2;
fid.getline(str,255);
while(strlen(str) < 1 && !fid.eof())
fid.getline(str,255);
}
if(!fid.eof())
{
for(int is =0; is < do_spin;is++)
{
float **kpt;
float **Egv;
float **Egvf;
kpt = new float*[NKPTS];
Egv = new float*[NKPTS];
for(unsigned int i1=0;i1<NKPTS;i1++)
{
kpt[i1] = new float[3];
Egv[i1] = new float[NBANDS];
}
for(int ik =0; ik < NKPTS; ik++)
{
sscanf(str,"%s%s%s%f%f%f",stmp,stmp,stmp,&kpt[ik][0],&kpt[ik][1],&kpt[ik][2]);
fid.getline(str,255);

for(int ib=0;ib<NBANDS;ib++)
{
int bandno,occ;
fid.getline(str,255);
sscanf(str,"%d%f%d",&bandno,&Egv[ik][ib],&occ);
}

fid.getline(str,255);
fid.getline(str,255);
}

// find those cross Ef
int Nf=0;
int *nstore;
nstore = new int[NBANDS];
for( int ib =0 ; ib < NBANDS; ib ++)
{
float vmin=10000, vmax=-10000;
for( int ik =0 ; ik < NKPTS ; ik++ )
{
if ( Egv[ik][ib] < vmin)
{ vmin = Egv[ik][ib];}
if( Egv[ik][ib] > vmax )
{ vmax= Egv[ik][ib];}
}
if( vmax > Ef && vmin < Ef)
{
nstore[Nf++]=ib;
}
}
if(run == 0)
printf("Number of Bands crossing E-Fermi: %d\n",Nf);
Egvf = new float*[NKPTS];
for(unsigned int i1=0;i1<NKPTS;i1++)
{
Egvf[i1] = new float[Nf];
}
int ib;
for(int i=0; i<Nf; i++)
{
ib=nstore;
for( int ik=0 ; ik < NKPTS ; ik ++ )
{
Egvf[ik]=Egv[ik][ib];
}
}

//save to bxsf

if(Nf > 0)
{
if(is == 0)
run++;
FILE *fid2;
char fname[200];
sprintf(fname,"FermiSurface_%d_%d.bxsf",run,is+1);
fid2=fopen(fname,"w");
if( fid2 == NULL)
{
printf( "can't create fout: FS.bxsf\n");
exit(0);
}

fprintf(fid2,"BEGIN_INFO\n");
fprintf(fid2," Fermi Energy: %7.5f\n",Ef);
fprintf(fid2,"END_INFO\n\n");

fprintf(fid2,"BEGIN_BLOCK_BANDGRID_3D\n");
fprintf(fid2," Bugfixed_by_mbB_28.06.2013\n");
fprintf(fid2," BEGIN_BANDGRID_3D\n");
fprintf(fid2," %d \n",Nf);
fprintf(fid2," %d %d %d \n",nkx,nky,nkz);
fprintf(fid2," 0 0 0\n"); // origin Gamma
fprintf(fid2," %f %f %f \n", B[0][0],B[0][1],B[0][2]);//Vec(1,:));
fprintf(fid2," %f %f %f \n", B[1][0],B[1][1],B[1][2]);
fprintf(fid2," %f %f %f \n", B[2][0],B[2][1],B[2][2]);

unsigned long ik = 0;
for(int i=0; i < Nf; i++)
{
ik = 0;
fprintf(fid2, " BAND: %d \n",nstore+1);
for(int ix=0 ; ix < nkx; ix ++)
{
for(int iy=0 ; iy < nky; iy ++)
{
for(int iz=0 ; iz <nkz; iz++)
{
fprintf(fid2," %f ",Egvf[ik++]);
}
fprintf(fid2,"\n");
}
fprintf(fid2,"\n");
}
}

fprintf(fid2," END_BANDGRID_3D\n");
fprintf(fid2,"END_BLOCK_BANDGRID_3D\n");
fclose(fid2);
delete kpt;
delete Egv;
delete Egvf;
if(run == 1 && is == 0)
printf("\nGenerated File FermiSurface_StepX_SpinY.bxsf. Now start\n");
printf("xcrysden --bxsf %s\n",fname);
}

if(is == 0 && do_spin ==2)
{
fid.getline(str,255);
fid.getline(str,255);
}
}
}
}
}while( !fid.eof()) ;

fid.close();
}
Last edited by weberju on Mon Sep 23, 2013 10:52 pm, edited 1 time in total.
Top
urpach
Newbie
Newbie
Posts: 1
Joined: Tue Jan 14, 2014 11:11 am
License Nr.: 424

Fermi surface

#12 Post by urpach » Tue Jan 14, 2014 11:45 am

Dear weberju and beck

Thanks for the updated code! However I still have some problems. By running "fermisurface" after the VASP calculation (VASP 4.6, ISYM=0, LSORBIT=.TRUE.) I got this:


open KPOINTS successfully!
KPOINTS found: 14 14 8
open OUTCAR successfully!
K-points found: 1568
nkx/y/z and total number of K-points consistent! GOOD!
Total number of BANDS found: 140733193388032

Realspace Cell:
3.684753 0.000000 0.000000
0.000000 3.684753 0.000000
0.000000 0.000000 6.141896
Brillouinzone:
0.271389 0.000000 0.000000
0.000000 0.271389 0.000000
0.000000 0.000000 0.162816

Fermienergy: 3.270600
Space Requirements: 37632 1765357177859473408 -> total (Byte): 1765357177859511040
terminate called after throwing an instance of 'std::bad_alloc'
what(): std::bad_alloc
Aborted (core dumped)


Thus the number of BANDS and therefore also the Space Requirements are far too large (I set NBANDS = 40 in the INCAR file). Do yo have any idea what could be the problem? Thanks in advance!
Last edited by urpach on Tue Jan 14, 2014 11:45 am, edited 1 time in total.
Top
Post Reply

Return to “Using VASP”