/************************************************/
/* 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() */