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();
}