Skip to content

Instantly share code, notes, and snippets.

@Ncmexp2717
Created Oct 31, 2021
Embed
What would you like to do?
Modified for a cube file including spaces in comment lines
/**********************************************************************
gcube2oned.c
gcube2oned.c is a program which transforms the data of 3D cube data
to an 1D data along a chosen direction by integrating over the
remaining 2D.
Usage:
./gcube2oned input.cube
Log of gcube2oned.c:
05/Oct./2019 Released by T. Ozaki
01/Nov./2021 Modified by N. Yamaguchi
***********************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#define fp_bsize 1048576 /* buffer size for setvbuf */
#define Hartree2eV 27.2113845
#define BohrR 0.529177249 /* Angstrom */
int main(int argc, char *argv[])
{
int i,j,itmp,n1,n2,n3,po,axis;
int atomnum1,atomnum2;
int Ngrid1_1,Ngrid1_2,Ngrid1_3;
int Ngrid2_1,Ngrid2_2,Ngrid2_3;
double **Gxyz1,**Gxyz2;
double ***CubeData1,***CubeData2;
double Grid_Origin1[4];
double Grid_Origin2[4];
double gtv1[4][4],gtv2[4][4];
double dtmp,step;
char ctmp[100];
char buf[1000],buf2[1000],*c;
FILE *fp1,*fp2;
char fp_buf[fp_bsize]; /* setvbuf */
if (argc!=3){
printf("Usage:\n");
printf(" ./gcube2oned input.cube axis_number\n");
exit(0);
}
axis = atoi(argv[2]);
/*******************************************
read the first file
*******************************************/
if ((fp1 = fopen(argv[1],"r")) != NULL){
setvbuf(fp1,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
/* scanf cube tile */
/* Disabled by N. Yamaguchi ***
fscanf(fp1,"%s",ctmp);
fscanf(fp1,"%s",ctmp);
* ***/
/* Added by N. Yamaguchi ***/
fgets(ctmp,100,fp1);
fgets(ctmp,100,fp1);
/* ***/
fscanf(fp1,"%d",&atomnum1);
fscanf(fp1,"%lf %lf %lf",&Grid_Origin1[1],&Grid_Origin1[2],&Grid_Origin1[3]);
fscanf(fp1,"%d",&Ngrid1_1);
fscanf(fp1,"%lf %lf %lf",&gtv1[1][1],&gtv1[1][2],&gtv1[1][3]);
fscanf(fp1,"%d",&Ngrid1_2);
fscanf(fp1,"%lf %lf %lf",&gtv1[2][1],&gtv1[2][2],&gtv1[2][3]);
fscanf(fp1,"%d",&Ngrid1_3);
fscanf(fp1,"%lf %lf %lf",&gtv1[3][1],&gtv1[3][2],&gtv1[3][3]);
/* allocation of arrays */
Gxyz1 = (double**)malloc(sizeof(double*)*(atomnum1+1));
for (i=0; i<(atomnum1+1); i++){
Gxyz1[i] = (double*)malloc(sizeof(double)*4);
}
CubeData1 = (double***)malloc(sizeof(double**)*Ngrid1_1);
for (i=0; i<Ngrid1_1; i++){
CubeData1[i] = (double**)malloc(sizeof(double*)*Ngrid1_2);
for (j=0; j<Ngrid1_2; j++){
CubeData1[i][j] = (double*)malloc(sizeof(double)*Ngrid1_3);
}
}
/* scanf xyz coordinates */
for (i=1; i<=atomnum1; i++){
fscanf(fp1,"%lf %lf %lf %lf %lf",&Gxyz1[i][0],&dtmp,&Gxyz1[i][1],&Gxyz1[i][2],&Gxyz1[i][3]);
}
/* scanf cube data */
for (n1=0; n1<Ngrid1_1; n1++){
for (n2=0; n2<Ngrid1_2; n2++){
for (n3=0; n3<Ngrid1_3; n3++){
fscanf(fp1,"%lf",&CubeData1[n1][n2][n3]);
}
}
}
fclose(fp1);
}
else{
printf("error in scanfing %s\n",argv[1]);
}
/*******************************************
make a file
*******************************************/
step = BohrR*sqrt(gtv1[axis][1]*gtv1[axis][1]+gtv1[axis][2]*gtv1[axis][2]+gtv1[axis][3]*gtv1[axis][3]);
if (axis==1){
for (n1=0; n1<Ngrid1_1; n1++){
double sum = 0.0;
for (n2=0; n2<Ngrid1_2; n2++){
for (n3=0; n3<Ngrid1_3; n3++){
sum += CubeData1[n1][n2][n3];
}
}
printf("%5d %18.15f %18.15f\n",n1,step*(double)n1,sum/(double)(n2*n3));
}
}
else if (axis==2){
for (n2=0; n2<Ngrid1_2; n2++){
double sum = 0.0;
for (n1=0; n1<Ngrid1_1; n1++){
for (n3=0; n3<Ngrid1_3; n3++){
sum += CubeData1[n1][n2][n3];
}
}
/* Disabled by N. Yamaguchi ***
printf("%5d %18.15f %18.15f\n",n1,step*(double)n2,sum/(double)(n1*n3));
* ***/
/* Added by N. Yamaguchi ***/
printf("%5d %18.15f %18.15f\n",n2,step*(double)n2,sum/(double)(n1*n3));
/* ***/
}
}
else if (axis==3){
for (n3=0; n3<Ngrid1_3; n3++){
double sum = 0.0;
for (n1=0; n1<Ngrid1_1; n1++){
for (n2=0; n2<Ngrid1_2; n2++){
sum += CubeData1[n1][n2][n3];
}
}
/* Disabled by N. Yamaguchi ***
printf("%5d %18.15f %18.15f\n",n1,step*(double)n3,sum/(double)(n1*n2));
* ***/
/* Added by N. Yamaguchi ***/
printf("%5d %18.15f %18.15f\n",n3,step*(double)n3,sum/(double)(n1*n2));
/* ***/
}
}
/* freeing of arrays */
for (i=0; i<(atomnum1+1); i++){
free(Gxyz1[i]);
}
free(Gxyz1);
for (i=0; i<Ngrid1_1; i++){
for (j=0; j<Ngrid1_2; j++){
free(CubeData1[i][j]);
}
free(CubeData1[i]);
}
free(CubeData1);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment