//
// 在参数文件中添加， WRITE_FORMATTED_ELF:TRUE
// Material Studio中CASTEP，计算完的ELF《电子局域化函数》数据转成xsf格式，可以通用程序打开
// Created by hgzhang on 2022/10/10.
//
#include <iostream>
#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <dirent.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
#define MAX_BUF_SIZE 255
void trim(char *strIn, char *strOut);
int GetEleNum(char *str);
typedef class CELL {
public:
float ra[3],rb[3],rc[3];//save for future; vector;
float a,b,c;//lattice a,b,c
char atoms[255][255];
float pos[255][3];
int num;
int ReadCell(char *cellfilename);
} myCELL;
int CELL::ReadCell( char *cellfilename )
{
FILE *fp=NULL;
float x,y,z,fa,fb,fc;//fractional
size_t len=0;
ssize_t read;
char lable1[]="%BLOCK LATTICE_CART";
char lable2[]="%BLOCK POSITIONS_FRAC";
fp=fopen(cellfilename,"r");
if (fp==NULL) return 0;
char *str=(char *)malloc(MAX_BUF_SIZE);
char *tstr=(char *)malloc(MAX_BUF_SIZE);
char *vstr=(char *)malloc(MAX_BUF_SIZE);
getline(&str,&len,fp);trim(str,tstr);
if (strncmp(tstr, lable1,(int)strlen(lable1))!=0)
{
printf("\nNot Find %s in this file\n",lable1);
fclose(fp);free(str);free(tstr);free(vstr);
return 0;
}
getline(&str ,&len,fp);
getline(&tstr,&len,fp);
getline(&vstr,&len,fp);
sscanf(str ,"%f %f %f",&ra[0],&ra[1],&ra[2]);//y,z for temp storage;
sscanf(tstr,"%f %f %f",&rb[0],&rb[1],&rb[2]);//y,z for temp storage;
sscanf(vstr,"%f %f %f",&rc[0],&rc[1],&rc[2]);//y,z for temp storage;
a=sqrt(ra[0]*ra[0]+ra[1]*ra[1]+ra[2]*ra[2]);
b=sqrt(rb[0]*rb[0]+rb[1]*rb[1]+rb[2]*rb[2]);
c=sqrt(rc[0]*rc[0]+rc[1]*rc[1]+rc[2]*rc[2]);
while((read=getline(&str,&len,fp))!=-1)
{
trim(str,tstr);//printf(tstr);
if (strncmp(tstr, lable2,(int)strlen(lable2))==0)
{
printf("\nFind %s\nStart to read POSITION atom x,y,z\n",lable2);
break;
}
}
int i=0;
while((read=getline(&str,&len,fp))!=-1)
{
if (str[0]=='%') break; //end of this atom position section
trim(str,tstr); //printf(tstr);
strcpy(atoms[i++],tstr);
}
num=i;
for(int k=0;k<i;k++)
{
sscanf(atoms[k]+2,"%f %f %f",&fa,&fb,&fc);//read fractional coordinates;
pos[k][0]=fa*ra[0]+fb*rb[0]+fc*rc[0];
pos[k][1]=fa*ra[1]+fb*rb[1]+fc*rc[1];
pos[k][2]=fa*ra[2]+fb*rb[2]+fc*rc[2];
//printf("%-2d %.7f %.7f %.7f\n",GetEleNum(atoms[k]),fa,fb,fc);
}
//fprintf(out,"BEGIN_BLOCK_DATAGRID_3D\nDensities\nBEGIN_DATAGRID_3D_Density\n");
free(str);free(tstr);free(vstr);
fclose(fp);
return 1;
}
typedef class ELF {
public:
int ma,mb,mc;
int *x,*y,*z;
float *F,*E;
ELF();
int ReadELFFile(char *filename);
~ELF();
} myELF;
ELF::ELF() { ma=mb=mc=0;
x=NULL;y=NULL;z=NULL;
F=NULL;E=NULL;
}
ELF::~ELF(){
if (x!=NULL) free(x);
if (y!=NULL) free(y);
if (z!=NULL) free(z);
if (F!=NULL) free(F);
if (E!=NULL) free(E);
}
//Read all information from *.elf_fmt
int ELF::ReadELFFile(char *filename){
char strmatrix[]="fine FFT grid along <a,b,c>";
int m,n;
size_t len=0;
ssize_t read;
FILE *fpELF=NULL;
fpELF=fopen(filename,"r");
if (fpELF==NULL) return 0;
char *str=(char *)malloc(MAX_BUF_SIZE);
char *tstr=(char *)malloc(MAX_BUF_SIZE);
//Read ma*mb*mc
while((read=getline(&str,&len,fpELF))!=-1)
{
trim(str,tstr);// printf(tstr);
m=(int)strlen(strmatrix);
n=(int)strlen(tstr);
if (strstr((tstr),strmatrix))
{
//printf(tstr);
sscanf(tstr,"%d %d %d",&ma,&mb,&mc);
ma++;mb++;mc++; //50 --->51 for completeness of unit cell
printf("%s\nTo complete data in unit cells\n Save the ELF data:<%d,%d,%d>",tstr,ma,mb,mc);
}
else if (strncmp(tstr, "END header:",11)==0)
{
printf("\nFind END header:\nStart to read <a,b,c> chi\n");
break;
}
}
//read the large matrix <a,b,c, chi>
int i=0; //i indicates how many data are read?
if (ma>0&&mb>0&&mc>0)
{
int j=0;
x=(int*)malloc(sizeof(int)*ma*mb*mc);
y=(int*)malloc(sizeof(int)*ma*mb*mc);
z=(int*)malloc(sizeof(int)*ma*mb*mc);
E=(float*)malloc(sizeof(float)*ma*mb*mc);
memset(E, 0, sizeof(float)*ma*mb*mc);
F=(float*)malloc(sizeof(float)*ma*mb*mc);
memset(F, 0, sizeof(float)*ma*mb*mc);
while((read=getline(&str,&len,fpELF))!=-1)
{
trim(str,tstr);//printf(tstr);
if (strlen(tstr)<=5) continue;
sscanf(tstr,"%d %d %d %f",&x[i],&y[i],&z[i],&E[i]);
//printf(tstr); printf("[%d,%d,%d]=%f\n",x[i],y[i],z[i],E[i]);
//complete the matrix by adding the data along the other three faces of an unit cell.
j=i;
if (x[j]==1)
{ //printf("[%d,%d,%d]=%f\n",x[i],y[i],z[i],E[i]);
i++; x[i]=ma;y[i]=y[j];z[i]=z[j];E[i]=E[j];
if (y[j]==1) {
i++; x[i]=ma;y[i]=mb;z[i]=z[j];E[i]=E[j];
i++; x[i]=1;y[i]=mb;z[i]=z[j];E[i]=E[j];
}
if (z[j]==1) {
i++; x[i]=ma;y[i]=y[j];z[i]=mc;E[i]=E[j];
i++; x[i]=1;y[i]=y[j];z[i]=mc;E[i]=E[j];
}
}
else if (y[j]==1)
{ //printf("[%d,%d,%d]=%f\n",x[i],y[i],z[i],E[i]);
i++; x[i]=x[j];y[i]=mb;z[i]=z[j];E[i]=E[j];//printf("[%d,%d,%d]=%f\n",x[i],y[i],z[i],E[i]);
if (z[j]==1) {
i++; x[i]=x[j];y[i]=mb;z[i]=mc;E[i]=E[j];
i++; x[i]=x[j];y[i]=1;z[i]=mc;E[i]=E[j];
}
}
else if (z[j]==1 )
{ //printf("[%d,%d,%d]=%f\n",x[i],y[i],z[i],E[i]);
i++; x[i]=x[j];y[i]=y[j];z[i]=mc;E[i]=E[j];//printf("[%d,%d,%d]=%f\n",x[i],y[i],z[i],E[i]);
}
//not big deal if ignore the datapoints along three edge???
//printf("\n%s%d,%d,%d,%f",tstr,x[i],y[i],z[i],E[i]);
i++;
if (i>ma*mb*mc) break; //read too many data, check the file
}
printf("\nRead %d row ELF data",i);
printf("\n%d datapoints are needed",ma*mb*mc);
for(int k=0;k<i;k++){//reorganize the matrix following the xsf file format
F[(x[k]-1)+(y[k]-1)*ma+(z[k]-1)*ma*mb]=E[k];
//printf("[%d,%d,%d]=%f\n",x[k],y[k],z[k],F[k]);
}
}
fclose(fpELF);
return 1; //
}
int main(int argc, const char * argv[])
{
char ELFfilename[255];//"/Users/zhg/TEMP/convertELF/c2e/c2e/Al2O3.pot_fmt";
char CELLfilename[255];//"/Users/zhg/TEMP/convertELF/c2e/c2e/Al2O3.cell";
char XSFfilename[255];//"/Users/zhg/TEMP/convertELF/c2e/c2e/Al2O3_elf.xsf";
char ShareFullPathFile[255];
myCELL cl;
myELF elf;
printf("\n*********************************************************************");
printf("\nWhen use MS CASTEP, be sure to add \n WRITE_FORMATTED_ELF:TRUE");
printf("\ninto the SEEDNAME.param file\n\n");
printf("\nIf *.elf_fmt is not found!!\nCheck \"retain server file\" before submit task \nAfter done, download them from the server by yourself");
printf("\n*********************************************************************\n\n\n");
if (argc==2) {
char curpath[255]; int len=0;
getcwd(curpath, 255);
len=strlen(curpath);
curpath[len]='/';curpath[len+1]='\0';
printf("Current working directory: %s\n", curpath);
strcpy(ShareFullPathFile,curpath);
strcat(ShareFullPathFile,argv[1]);
strcpy(ELFfilename,ShareFullPathFile);strcpy(CELLfilename,ShareFullPathFile);strcpy(XSFfilename,ShareFullPathFile);
strcat(ELFfilename,".elf_fmt");strcat(CELLfilename,".cell");strcat(XSFfilename,"_ELF.xsf");
printf("\nConvert %s to\n %s \n using %s\n",ELFfilename,XSFfilename,CELLfilename);
if (cl.ReadCell(CELLfilename)==0) {printf("\nRead %s error!",CELLfilename);return 0;}
if (elf.ReadELFFile(ELFfilename)==0) {printf("\nRead %s error!",ELFfilename);return 0;}
}
else {
printf("\nPlease provide SEEDNAME");
printf("\nFor example: c2e Al2O3 ");
printf("\nwhich needs two file: SEEDNAME.elf_fmt and SEEDNAME.cell");;
printf("\nc2e iwll convert Al2O3.elf_fmt to Al2O3_elf.xsf");
printf("\nYou can use VESTA or Xcrysden to open SEEDNAME.xsf\n\n");
return 0;
}
FILE *fwELF=NULL;
fwELF=fopen(XSFfilename,"w");
if (fwELF==NULL) {printf("\n Can not open %s",XSFfilename);return 0;}
printf("\n Successfully open %s\n",XSFfilename);
//GetCell(fwELF, "/Users/zhg/TEMP/ELF/ELF/TiO2_rutile.cell",abc);
fprintf(fwELF,"#\n# Electron localization function\n#\nCRYSTAL\nPRIMVEC\n");
//printf("#\n# Electron localization function\n#\nCRYSTAL\nPRIMVEC\n");
fprintf(fwELF,"%15.7f%15.7f%15.7f\n",cl.ra[0],cl.ra[1],cl.ra[2]);
fprintf(fwELF,"%15.7f%15.7f%15.7f\n",cl.rb[0],cl.rb[1],cl.rb[2]);
fprintf(fwELF,"%15.7f%15.7f%15.7f\n",cl.rc[0],cl.rc[1],cl.rc[2]);
fprintf(fwELF,"CONVVEC\n");
fprintf(fwELF,"%15.7f%15.7f%15.7f\n",cl.ra[0],cl.ra[1],cl.ra[2]);
fprintf(fwELF,"%15.7f%15.7f%15.7f\n",cl.rb[0],cl.rb[1],cl.rb[2]);
fprintf(fwELF,"%15.7f%15.7f%15.7f\n",cl.rc[0],cl.rc[1],cl.rc[2]);
//fprintf(fwELF," %.7f %.7f %.7f\n",cl.a,0.0,0.0);
//fprintf(fwELF," %.7f %.7f %.7f\n",0.0,cl.b,0.0);
//fprintf(fwELF," %.7f %.7f %.7f\n",0.0,0.0,cl.c);
fprintf(fwELF,"PRIMCOORD\n");
//abc[0]=a;abc[1]=b;abc[2]=c;
fprintf(fwELF,"%d 1\n",cl.num);
for(int k=0;k<cl.num;k++)
fprintf(fwELF,"%-3d%15.7f%15.7f%15.7f\n",GetEleNum(cl.atoms[k]),cl.pos[k][0],cl.pos[k][1],cl.pos[k][2]);
fprintf(fwELF,"BEGIN_BLOCK_DATAGRID_3D\nELF\nBEGIN_DATAGRID_3D_Density\n");
fprintf(fwELF,"%d %d %d\n0.0 0.0 0.0\n",elf.ma,elf.mb,elf.mc);
fprintf(fwELF,"%15.7f%15.7f%15.7f\n",cl.ra[0],cl.ra[1],cl.ra[2]);
fprintf(fwELF,"%15.7f%15.7f%15.7f\n",cl.rb[0],cl.rb[1],cl.rb[2]);
fprintf(fwELF,"%15.7f%15.7f%15.7f\n",cl.rc[0],cl.rc[1],cl.rc[2]);
float max=-1E99,min=1E99;
for(int k=0;k<elf.ma*elf.mb*elf.mc;k++) {
fprintf(fwELF,"%15.7f\n",elf.F[k]);
if (elf.F[k]>max) max=elf.F[k];
if (elf.F[k]<min) min=elf.F[k];
}
printf("\nWrite %d datapoints,[min,max]=[%.7f,%.7f]\n",elf.ma*elf.mb*elf.mc,min,max);
fprintf(fwELF,"END_DATAGRID_3D\nEND_BLOCK_DATAGRID_3D");
fclose(fwELF);
printf("\n Successfully write %s\n",XSFfilename);
}// //end of main function;
//all function list
void trim(char *strIn, char *strOut)
{
int i, j ;
i = 0;
j = (int)strlen(strIn) - 1;
while(strIn[i] == ' ') ++i;
while(strIn[j] == ' ') --j;
strncpy(strOut, strIn + i , j - i + 1);
strOut[j - i + 1] = '\0';
}
int GetEleNum(char *str)
{
if (str[0]=='H' &&str[1]==' ') return 1;
if (str[0]=='L' &&str[1]=='i') return 3;
if (str[0]=='B' &&str[1]==' ') return 5;
if (str[0]=='C' &&str[1]==' ') return 6;
if (str[0]=='N' &&str[1]==' ') return 7;
if (str[0]=='O' &&str[1]==' ') return 8;
if (str[0]=='F' &&str[1]==' ') return 9;
if (str[0]=='N' &&str[1]=='a') return 11;
if (str[0]=='M' &&str[1]=='g') return 12;
if (str[0]=='A' &&str[1]=='l') return 13;
if (str[0]=='S' &&str[1]=='i') return 14;
if (str[0]=='P' &&str[1]==' ') return 15;
if (str[0]=='S' &&str[1]==' ') return 16;
if (str[0]=='C' &&str[1]=='l') return 17;
if (str[0]=='K' &&str[1]==' ') return 19;
if (str[0]=='C' &&str[1]=='a') return 20;
if (str[0]=='T' &&str[1]=='i') return 22;
if (str[0]=='V' &&str[1]==' ') return 23;
if (str[0]=='C' &&str[1]=='r') return 24;
if (str[0]=='M' &&str[1]=='n') return 25;
if (str[0]=='F' &&str[1]=='e') return 26;
if (str[0]=='C' &&str[1]=='o') return 27;
if (str[0]=='N' &&str[1]=='i') return 28;
if (str[0]=='C' &&str[1]=='u') return 29;
if (str[0]=='Z' &&str[1]=='n') return 30;
if (str[0]=='G' &&str[1]=='a') return 31;
if (str[0]=='G' &&str[1]=='e') return 32;
if (str[0]=='A' &&str[1]=='s') return 33;
if (str[0]=='S' &&str[1]=='e') return 34;
if (str[0]=='B' &&str[1]=='r') return 35;
if (str[0]=='R' &&str[1]=='b') return 37;
if (str[0]=='S' &&str[1]=='r') return 38;
if (str[0]=='Y' &&str[1]==' ') return 39;
if (str[0]=='Z' &&str[1]=='r') return 40;
if (str[0]=='N' &&str[1]=='b') return 41;
if (str[0]=='M' &&str[1]=='o') return 42;
if (str[0]=='R' &&str[1]=='u') return 44;
if (str[0]=='R' &&str[1]=='h') return 45;
if (str[0]=='P' &&str[1]=='d') return 46;
if (str[0]=='A' &&str[1]=='g') return 47;
if (str[0]=='C' &&str[1]=='d') return 48;
if (str[0]=='I' &&str[1]=='n') return 49;
if (str[0]=='S' &&str[1]=='n') return 50;
if (str[0]=='S' &&str[1]=='b') return 51;
if (str[0]=='T' &&str[1]=='e') return 52;
if (str[0]=='I' &&str[1]==' ') return 53;
if (str[0]=='C' &&str[1]=='s') return 55;
if (str[0]=='B' &&str[1]=='a') return 56;
if (str[0]=='C' &&str[1]=='e') return 58;
if (str[0]=='L' &&str[1]=='a') return 57;
if (str[0]=='H' &&str[1]=='f') return 72;
if (str[0]=='T' &&str[1]=='a') return 73;
if (str[0]=='W' &&str[1]==' ') return 74;
if (str[0]=='R' &&str[1]=='e') return 75;
if (str[0]=='O' &&str[1]=='s') return 76;
if (str[0]=='I' &&str[1]=='r') return 77;
if (str[0]=='P' &&str[1]=='t') return 78;
if (str[0]=='A' &&str[1]=='u') return 79;
if (str[0]=='H' &&str[1]=='h') return 80;
if (str[0]=='T' &&str[1]=='l') return 81;
if (str[0]=='P' &&str[1]=='b') return 82;
if (str[0]=='B' &&str[1]=='i') return 83;
if (str[0]=='P' &&str[1]=='o') return 84;
if (str[0]=='A' &&str[1]=='t') return 85;
return 0;
}