/************************************************/ /* C2P Mar 19 1999 <liju99@mit.edu> */ /* Convert configuration file of Ju Li' format */ /* to .pdb file which can be viewed by Rasmol. */ /* cc -o c2p c2p.c; c2p -f config.amorphous.864 */ /************************************************/ #include <stdio.h> #define MAX_STRLEN 128 #define TP_STRLEN 2 #define MAX_TP 65536 #define STATS_FILE_HEADER "stats.c2p" struct ATOMTYPE { int occurrence; /* number of atoms of such type */ char name[TP_STRLEN+1]; /* type name */ double mass; /* atomic mass */ }; /* see if a file exists, and if so ask */ /* the user whether to overwrite it. */ int freetowrite (char filename[]) { FILE *file; char c; file = fopen (filename, "r"); if (file != NULL) { fclose (file); printf ("file \"%s\" exists, overwrite (y/n)? ", filename); c = getc(stdin); if (c!='\n') while(getc(stdin)!='\n'); if (!((c=='y')||(c=='Y'))) return(0); } return(1); } /* end freetowrite() */ double matinv(double A[3][3], double B[3][3]) { double D11,D22,D33,D12,D21,D13,D31,D32,D23,determinant; D11=A[1][1]*A[2][2]-A[1][2]*A[2][1]; D22=A[2][2]*A[0][0]-A[2][0]*A[0][2]; D33=A[0][0]*A[1][1]-A[0][1]*A[1][0]; D12=A[1][2]*A[2][0]-A[1][0]*A[2][2]; D23=A[2][0]*A[0][1]-A[2][1]*A[0][0]; D31=A[0][1]*A[1][2]-A[0][2]*A[1][1]; D13=A[1][0]*A[2][1]-A[2][0]*A[1][1]; D21=A[2][1]*A[0][2]-A[0][1]*A[2][2]; D32=A[0][2]*A[1][0]-A[1][2]*A[0][0]; determinant=A[0][0]*D11+A[0][1]*D12+A[0][2]*D13; B[0][0]=D11/determinant; B[1][1]=D22/determinant; B[2][2]=D33/determinant; B[0][1]=D21/determinant; B[1][2]=D32/determinant; B[2][0]=D13/determinant; B[1][0]=D12/determinant; B[2][1]=D23/determinant; B[0][2]=D31/determinant; return (determinant); } /* end matinv() */ #define max(x,y) ((x)>(y)?(x):(y)) int main (int argc, char *argv[]) { int i, j, k, np, nc[3], ntp, to_stdout = 0, force_write = 0, save_stats = 0; char buf[MAX_STRLEN], read[MAX_STRLEN], write[2*MAX_STRLEN], tp[TP_STRLEN+1], tps[MAX_TP][TP_STRLEN+1], *r, *p, *q; struct ATOMTYPE at[MAX_TP]; double H[3][3], HI[3][3], mass, s[3], x[3]; FILE *in, *out, *stats; for (k=1; k<max(argc,2); k++) { if (argc == 1) { /* no input argument */ printf ("Convert config file of JL's format: "); for (r=read; r<read+MAX_STRLEN-1; r++) if ((*r=(char)getc(stdin))=='\n') break; *r = (char)0; } else if (argv[k][0] == '-') { /* command line options */ for (p=argv[k]+1; *p!=(char)0; p++) switch (*p) { case 'f': /* forced write like % mv -f */ force_write = 1; break; case 'o': to_stdout = 1; break; case 's': save_stats = 1; break; case 'h': printf ("-------------------------------------------\n" "Convert configuration file of Ju Li' format\n" "to .pdb file which can be viewed by Rasmol.\n" "-------------------------------------------\n" "%% c2p -f config.amorphous config.xtal\n" "-f: forced overwrite existent .pdb file\n" "-o: send converted .pdb to stdout\n" "-s: save stats to \"%s.*\" instead of stdout\n", STATS_FILE_HEADER); break; default: printf ("c2p: unknown option \"-%c\".\n", *p); exit(1); } continue; } else /* configuration filename to be processed */ sprintf (read, "%s", argv[k]); /* my handling of string input is more */ /* robust than scanf(): */ for (r=read; (*r==' ')||(*r=='\t'); r++); for (q=r; *q!=(char)0; q++) if ((*q==' ')||(*q=='\t')) *q=(char)0; /* open that configuration file */ if ( (in=fopen(r,"r")) == NULL ) { printf("c2p: cannot open config file \"%s\".\n", r); continue; } fscanf (in, "Number of particles = %d\n\n", &np); fscanf (in, "H(1,1) = %lf A\n", &H[0][0]); fscanf (in, "H(1,2) = %lf A\n", &H[0][1]); fscanf (in, "H(1,3) = %lf A\n", &H[0][2]); fscanf (in, "nc(1) = %d\n\n", &nc[0]); fscanf (in, "H(2,1) = %lf A\n", &H[1][0]); fscanf (in, "H(2,2) = %lf A\n", &H[1][1]); fscanf (in, "H(2,3) = %lf A\n", &H[1][2]); fscanf (in, "nc(2) = %d\n\n", &nc[1]); fscanf (in, "H(3,1) = %lf A\n", &H[2][0]); fscanf (in, "H(3,2) = %lf A\n", &H[2][1]); fscanf (in, "H(3,3) = %lf A\n", &H[2][2]); fscanf (in, "nc(3) = %d\n\n", &nc[2]); /* descriptor line */ fscanf (in, "TP Mass(amu) sx sy sz"); /* absorb rest of the stuff on this line */ while (((buf[0]=(char)getc(in))!='\n')&&(buf[0]!=(char)EOF)); if ( !to_stdout ) { sprintf (write, "%s.pdb", r); if (argc == 1) { /* interactive mode */ printf("Write to (default=\"%s\"): ", write); q = write + strlen(write); for (p=q; p<write+2*MAX_STRLEN-5; p++) if ((*p=(char)getc(stdin))=='\n') break; if (p==q) q = write; *p = (char)0; /* add file suffix if without */ if (strcasecmp(p-4,".pdb")) sprintf(p,".pdb"); } else q = write; for (; (*q==' ')||(*q=='\t'); q++); if ( force_write || freetowrite(q) ) out = fopen (q, "w"); else continue; } else out = stdout; /* .pdb header: */ fprintf (out, "HEADER Converted from \"%s\" by c2p.\n", r); /* H matrix in footnote -1 and 0: */ fprintf (out, "FTNOTE 1 Supercell Parallelepiped (H) in Angstroms:\n" "FTNOTE 2 a=%d*[%.3f %.3f %.3f]=(%.3f %.3f %.3f)\n" "FTNOTE 3 b=%d*[%.3f %.3f %.3f]=(%.3f %.3f %.3f)\n" "FTNOTE 4 c=%d*[%.3f %.3f %.3f]=(%.3f %.3f %.3f)\n", nc[0], H[0][0]/nc[0], H[0][1]/nc[0], H[0][2]/nc[0], H[0][0], H[0][1], H[0][2], nc[1], H[1][0]/nc[1], H[1][1]/nc[1], H[1][2]/nc[1], H[1][0], H[1][1], H[1][2], nc[2], H[2][0]/nc[2], H[2][1]/nc[2], H[2][2]/nc[2], H[2][0], H[2][1], H[2][2]); /* crystal dimensions: */ fprintf (out, "CRYST1 %8.3f %8.3f %8.3f " "90.000 90.000 90.000 P 1 1\n", H[0][0], H[1][1], H[2][2]); for (ntp=i=0; i<np; i++) { /* atom types and positions: */ fscanf (in, "%2s %lf %lf %lf %lf", tp, &mass, &s[0], &s[1], &s[2]); for (j=0; j<ntp; j++) if ( (!strcmp(tp, at[j].name)) && (mass==at[j].mass) ) { at[j].occurrence++; break; } if (j==ntp) /* found a new atom type */ if (ntp<MAX_TP) { sprintf(at[ntp].name, "%2s", tp); at[ntp].mass = mass; at[ntp++].occurrence = 1; } else { printf ("v: MAX_TP = %d exceeded.\n", MAX_TP); exit(1); } while (((buf[0]=(char)getc(in))!='\n')&&(buf[0]!=(char)EOF)); if (buf[0] == (char)EOF) { printf ("c2p: premature end of config file \"%s\":", r); printf ("claimed number of particles = %d, actual = %d\n", np, i+1); goto finish; } x[0] = s[0]*H[0][0]+s[1]*H[1][0]+s[2]*H[2][0]; x[1] = s[0]*H[0][1]+s[1]*H[1][1]+s[2]*H[2][1]; x[2] = s[0]*H[0][2]+s[1]*H[1][2]+s[2]*H[2][2]; fprintf (out, "ATOM %5d %2s %2s " "%8.3f%8.3f%8.3f\n", i+1, tp, " 1", x[0], x[1], x[2]); } finish: /* print statistics */ close(in); if ( save_stats ) { /* look for true file name without "/" */ for (p=r+strlen(r)-1; p>=r; p--) if (*p=='/') break; sprintf (buf, "%s.%s", STATS_FILE_HEADER, p+1); stats = fopen (buf, "w"); } else stats = stdout; if ( (!to_stdout) || save_stats ) { fprintf (stats, "\nIn \"%s\", we found\n", r); fprintf (stats, "--------------------------------------\n"); fprintf (stats, " Occurrences Type Mass\n"); for (j=0; j<ntp; j++) fprintf (stats, " %7d %2s %.3f\n", at[j].occurrence, at[j].name, at[j].mass); fprintf (stats, "--------------------------------------\n"); if (!to_stdout) fprintf (stats, "saved to \"%s\".\n\n", q); close(stats); if (out!=NULL) close(out); } } /* loop over files */ return(0); } /* end main() */