/*
************************************************************************
*           All Copyright Reserved!                                    *
*                                                                      *
*  Prog:    respgen                                                    *
*  Version: version 1.0                                                *
*  Author:  Junmei Wang                                                *
*                                                                      *
*  Department of Pharmaceutical Chemistry                              *
*  School of Pharmacy                                                  *
*  University of California                                            *
*  San Francisco   CA 94143                                            *
*  Octomber, 2001                                                      *
************************************************************************
*/
# include "common.h"
# include "define.h"
# include "atom.h"
# include "utility.c"
# include "common.c"
# include "ac.c"
# define MAXCONF 20
# define MAXCHARGE 20

typedef struct {
	char name[20];
	double charge;
	int id;
} CHARGE; 

typedef struct {
	char resat[20];
	char capat[20];
	int  resid;
	int  capid;
} SEPBOND; 

ATOM *atom;
BOND *bond;
int atomnum = 0;
int bondnum = 0;
int *flag;
CONTROLINFO cinfo;
MOLINFO minfo;
int i, j, k;
FILE *fpin;
FILE *fpout;
char line[MAXCHAR];
char ifilename[MAXCHAR];

int overflow_flag = 0;

int confnum = 1;
SEPBOND sepbond[2];
int sepbondnum;
CHARGE charge[MAXCHARGE];
int chargenum = 0;
int residueatnum ;
double netcharge = 0;
char inputfile[MAXCHAR];
char espfile[MAXCHAR];
char prepfile[MAXCHAR];
char residue_symbol[MAXCHAR] = "MOL";
char residue_file_name[MAXCHAR] = "MOL.res";

int *selectelement;
int selectnum ;
int selectnum2;

void findpath(ATOM atm[], int selectnum, int startnum) {
int i, j;
int start;
int resetindex = -1;
start = -1;
selectelement[selectnum++] = startnum;

for (i = 0; i < 6; i++) {
        if(atm[startnum].con[i] == -1) return;
	if(atm[startnum].con[i] == sepbond[1].resid) {
		selectelement[selectnum++] = sepbond[1].resid;
		selectnum2 = selectnum ;
		return;
	}
        start = atm[startnum].con[i];
        for (j = 0; j < selectnum; j++)
                if (start == selectelement[j]) {
                        resetindex = 1;
                        break;
                }
	if(resetindex == 1) {
		resetindex = -1;
		continue;
	}
        if (start == -1) return;
        findpath(atm, selectnum, start);
}
}

void group_atom(int index) {
        int i, j, num, num_old;
        int id;
        num = 1;
        num_old = 0;

	if(index == 0) {
        	flag[sepbond[0].resid] = 1;
        	flag[sepbond[1].resid] = 1;
	}
	else if(index == 1) 
        	flag[sepbond[0].capid] = 1;
	else if(index == 2) 
        	flag[sepbond[1].capid] = 1;
		

        while (num_old < num) {
                num_old = num;
                for (i = 0; i<atomnum; i++)
                        if(flag [i] != -1)
                        for(j =0 ; j<6; j++){
                                id = atom[i].con[j];
                                if(id == -1) continue;
                                if(flag[id] != -1) continue;
				if(index == 0) {
                                	if(id == sepbond[0].capid) continue;
                                	if(id == sepbond[1].capid) continue;
				}
				else if (index == 1) {
                                	if(id == sepbond[0].resid) continue;
				}
				else if (index == 2) {
                                	if(id == sepbond[1].resid) continue;
				}
                                flag[id] = 1;
                                num ++;
                        }
                }
}

void proceed1() {
int i;
int length;
int index;
char command[MAXCHAR];
char tmpchar[MAXCHAR];
	printf("\n----- This is an amino-acid-like or nucleic-acid-like residue ---\n");

	selectelement = (int *) malloc(sizeof(int) * atomnum);
	if (selectelement == NULL) {
       		fprintf(stderr, "memory allocation error for *selectelement\n");
       		exit(1);
	}

	flag = (int *) malloc(sizeof(int) * atomnum);
	if (flag == NULL) {
        	fprintf(stderr, "memory allocation error for *flag\n");
        	exit(1);
	}
	for(i =0; i< chargenum; i++) { 
		charge[i].id = -1;
		index = 0;
		length = strlen(charge[i].name);
		for(j =0; j< atomnum; j++) 
			if(strncmp(atom[j].name, charge[i].name,length) == 0) {
				charge[i].id = j;
				index = 1;
				break;
			}
		if(index == 0) {
			fprintf(stderr, "The charge name, %s does not show up in the molecule\n", charge[i].name);
			exit(1);
		}
	}
	sepbond[0].resid = -1;
	sepbond[1].resid = -1;
	sepbond[0].capid = -1;
	sepbond[1].capid = -1;
	for(i =0; i< atomnum; i++) {
		if(sepbond[0].resid != -1 && sepbond[0].capid != -1 &&
		   sepbond[1].resid != -1 && sepbond[1].capid != -1) break;
		length = 0;
		length = strlen(sepbond[0].resat);			
		if(strncmp(atom[i].name, sepbond[0].resat, length) == 0) {
			sepbond[0].resid = i; 
			continue;
		}
		length = strlen(sepbond[0].capat);			
		if(strncmp(atom[i].name, sepbond[0].capat, length) == 0) { 
			sepbond[0].capid = i; 
			continue;
		}
		length = strlen(sepbond[1].resat);			
		if(strncmp(atom[i].name, sepbond[1].resat, length) == 0) {
			sepbond[1].resid = i; 
			continue;
		}
		length = strlen(sepbond[1].capat);			
		if(strncmp(atom[i].name, sepbond[1].capat, length) == 0) {
			sepbond[1].capid = i; 
			continue;
		}
	}	

	printf("\n\nStep 1:  Find residue atoms and cap atoms ... ");
/*	for the first cap	*/
	printf("\n\nFirst Cap Atoms");
        for (i = 0; i<atomnum; i++)
                flag[i] = -1;
	group_atom(1);
        for (i = 0; i<atomnum; i++)
                if(flag[i] == 1) 
			printf("\nATOM\t %d\t%s", i+1, atom[i].name);

/*	for the second cap	*/
	printf("\n\nSecond Cap Atoms");
        for (i = 0; i<atomnum; i++)
                flag[i] = -1;
	group_atom(2);
        for (i = 0; i<atomnum; i++)
                if(flag[i] == 1) 
			printf("\nATOM\t %d\t%s", i+1, atom[i].name);

/*	for residue	*/
	printf("\n\nResidue Atoms");
        for (i = 0; i<atomnum; i++)
                flag[i] = -1;
	group_atom(0);
	residueatnum = 0;
        for (i = 0; i<atomnum; i++)
                if(flag[i] == 1) {
			residueatnum ++;
			printf("\nATOM\t %d\t%s", i+1, atom[i].name);
		}

/*	for main chain	*/
	selectnum = 0;
        for (i = 0; i<atomnum; i++)
                selectelement[i] = -1;
	findpath(atom, 0, sepbond[0].resid); 
	printf("\n\nMain chain Atoms");
        	for (j = 0; j < selectnum2 ; j++)
			printf("\nATOM\t%d\t%s", j+1, atom[selectelement[j]].name);

	printf("\n\nStep 2:  Generate RESP input files ... \n");
	printf("\nGenerate charge input file for respgen ... \n");
       	if ((fpout = fopen("RESIDUE_GEN_RESPGEN.DAT", "w")) == NULL) {
               	fprintf(stderr, "Cannot open input file RESIDUE_GEN_RESPGEN.DAT, exit\n");
               	exit(1);
       	}
	fprintf(fpout, "//predefined charges in a format of (CHARGE partical_charge atom_ID atom_name)");
	for(i=0;i<chargenum;i++)
		fprintf(fpout, "\nCHARGE	%9.6lf\t%d\t%s", charge[i].charge, charge[i].id + 1, charge[i].name);		

	fprintf(fpout, "\n//charge groups in a format of (GROUP num_atom net_charge)");		
	fprintf(fpout, "\nGROUP\t%d\t%9.5lf\n", residueatnum, netcharge);		
	fprintf(fpout, "//atoms in the group in a format of (ATOM atom_ID atom_name)");		
        for (i = 0; i<atomnum; i++)
                if(flag[i] == 1) 
			fprintf(fpout, "\nATOM\t%d\t%s", i+1, atom[i].name);
	fclose(fpout);

	printf("\nRun respgen to generate first stage resp input file ... \n");
	strcpy(command, BINDIR);
	strcat(command, "/respgen -i ");
	strcat(command, inputfile);
	sprintf(tmpchar, "%d", confnum);
	strcat(command, " -o RESIDUE_GEN_RESP.INPUT1 -f resp1 -a RESIDUE_GEN_RESPGEN.DAT -n ");
	strcat(command, tmpchar);
	printf("\nCommand: %s\n", command);
	system(command);

	printf("\nRun respgen to generate first stage resp input file ... \n");
	strcpy(command, BINDIR);
	strcat(command, "/respgen -i ");
	strcat(command, inputfile);
	sprintf(tmpchar, "%d", confnum);
	strcat(command, " -o RESIDUE_GEN_RESP.INPUT2 -f resp2 -a RESIDUE_GEN_RESPGEN.DAT -n ");
	strcat(command, tmpchar);
	printf("\nCommand: %s\n", command);
	system(command);

	strcpy(command, BINDIR);
	strcat(command, "/resp -O -i RESIDUE_GEN_RESP.INPUT1 -o RESIDUE_GEN_RESP.OUTPUT1 -q QIN -e ");
	strcat(command, espfile);
	printf("\nCommand: %s\n", command);
	system(command);

	strcpy(command, BINDIR);
	strcat(command, "/resp -O -i RESIDUE_GEN_RESP.INPUT2 -o RESIDUE_GEN_RESP.OUTPUT2 -q qout -e ");
	strcat(command, espfile);
	printf("\nCommand: %s\n", command);
	system(command);

	printf("\n\nStep 3:  Read RESP charges ... \n");
	strcpy(command, BINDIR);
	strcat(command, "/antechamber -fi ac -fo ac -i ");
	strcat(command, inputfile);
	strcat(command, " -o RESIDUE_GEN.AC -c rc -cf qout -at amber");
	printf("\nCommand: %s\n", command);
	system(command);

	printf("\n\nStep 4:  Generate mainchain file to be read by prepgen ... \n");

        if ((fpout = fopen("RESIDUE_GEN_MAINCHAIN.DAT", "w")) == NULL) {
                fprintf(stderr, "Cannot open input file RESIDUE_GEN_MAINCHAIN.DAT, exit\n");
                exit(1);
        }
	fprintf(fpout, "HEAD_NAME\t%s\n", sepbond[0].resat);
	fprintf(fpout, "TAIL_NAME\t%s\n", sepbond[1].resat);
       	for (i = 0; i < selectnum2 ; i++) {
		length = strlen(sepbond[0].resat);			
		if(strncmp(atom[selectelement[i]].name, sepbond[0].resat, length) == 0) continue;
		length = strlen(sepbond[1].resat);			
		if(strncmp(atom[selectelement[i]].name, sepbond[1].resat, length) == 0) continue;
		fprintf(fpout, "MAIN_CHAIN\t%s\n", atom[selectelement[i]].name);
	}
        for (i = 0; i<atomnum; i++) 
                if(flag[i] == 1) 
			continue; 
		else
			fprintf(fpout, "OMIT_NAME\t%s\n", atom[i].name);
	fprintf(fpout, "PRE_HEAD_TYPE\t%s\n", atom[sepbond[0].capid].ambername);
	fprintf(fpout, "POST_TAIL_TYPE\t%s\n", atom[sepbond[1].capid].ambername);
	fprintf(fpout, "CHARGE\t%9.5lf\n", netcharge);
	fclose(fpout);

	printf("\n\nStep 5:  Generate prep input file ... \n");
	strcpy(command, BINDIR);
	strcat(command, "/prepgen -i RESIDUE_GEN.AC -o ");
	strcat(command, prepfile);
	strcat(command, " -m RESIDUE_GEN_MAINCHAIN.DAT -rn "); 
	strcat(command, residue_file_name);
	strcat(command, " -rf "); 
	strcat(command, residue_symbol);
	printf("\nCommand: %s\n", command);
	system(command);
	printf("\n\nAll Done !\n");
}

void proceed2() {
int index;
int length;
char command[MAXCHAR];
char tmpchar[MAXCHAR];
	printf("\n----- This is a normal small molecule ---\n");
	for(i =0; i< chargenum; i++) { 
		charge[i].id = -1;
		index = 0;
		length = strlen(charge[i].name);
		for(j =0; j< atomnum; j++) 
			if(strncmp(atom[j].name, charge[i].name,length) == 0) {
				charge[i].id = j;
				index = 1;
				break;
			}
		if(index == 0) {
			fprintf(stderr, "The charge name, %s does not show up in the molecule\n", charge[i].name);
			exit(1);
		}
	}

	printf("\n\nStep 1:  Generate RESP input files ... \n");
	printf("\nGenerate charge input file for respgen ... \n");
       	if ((fpout = fopen("RESIDUE_GEN_RESPGEN.DAT", "w")) == NULL) {
               	fprintf(stderr, "Cannot open input file RESIDUE_GEN_RESPGEN.DAT, exit\n");
               	exit(1);
       	}
        fprintf(fpout, "//predefined charges in a format of (CHARGE partical_charge atom_ID atom_name)");
        for(i=0;i<chargenum;i++)
                fprintf(fpout, "\nCHARGE        %9.6lf\t%d\t%s", charge[i].charge, charge[i].id + 1, charge[i].name);
        fclose(fpout);

	printf("\nRun respgen to generate first stage resp input file ... \n");
	strcpy(command, BINDIR);
	strcat(command, "/respgen -i ");
	strcat(command, inputfile);
	sprintf(tmpchar, "%d", confnum);
	strcat(command, " -o RESIDUE_GEN_RESP.INPUT1 -f resp1 -a RESIDUE_GEN_RESPGEN.DAT -n ");
	strcat(command, tmpchar);
	printf("\nCommand: %s\n", command);
	system(command);

	printf("\nRun respgen to generate first stage resp input file ... \n");
	strcpy(command, BINDIR);
	strcat(command, "/respgen -i ");
	strcat(command, inputfile);
	sprintf(tmpchar, "%d", confnum);
	strcat(command, " -o RESIDUE_GEN_RESP.INPUT2 -f resp2 -a RESIDUE_GEN_RESPGEN.DAT -n ");
	strcat(command, tmpchar);
	printf("\nCommand: %s\n", command);
	system(command);

	strcpy(command, BINDIR);
	strcat(command, "/resp -O -i RESIDUE_GEN_RESP.INPUT1 -o RESIDUE_GEN_RESP.OUTPUT1 -q QIN -e ");
	strcat(command, espfile);
	printf("\nCommand: %s\n", command);
	system(command);

	strcpy(command, BINDIR);
	strcat(command, "/resp -O -i RESIDUE_GEN_RESP.INPUT2 -o RESIDUE_GEN_RESP.OUTPUT2 -q qout -e ");
	strcat(command, espfile);
	printf("\nCommand: %s\n", command);
	system(command);

	printf("\n\nStep 2:  Read RESP charges ... \n");
	strcpy(command, BINDIR);
	strcat(command, "/antechamber -fi ac -fo ac -i ");
	strcat(command, inputfile);
	strcat(command, " -o RESIDUE_GEN.AC -c rc -cf qout -at amber");
	printf("\nCommand: %s\n", command);
	system(command);

	printf("\n\nStep 3:  Generate prep input file ... \n");
	strcpy(command, BINDIR);
	strcat(command, "/prepgen -i RESIDUE_GEN.AC -o ");
	strcat(command, prepfile);
	strcat(command, " -rn "); 
	strcat(command, residue_file_name);
	strcat(command, " -rf "); 
	strcat(command, residue_symbol);
	printf("\nCommand: %s\n", command);
	system(command);
	printf("\n\nAll Done !\n");
}

int main(int argc, char *argv[])
{
	if (strcmp(COLORTEXT, "YES") == 0 || strcmp(COLORTEXT, "yes") == 0) {
		if (argc == 2
			&& (strcmp(argv[1], "-h") == 0
				|| strcmp(argv[1], "-H") == 0)) {
			printf("Usage: residuegen  input_file \n");
			exit(1);
		}
		if (argc != 2) {
			printf("Usage: residuegen  input_file \n");
			exit(1);
		}
	} else {
		if (argc == 2
			&& (strcmp(argv[1], "-h") == 0
				|| strcmp(argv[1], "-H") == 0)) {
			printf("Usage residuegen  input_file\n");
			exit(1);
		}
		if (argc != 2) {
			printf("Usage residuegen  input_file\n");
			exit(1);
		}
	}
        if ((fpin = fopen(argv[1], "r")) == NULL) {
                fprintf(stderr, "Cannot open input file: %s, exit\n", argv[1]);
                exit(1);
        }
	
	sepbondnum = 0;
        for (;;) {
                if (fgets(line, MAXCHAR, fpin) == NULL) break;
                if (strncmp("ESP_FILE", line, 8) == 0) {
			sscanf(&line[9], "%s", espfile); 
                        continue;
                }
                if (strncmp("INPUT_FILE", line, 10) == 0) {
			sscanf(&line[11], "%s", inputfile); 
                        continue;
                }
                if (strncmp("CONF_NUM", line, 8) == 0) {
			sscanf(&line[9], "%d", &confnum); 
                        continue;
                }
                if (strncmp("SEP_BOND", line, 8) == 0) {
			sscanf(&line[9], "%s%s", sepbond[sepbondnum].resat, sepbond[sepbondnum].capat); 
			sepbondnum++;
			if(sepbondnum > 2) {
				fprintf(stderr, "\nToo many separative bonds, exit"); 
				exit(1);
			}
                        continue;
                }
                if (strncmp("NET_CHARGE", line, 10) == 0) {
			sscanf(&line[11], "%lf", &netcharge); 
                        continue;
                }
                if (strncmp("PREP_FILE", line, 9) == 0) {
			sscanf(&line[10], "%s", prepfile); 
                        continue;
                }
                if (strncmp("RESIDUE_FILE_NAME", line, 17) == 0) {
			sscanf(&line[18], "%s", residue_file_name); 
                        continue;
                }
                if (strncmp("RESIDUE_SYMBOL", line, 14) == 0) {
			sscanf(&line[15], "%s", residue_symbol); 
                        continue;
                }
                if (strncmp("ATOM_CHARGE", line, 11) == 0) {
			sscanf(&line[12], "%s%lf", charge[chargenum].name, &charge[chargenum].charge); 
			chargenum++;
                        continue;
                }
	}
	fclose(fpin);
	if(sepbondnum != 2 && sepbondnum != 0) {
		printf("There must be two separative bonds or no separative bond at all");
		exit(1);
	}
	printf("\nINPUT_FILE:\t\t%s", inputfile);
	printf("\nESP_FILE:\t\t%s", espfile);
	printf("\nCONF_NUM:\t\t%d", confnum);
	printf("\nNET_CHARGE:\t\t%-9.5lf", netcharge);
	for(i=0;i<sepbondnum;i++)
		printf("\nSEPBOND\t\t\t%d\t%s\t%s", i+1, sepbond[i].resat, sepbond[i].capat);
	printf("\nThere are %d predefined charges", chargenum);
	for(i=0;i<chargenum;i++) 
			printf("\nATOM_CHARGE\t\t%d\t%s\t%9.5lf", i+1, charge[i].name, charge[i].charge);
	printf("\nPREP_FILE: %s", prepfile);

	default_minfo(&minfo);
	default_cinfo(&cinfo);

	atom = (ATOM *) malloc(sizeof(ATOM) * cinfo.maxatom);
	if (atom == NULL) {
		fprintf(stderr, "memory allocation error for *atom\n");
		exit(1);
	}

	bond = (BOND *) malloc(sizeof(BOND) * cinfo.maxbond);

	if (bond == NULL) {
		fprintf(stderr, "memory allocation error for *bond\n");
		exit(1);
	}

	overflow_flag =
		rac(inputfile, &atomnum, atom, &bondnum, bond, &cinfo, &minfo);

	if (overflow_flag) {
		cinfo.maxatom = atomnum + 10;
		cinfo.maxbond = bondnum + 10;
		free(atom);
		free(bond);
		atom = (ATOM *) malloc(sizeof(ATOM) * cinfo.maxatom);
		if (atom == NULL) {
			fprintf(stderr, "memory allocation error for *atom\n");
			exit(1);
		}
		bond = (BOND *) malloc(sizeof(BOND) * cinfo.maxbond);
		if (bond == NULL) {
			fprintf(stderr, "memory allocation error for *bond\n");
			exit(1);
		}
		overflow_flag =
			rac(ifilename, &atomnum, atom, &bondnum, bond, &cinfo, &minfo);
	}
	atomicnum(atomnum, atom);
	if(sepbondnum == 2) proceed1();
	if(sepbondnum == 0) proceed2();
	printf("\n");
	return (0);

}