#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <time.h>
#include <sys/types.h>

#define  NMISC 29 
extern FILE *fout;

#include "nrutil.h"
#include "util.h"
#include "match.h"



long number_of_atoms(char *pdbfile)
/* number of atom records in a PDB file */
{
    char str[512];
    long n =0, nlen=0;
    FILE *fp;

    if((fp = fopen(pdbfile, "r"))==NULL){
        printf("Can not open the file %s (routine: number_of_atoms)\n",
               pdbfile);
        return 0;
    }
    while (fgets(str, sizeof str, fp) != NULL) {
        nlen = upperstr(str);
        if (!strncmp(str+17, "HOH", 3) || !strncmp(str+17, "WAT", 3)
            || !strncmp(str+17, "SOL", 3) || !strncmp(str+17, "WO2", 3) )
            continue; /*get ride of WATER */  
        if (str[13] == 'H')
            continue;  /*get ride of H atoms */
		
        if (nlen >= 54 && !strncmp(str, "ATOM", 4) 
            || !strncmp(str, "HETATM", 6)  )
            n++;
    }
    fclose(fp);
    return n;
}

long read_pdb(char *pdbfile, char **AtomName, char **ResName, char *ChainID,
              long *ResSeq, double **xyz, char **Miscs, char *ALT_LIST)
/* read in a PDB file and do some processing
 * Miscs[][NMISC]: H/A, altLoc, iCode, occ./tempFac./segID/element/charge
 *           col#   0     1       2       3-28 [combined together]
 */
{
    char str[512], temp[512], *pchar, str_id[20], str_id0[20];
    long i=0, n=0, nlen=0;
    char atomname[5], resname[4], chainid, resseq[5];
    FILE *fp;
          
    if((fp = fopen(pdbfile, "r"))==NULL) {        
        printf("Can not open the file %s (routine: read_pdb)\n",pdbfile);
        return 0;
    }
       
    n=1;
    while (fgets(str, sizeof str, fp) != NULL) {
        nlen = upperstr(str);
        if (!strncmp(str+17, "HOH", 3) || !strncmp(str+17, "WAT", 3)
            || !strncmp(str+17, "SOL", 3) || !strncmp(str+17, "WO2", 3))
            continue; /*get ride of WATER */
        
        if (str[13] == 'H')
            continue;  /*get ride of H atoms */
        
        if (nlen >= 54 && (!strncmp(str, "ATOM", 4) 
            || !strncmp(str, "HETATM", 6) ) ){

            strncpy(atomname, str + 12, 4);
            atomname[4] = '\0';

            strncpy(resname, str + 17, 3);        /* residue name */
            resname[3] = '\0';
            /* delete ending spaces as in "C  " */
            for (i = 1; i <= 2; i++)
                if (resname[2] == ' ') {
                    resname[2] = resname[1];
                    resname[1] = resname[0];
                    resname[0] = ' ';
                }
            if (resname[2] == ' ') {
                printf( "%s\n", str);
                nrerror("==> residue name field empty <==");
            }

            chainid=str[21];

            strncpy(resseq, str + 22, 4);    /* residue sequence */
            resseq[4] = '\0';

            sprintf(str_id, "%s%s%c%s",atomname, resname, chainid, resseq); 


            if(!strcmp(str_id, str_id0) && n>1) continue; /*rid of alternate*/

            strcpy(AtomName[n], atomname);
            strcpy(ResName[n], resname);
            ChainID[n] = chainid;
            if (sscanf(resseq, "%4ld", &ResSeq[n]) != 1) {
                printf( "residue #? ==> %.54s\n", str);
                ResSeq[n] = 9999;
            }

            strncpy(temp, str + 30, 25);           /* xyz */
            temp[25] = '\0';
            if (sscanf(temp,"%8lf%8lf%8lf",
                       &xyz[n][1],&xyz[n][2],&xyz[n][3])!=3)
                nrerror("error reading xyz-coordinate");
            
            Miscs[n][0] = str[0];        /* H for HETATM, A for ATOM */
            Miscs[n][1] = str[16];        /* alternative location indicator */
            Miscs[n][2] = str[26];        /* code of insertion residues */
            strncpy(Miscs[n] + 3, str + 54, NMISC - 3);
            if ((pchar = strrchr(Miscs[n], '\n')) != NULL)
                Miscs[n][pchar - Miscs[n]] = '\0';
            else
                Miscs[n][NMISC] = '\0';                /* just to make sure */
			            
			
            if (AtomName[n][3] == '*')        /* * to ' */
                AtomName[n][3] = '\'';
            if (!strcmp(AtomName[n], " O1'")){        /* O1' to O4' */
                strcpy(AtomName[n], " O4'");
            }else if (!strcmp(AtomName[n], " OL ")) {       /* OL to O1P */
                strcpy(AtomName[n], " O1P");
            }else if (!strcmp(AtomName[n], " OR ")) {       /* OR to O2P */
                strcpy(AtomName[n], " O2P");
            }else if  (!strcmp(AtomName[n], " C5A")){        /* C5A to C5M */
                strcpy(AtomName[n], " C5M");
            }else if (!strcmp(AtomName[n], " O5T")){        /* terminal O5' */
                strcpy(AtomName[n], " O5'");
            }else if (!strcmp(AtomName[n], " O3T")){        /* terminal O3' */
                strcpy(AtomName[n], " O3'");
            }
/*            
         printf("%s%5ld %4s%c%3s %c%4ld%c   %8.3lf%8.3lf%8.3lf\n", 
                 "ATOM  ", n, AtomName[n], Miscs[n][1],
                ResName[n], ChainID[n], ResSeq[n], Miscs[n][2], xyz[n][1],
                xyz[n][2], xyz[n][3]);
            printf("WHAT: %6d  atom=%s; res=%s; id=%c; nseq=%s\n",n,
                   AtomName[n],ResName[n] ,chainid, resseq);
*/
           

            n++;
            sprintf(str_id0, "%s%s%c%s",atomname, resname, chainid, resseq); 
        }
    }

    fclose(fp);
    return n-1;
}

long **residue_idx(long num, long *ResSeq, char **Miscs, char *ChainID,
                   char **ResName, long *num_residue)
/* number of residues, and starting-ending indexes for each */
{
    char iCode;
    char **bidx;
    long i=0, n=0, **seidx, *temp;

    bidx = cmatrix(1, num, 0, 12);        /* normally 9 */
    temp = lvector(1, num);

    for (i = 1; i <= num; i++) {
        iCode = (Miscs == NULL) ? ' ' : Miscs[i][2];
        sprintf(bidx[i],"%3s%c%4ld%c",ResName[i],ChainID[i],ResSeq[i], iCode);
    }
    for (i = 1; i < num; i++)
        temp[i] = strcmp(bidx[i + 1], bidx[i]) ? 1 : 0;
    temp[num] = 1;

    n = 0;                        /* get number of residues */
    for (i = 1; i <= num; i++)
        if (temp[i])
            ++n;

    seidx = lmatrix(0, n, 0, 2);        /* allocate spaces */
    n = 0;
    for (i = 1; i <= num; i++)
        if (temp[i])
            seidx[++n][2] = i;
    for (i = 2; i <= n; i++)
        seidx[i][1] = seidx[i - 1][2] + 1;
    seidx[1][1] = 1;

    *num_residue = n;

    free_cmatrix(bidx, 1, num, 0, 12);
    free_lvector(temp, 1, num);
/*    printf("!!num_residue= %5d %5d \n",*num_residue, num);*/

    return seidx;
}


void get_chain_idx(long num_residue, long **seidx, char *ChainID, long *RY,
                   long *nchain, long **chain_idx)
/* get chain index */
{

    long n=0, j=0, k1=0, k2=0;
    
    chain_idx[1][1] = 1;
    n=1;
    for (j=2; j<=num_residue; j++){
        k1 = seidx[j-1][1];
        k2 = seidx[j][1];
        chain_idx[n][2] = j-1;
        if((ChainID[k1] != ChainID[k2])){
            n++;
            chain_idx[n][1] = j;
        }
    }
    chain_idx[n][2] = num_residue;
    *nchain = n;
}


void get_one_lett_seq(long num_residue, long **seidx, char **AtomName,
                      char **ResName, char *ChainID, long *ResSeq, 
                      char **Miscs,double **xyz, char *bseq, long *RY)
/* Get one letter sequnce. The modified give X or ?  */
{
    
    long i=0, ib=0, ie=0;
    
    for (i = 1; i <= num_residue; i++) {
        ib = seidx[i][1];
        ie = seidx[i][2];
        RY[i] = residue_ident(AtomName, xyz, ib, ie);
        if (RY[i] >= 0) { /* nuclear acids */
            
            if      (!strcmp(ResName[ib],"  A") || !strcmp(ResName[ib],"ADE"))
                bseq[i] = 'A';
            else if (!strcmp(ResName[ib],"  G") || !strcmp(ResName[ib],"GUA"))
                bseq[i] = 'G';
            else if (!strcmp(ResName[ib],"  U") || !strcmp(ResName[ib],"URA"))
                bseq[i] = 'U';
            else if (!strcmp(ResName[ib],"  C") || !strcmp(ResName[ib],"CYT"))
                bseq[i] = 'C';
            else if (!strcmp(ResName[ib],"  T") || !strcmp(ResName[ib],"THY"))
                bseq[i] = 'T';

            else 
                bseq[i] = 'X';

        }else if (RY[i] == -1){ /* amino acid */
            
/* Hydrophobic  amino acids */
            if      (!strcmp(ResName[ib],"  L") || !strcmp(ResName[ib],"LEU"))
                bseq[i] = 'L';
            else if (!strcmp(ResName[ib],"  I") || !strcmp(ResName[ib],"ILE"))
                bseq[i] = 'I';
            else if (!strcmp(ResName[ib],"  A") || !strcmp(ResName[ib],"ALA"))
                bseq[i] = 'A';
            else if (!strcmp(ResName[ib],"  V") || !strcmp(ResName[ib],"VAL"))
                bseq[i] = 'V';
            else if (!strcmp(ResName[ib],"  F") || !strcmp(ResName[ib],"PHE"))
               bseq[i] = 'F';
            else if (!strcmp(ResName[ib],"  P") || !strcmp(ResName[ib],"PRO"))
               bseq[i] = 'P';
            else if (!strcmp(ResName[ib],"  M") || !strcmp(ResName[ib],"MET"))
               bseq[i] = 'M';

/* chared amino acids */
            else if (!strcmp(ResName[ib],"  D") || !strcmp(ResName[ib],"ASP"))
               bseq[i] = 'D';
            else if (!strcmp(ResName[ib],"  E") || !strcmp(ResName[ib],"GLU"))
               bseq[i] = 'E';
            else if (!strcmp(ResName[ib],"  K") || !strcmp(ResName[ib],"LYS"))
               bseq[i] = 'K';
            else if (!strcmp(ResName[ib],"  R") || !strcmp(ResName[ib],"ARG"))
               bseq[i] = 'R';

/* Polar amino acids */
            else if (!strcmp(ResName[ib],"  S") || !strcmp(ResName[ib],"SER"))
               bseq[i] = 'S';
            else if (!strcmp(ResName[ib],"  T") || !strcmp(ResName[ib],"THR"))
               bseq[i] = 'T';
            else if (!strcmp(ResName[ib],"  Y") || !strcmp(ResName[ib],"TYR"))
               bseq[i] = 'Y';
            else if (!strcmp(ResName[ib],"  H") || !strcmp(ResName[ib],"HIS"))
               bseq[i] = 'H';

            else if (!strcmp(ResName[ib],"  C") || !strcmp(ResName[ib],"CYS"))
               bseq[i] = 'C';
            else if (!strcmp(ResName[ib],"  N") || !strcmp(ResName[ib],"ASN"))
               bseq[i] = 'N';
            else if (!strcmp(ResName[ib],"  Q") || !strcmp(ResName[ib],"GLN"))
               bseq[i] = 'Q';
            else if (!strcmp(ResName[ib],"  W") || !strcmp(ResName[ib],"TRP"))
               bseq[i] = 'W';

/* Glycine */
            else if (!strcmp(ResName[ib],"  G") || !strcmp(ResName[ib],"GLY"))
                bseq[i] = 'G';
            else 
                bseq[i] = 'X';
        }else if (RY[i] == -2){ /* others */
            bseq[i] = '?';
        }

    }
}


long residue_ident(char **AtomName, double **xyz, long ib, long ie)
/* identifying a residue as follows:
 *  R-base  Y-base  amino-acid, others [default]
 *   +1        0        -1        -2 [default]
 */
{
    double d1, d2, d3, dcrt = 2.0, dcrt2 = 3.0, temp[4];
    long i=0, id = -2;
    long CA=0, C=0, N=0, N1=0, C2=0, C6=0, N9=0, P=0, C4p=0,C5p=0, O5p=0;
    

    N9 = find_1st_atom(" N9 ", AtomName, ib, ie, "");
    N1 = find_1st_atom(" N1 ", AtomName, ib, ie, "");
    C2 = find_1st_atom(" C2 ", AtomName, ib, ie, "");
    C6 = find_1st_atom(" C6 ", AtomName, ib, ie, "");
    
    if (N1 && C2 && C6) { /* base */
        for (i = 1; i <= 3; i++)
            temp[i] = xyz[N1][i] - xyz[C2][i];
        d1 = veclen(temp);
        for (i = 1; i <= 3; i++)
            temp[i] = xyz[N1][i] - xyz[C6][i];
        d2 = veclen(temp);
        for (i = 1; i <= 3; i++)
            temp[i] = xyz[C2][i] - xyz[C6][i];
        d3 = veclen(temp);
        if (d1 <= dcrt && d2 <= dcrt && d3 <= dcrt2) {
            id = 0;
            if (N9) {
                for (i = 1; i <= 3; i++)
                    temp[i] = xyz[N1][i] - xyz[N9][i];
                d3 = veclen(temp);
                if (d3 >= 3.5 && d3 <= 4.5)        /* ~4.0 */
                    id = 1;
            }
            return id;
        }
    }else { 
        C4p = find_1st_atom(" C4'", AtomName, ib, ie, "");
        C5p = find_1st_atom(" C5'", AtomName, ib, ie, "");
        O5p = find_1st_atom(" O5'", AtomName, ib, ie, "");
        if(C4p && C5p && O5p){ /* backbone only */
            for (i = 1; i <= 3; i++)
                temp[i] = xyz[C5p][i] - xyz[C4p][i];
            d1 = veclen(temp);
            for (i = 1; i <= 3; i++)
                temp[i] = xyz[C5p][i] - xyz[O5p][i];
            d2 = veclen(temp);
            for (i = 1; i <= 3; i++)
                temp[i] = xyz[C4p][i] - xyz[O5p][i];
            d3 = veclen(temp);
            if (d1 <= 1.7 && d2 <= 1.7 && d3 <= 3.0) {
                id = 0;
                return id;
            }
        }else{
            P = find_1st_atom(" P  ", AtomName, ib, ie, "");
            if((ie-ib)==0 && P){ /* P only */
                id = 0;
                return id;
            }
        }
    }
    
    
    CA = find_1st_atom(" CA ", AtomName, ib, ie, "");
    C = find_1st_atom(" C  ", AtomName, ib, ie, "");
    N = find_1st_atom(" N  ", AtomName, ib, ie, "");
    if((ie-ib)==0 && CA && !C && !N){ /* only CA atoms (ie-ib)==1 && */
        id = -1;
        return id;
    }
    
        
    if (!C)                        /* if C does not exist, use N */
        C = find_1st_atom(" N  ", AtomName, ib, ie, "");
    if (CA && C) {
        for (i = 1; i <= 3; i++)
            temp[i] = xyz[CA][i] - xyz[C][i];
        if (veclen(temp) <= dcrt){
            id = -1;
        }
        return id;
        
    }
    
    return id;                        /* other cases */
}


double veclen(double *va)
/* length (magnitude) of a 1-by-3 vector */
{
    
    return sqrt(dot(va, va));
}

long find_1st_atom(const char *str,  char **strmat, long nb, long ne,const  char *idmsg)
/* return index of the first match, or 0 for no-match */
{
    long i=0, num=0;

    num = num_strmatch((char *)str, strmat, nb, ne);

    if (!num) {
        if (strcmp(idmsg, ""))
            printf( "missing \"%s\" atom %s\n", str, idmsg);
        return 0;
    }
    if (num > 1 && strcmp(idmsg, "")) {
        printf( "more than one %s atoms %s\n", str, idmsg);
        printf( "   *****the first atom is used*****\n");
    }
    for (i = nb; i <= ne; i++)
        if (!strcmp(str, strmat[i]))
            break;
    return i;
}

long num_strmatch(char *str, char **strmat, long nb, long ne)
/*  return number of matchs of str in strmat */
{
    long i=0, num = 0;

    for (i = nb; i <= ne; i++)
        if (!strcmp(str, strmat[i]))
            num++;

    return num;
}


void seq_alignment(long i,long j, char **chain_res, double *percent,
                   long *match, long **seq1_match_acc, long **seq2_match_acc)
/* align sequence (take care of several broken)*/
{
    long len1=0, len2=0,  large=0 , sum=0, nmtach=0;
    char *str1_left, *str2_left, *str1_right, *str2_right;
    char *str01, *str02, *str03, *str04;
    long seq1_match[3], seq2_match[3], seq1_match_left[3], seq2_match_left[3];
    long seq1_match_right[3], seq2_match_right[3], right_tmp1=0,right_tmp2=0,m;
    long seq01_match[3], seq02_match[3],seq03_match[3], seq04_match[3];
    char *str_left1, *str_left2, *str_right1, *str_right2 ;
    char *str_left3, *str_left4, *str_right3, *str_right4 ;
    
    
    len1=strlen(chain_res[i]);
    len2=strlen(chain_res[j]);

    if(len1<=len2){  /* use the larger as memory allocation */
        large=len2;
    }else{
        large=len1;
    }
    
        
    str01=cvector(0, large);
    str02=cvector(0, large);
    str03=cvector(0, large);
    str04=cvector(0, large);
    
    
    str1_left=cvector(0, large);
    str2_left=cvector(0, large);
    str1_right=cvector(0, large);
    str2_right=cvector(0, large);
    
    str_left1=cvector(0, large);
    str_left2=cvector(0, large);
    str_left3=cvector(0, large);
    str_left4=cvector(0, large);
    str_right1=cvector(0, large);
    str_right2=cvector(0, large);
    str_right3=cvector(0, large);
    str_right4=cvector(0, large);

    strcpy(str01, chain_res[i]);
    strcpy(str02, chain_res[j]);
    
    get_left_right_str(str01, str02, str1_left, str2_left,
                       str1_right, str2_right, seq1_match, seq2_match);
/*
    fprintf(fout,"\nfirst seq align %d %d: %4d %4d: %4d %4d:\n", i,j,
            seq1_match[1],seq1_match[2],seq2_match[1],seq2_match[2]);
    fprintf(fout, "str1_left=%s; str2_left=%s;\n str1_right=%s; str2_right=%s;\n\n",
            str1_left, str2_left, str1_right, str2_right);
*/
    
    sum= seq1_match[2]-seq1_match[1]+1;
    
    nmtach = 1;
    seq1_match_acc[nmtach][1] = seq1_match[1];
    seq1_match_acc[nmtach][2] = seq1_match[2];
    seq2_match_acc[nmtach][1] = seq2_match[1];
    seq2_match_acc[nmtach][2] = seq2_match[2];
    
        /*--------max 7  broken------------*/

    strcpy(str01,"");
    strcpy(str02,"");
    strcpy(str03,"");
    strcpy(str04,"");
    
    strcpy(str01,str1_left);
    strcpy(str02,str2_left);
    strcpy(str03,str1_right);
    strcpy(str04,str2_right);
    
/*    fprintf(fout, "new str left : str01=%s; str02=%s;\n\n", str01,str02);*/

    if(strlen(str01)>0 && strlen(str02)>0){
        get_left_right_str(str01, str02, str1_left, str2_left,str1_right,
                           str2_right, seq01_match, seq02_match);
 
        strcpy(str_left1,"");
        strcpy(str_left2,"");
        strcpy(str_right1,"");
        strcpy(str_right2,"");
        
        strcpy(str_left1,  str1_left);
        strcpy(str_left2,  str2_left);
        strcpy(str_right1,  str1_right);
        strcpy(str_right2,  str2_right);
        
        m=seq01_match[2]- seq01_match[1];
        
        if(m>0){
            nmtach++;
            seq1_match_acc[nmtach][1] = seq01_match[1];
            seq1_match_acc[nmtach][2] = seq01_match[2];
            seq2_match_acc[nmtach][1] = seq02_match[1];
            seq2_match_acc[nmtach][2] = seq02_match[2];
            sum=sum + m + 1;
        }
/*        
        fprintf(fout, "seq align : str1_left=%s; str2_left=%s;\n str1_right=%s; str2_right=%s;\n\n",
                str1_left, str2_left, str1_right, str2_right);
        fprintf(fout,"\nseq align %d %d: %4d %4d: %4d %4d: sum=%4d \n", i,j,
                seq01_match[1],seq01_match[2],seq02_match[1],seq02_match[2], sum);
*/  

        if(strlen(str_left1)>0 && strlen(str_left2)>0){
            get_left_right_str(str_left1, str_left2, str1_left, str2_left,str1_right,
                               str2_right, seq1_match_left, seq2_match_left);

            m=seq1_match_left[2] - seq1_match_left[1];
            if(m>0){
                nmtach++;
                seq1_match_acc[nmtach][1] = seq1_match_left[1];
                seq1_match_acc[nmtach][2] = seq1_match_left[2];
                seq2_match_acc[nmtach][1] = seq2_match_left[1];
                seq2_match_acc[nmtach][2] = seq2_match_left[2];
                sum=sum + m  +1;
            }
                /*
            fprintf(fout,"\nseq align %d %d: %4d %4d: %4d %4d: sum=%4d \n", i,j,
                seq01_match[1],seq01_match[2],seq02_match[1],seq02_match[2], sum);
                */
        }
        
        if(strlen(str_right1)>0 && strlen(str_right2)>0){
            get_left_right_str(str_right1, str_right2, str1_left, str2_left,str1_right,
                               str2_right, seq1_match_right, seq2_match_right);
            m=seq1_match_right[2] - seq1_match_right[1];
            if(m>0){
                nmtach++;
                seq1_match_acc[nmtach][1] = seq1_match_right[1]+seq01_match[2]+1;
                seq1_match_acc[nmtach][2] = seq1_match_right[2]+seq01_match[2]+1;
                seq2_match_acc[nmtach][1] = seq2_match_right[1]+seq02_match[2]+1;
                seq2_match_acc[nmtach][2] = seq2_match_right[2]+seq02_match[2]+1;
                sum=sum + m +1;
            }
/*            
            fprintf(fout,"\nseq align %d %d: %4d %4d: %4d %4d: sum=%4d \n", i,j,
                seq01_match[1],seq01_match[2],seq02_match[1],seq02_match[2], sum);
*/
        }
        
    }
    
/*    fprintf(fout, "new str : str03=%s; str04=%s;\n\n", str03,str04);*/
    
    if(strlen(str03)>0 && strlen(str04)>0){ /*right side*/
        get_left_right_str(str03, str04, str1_left, str2_left,
                       str1_right, str2_right, seq03_match, seq04_match);
        
        strcpy(str_left3,"");
        strcpy(str_left4,"");
        strcpy(str_right3,"");
        strcpy(str_right4,"");
        
        strcpy(str_left3,  str1_left);
        strcpy(str_left4,  str2_left);
        strcpy(str_right3,  str1_right);
        strcpy(str_right4,  str2_right);
        
        m=seq03_match[2] - seq03_match[1] ;
        
        if(m>0){
            nmtach++;
            seq1_match_acc[nmtach][1] = seq03_match[1]+seq1_match[2]+1;
            seq1_match_acc[nmtach][2] = seq03_match[2]+seq1_match[2]+1;
            seq2_match_acc[nmtach][1] = seq04_match[1]+seq2_match[2]+1;
            seq2_match_acc[nmtach][2] = seq04_match[2]+seq2_match[2]+1;
            right_tmp1=seq1_match_acc[nmtach][2];
            right_tmp2=seq2_match_acc[nmtach][2];
            sum=sum +m + 1;
        }
            /*        
        fprintf(fout,"\nseq align %d %d: %4d %4d: %4d %4d: sum=%4d \n", i,j,
                seq03_match[1],seq03_match[2], seq04_match[1],seq04_match[2],sum );
        fprintf(fout, "seq align : str1_left=%s; str2_left=%s;\n str1_right=%s; str2_right=%s;\n\n",
                str1_left, str2_left, str1_right, str2_right);
            */

        if(strlen(str_left3)>0 && strlen(str_left4)>0){
            get_left_right_str(str_left3, str_left4, str1_left, str2_left,str1_right,
                               str2_right, seq1_match_left, seq2_match_left);
            m=seq1_match_left[2] - seq1_match_left[1];
            
            if(m>0){
                nmtach++;
                seq1_match_acc[nmtach][1] = seq1_match_left[1]+seq1_match[2]+1;
                seq1_match_acc[nmtach][2] = seq1_match_left[2]+seq1_match[2]+1;
                seq2_match_acc[nmtach][1] = seq2_match_left[1]+seq2_match[2]+1;
                seq2_match_acc[nmtach][2] = seq2_match_left[2]+seq2_match[2]+1;
                sum=sum + m +1;
            }

   /*     fprintf(fout,"\nseq align left03 %4d %4d  sum=%4d \n", i,j,sum );*/
            
        }
        if(strlen(str_right3)>0 && strlen(str_right4)>0){
            get_left_right_str(str_right3, str_right4, str1_left, str2_left,str1_right,
                               str2_right, seq1_match_right, seq2_match_right);
            
            m= seq1_match_right[2]- seq1_match_right[1];
            if(m>0){
                nmtach++;
                seq1_match_acc[nmtach][1] = seq1_match_right[1]+right_tmp1+1;
                seq1_match_acc[nmtach][2] = seq1_match_right[2]+right_tmp1+1;
                seq2_match_acc[nmtach][1] = seq2_match_right[1]+right_tmp2+1;
                seq2_match_acc[nmtach][2] = seq2_match_right[2]+right_tmp2+1;
                sum=sum + m+1;
            }
/*            
            fprintf(fout,"\nseq align right03 %4d %4d  sum=%4d \n", i,j,sum );
            printf("\nseq align right03 %4d %4d  sum=%4d \n", i,j,sum );
*/
        }
    
    }
    *match=nmtach;
    *percent=100.0*sum/large;
    printf("seq align percent = %ld %ld   %7.2lf nmatch=%ld  \n", i,j,  *percent,  nmtach);

        /*
    for(m=1; m<=nmtach; m++){
        fprintf(fout,"seq align final  %4d %4d %4d %4d %4d\n", m,
                seq1_match_acc[m][1],seq1_match_acc[m][2],seq2_match_acc[m][1],
                seq2_match_acc[m][2]);
    }
        */
    free_cvector(str01,0, large);
    free_cvector(str02,0, large);
    free_cvector(str03,0, large);
    free_cvector(str04,0, large);
    
    
    free_cvector(str1_left,0, large);
    free_cvector(str2_left,0, large);
    free_cvector(str1_right,0, large);
    free_cvector(str2_right,0, large);
    
    free_cvector(str_left1,0, large);
    free_cvector(str_left2,0, large);
    free_cvector(str_left3,0, large);
    free_cvector(str_left4,0, large);
    free_cvector(str_right1,0, large);
    free_cvector(str_right2,0, large);
    free_cvector(str_right3,0, large);
    free_cvector(str_right4,0, large);

    
}

void get_left_right_str(char *str1, char *str2,  char *str1_left,
                        char *str2_left,char *str1_right, char *str2_right,
                        long *seq1_match,long *seq2_match)
/* get the max match, return new strings (left & right) and sequence number*/
{
    long len=0,len1=0, len2=0, k1=0, k2=0, n=0, n1=0, n2=0,  nmatch=0,  *match;
    long **seq1, **seq2, nseg,  max, large;
    char *sub_str;

    len1=strlen(str1); /*reference*/
    len2=strlen(str2);

    if(len1<=len2){  
        large=len2;
    }else{
        large=len1;
    }

    sub_str=cvector(0, large);
    match=lvector(0, large);
    seq1=lmatrix(0, large, 1,2);
    seq2=lmatrix(0, large, 1,2);

/*    fprintf(fout, "\nrout: input: str1=%s;\n str2=%s;\n", str1,  str2);*/

    k2=0;
    nseg=0;
    for(k1=k2; k1<len1; k1++) {
        len=0;
        get_sub_str(k2, k1, str1, sub_str);

        if(strstr(str2, sub_str)){
            len=strlen(sub_str)-1;  /* take off '\0' */
            n1= strstr(str2, sub_str) - str2; 
            n2=n1+len;
            
            seq1[nseg][1]=k2;
            seq1[nseg][2]=k2+len;
            
            seq2[nseg][1]=n1;
            seq2[nseg][2]=n2;
/*            
            fprintf(fout, "if: %5d %5d ; %5d %5d \n",
                    seq1[nseg][1],seq1[nseg][2],seq2[nseg][1],seq2[nseg][2]);
*/
        }else {
            nseg++;
            k2=k1;
        }
    }
    
    nmatch=1;
    max=0;
    for(k1=0; k1<=nseg; k1++) { /* get the largest match */
        len = seq2[k1][2] - seq2[k1][1] + 1;
        if(len>nmatch) {
            nmatch=len;
            max=k1;
        }
    }
    
    seq1_match[1] = seq1[max][1];
    seq1_match[2] = seq1[max][2];
    seq2_match[1] = seq2[max][1];
    seq2_match[2] = seq2[max][2];

        /*
    for(k1=0; k1<=nseg; k1++) {			
        fprintf(fout, "rout: %5d  : %5d %5d: %5d %5d: \n",   k1, 
                 seq1[k1][1],seq1[k1][2],seq2[k1][1],seq2[k1][2]);
    }
    fprintf(fout, "rout: MAX  : %5d %5d: %5d %5d: \n",    
                 seq1[max][1],seq1[max][2],seq2[max][1],seq2[max][2]);
        */

    n=seq1[max][1];
    new_str(str1, 0 ,  n, str1_left);

    n=seq2[max][1];
    new_str(str2, 0 ,  n, str2_left);

    n=seq1[max][2] + 1;
    new_str(str1, n ,  len1, str1_right);
    
    n=seq2[max][2] + 1;
    new_str(str2, n ,  len2, str2_right);


    free_cvector(sub_str, 0, len1);
    free_lvector(match, 0, len1);
    free_lmatrix(seq1, 0, len1, 1,2);
    free_lmatrix(seq2, 0, len2, 1,2);


}
void new_str(char *str, long n1, long n2, char *newstr)
{
    long i=0, j=0;
    for (i=n1; i<n2; i++){
        if(isspace(str[i])) continue;
        newstr[j]=str[i];
        j++;
    }
    newstr[j]='\0';
}

        

void get_sub_str(long k2, long k1, char *str_small, char *sub_str)
{
	
    long i=0, j=0;
    
    for(i=k2; i<=k1; i++){
        sub_str[j]=str_small[i];
        j++;
    }
    sub_str[j]='\0';
}


void fasta_seq_align(long i,long j, char **chain_res, double *score)
/* align sequence using fasta*/
{
    long  k=0;
    char str[200], tmp[20];
    FILE *fp1, *fp2, *fo;
    fp1=fopen("fasta_inp.aa", "w");
    fp2=fopen("fasta_inp.lib", "w");

    fprintf(fp1,">fast_aa\n%s\n", chain_res[i]);
    fprintf(fp2,">fast_lib\n%s\n", chain_res[j]);

    fclose(fp1);
    fclose(fp2);
    
    sprintf(str, "/home/hyang/struct-align/fasta/bin/fasta3 -q -E 0.02 -O fasta_align.out fasta_inp.aa fasta_inp.lib");
    system(str);

    
    if((fo=fopen("fasta_align.out", "r"))==NULL){
        printf("Can not open file: fasta_align.out (routine: fasta_seq_align)\n");
        return;
    }
    
    while(fgets(str, sizeof str, fo)!=NULL){
        if(strstr(str,"Smith-Waterman score:") && strstr(str,"identity") && strstr(str,";")){
            strcpy(str, strchr(str,';'));
            if(sscanf(str+1, "%s",tmp) !=1 ) continue;
            k=strlen(tmp);
            tmp[k-1]='\0';
            *score=atof(tmp);
            
            break;
            
        }else if(!strstr(str,"Smith-Waterman score:") && strstr(str,"identity")){
            if(sscanf(str+1, "%s",tmp) !=1 ) continue;
            k=strlen(tmp);
            tmp[k-1]='\0';
            *score=atof(tmp);            
            break;
        }
        
    }
    
    if(*score>50.0)
        printf("fastaa seq align %4ld  %4ld %7.2lf \n", i,j,   *score );

    fclose(fo);
/*    system("rm -f fasta_align.out fasta_inp.aa fasta_inp.lib");*/
    
    
}



void struct_3d_align(long i,long j,long larg,long **chain_idx,char **AtomName, 
                    long **seidx, double **xyz, long *RY,  double *rmsd)
/* 3D struct-alignment by shifting the shorter chain (i) along the longer one
   (j). If there is a break in the sequnence, this method may not be good.
*/
{
    long k=0, m=0, n1=0, nmatch=0, npair=0, small=0;
    long large=0, m1=0, m2=0, k1=0, nterm=0;
    double orgi[4], *rmsd_match, rmsd_fit=0, svalue;
    double **e_xyz, **s_xyz, **fitted_xyz, **tmpe_xyz,**tmps_xyz, **R;


    e_xyz=dmatrix(0, larg, 0, 4);
    s_xyz=dmatrix(0, larg, 0, 4);
    tmpe_xyz = dmatrix(0, larg, 0, 4);
    tmps_xyz = dmatrix(0, larg, 0, 4);
    fitted_xyz = dmatrix(0, larg, 0, 4);
    rmsd_match = dvector(0, larg );
    R = dmatrix(1, 4, 1, 4);

    get_es_xyz(i,chain_idx, AtomName, seidx, xyz, RY, tmpe_xyz);
    get_es_xyz(j,chain_idx, AtomName, seidx, xyz, RY, tmps_xyz);

    small=chain_idx[i][2] - chain_idx[i][1] + 1;  
    large=chain_idx[j][2] - chain_idx[j][1] + 1; 
	
    nterm=(long)( small*0.05); /* 5% shift in front */

    nmatch=0;
    for(k=nterm; k>=2 ; k--) {
        m1=0;
        for(k1=k; k1<=small ; k1++) {
            m1++;
            for(n1=1; n1<=3; n1++)
                e_xyz[m1][n1]=tmpe_xyz[k1][n1];
        }
        npair=m1;
        m2=0;
        for(k1=1; k1<=npair ; k1++) {
            m2++;
            for(n1=1; n1<=3; n1++)
                s_xyz[m2][n1]=tmps_xyz[k1][n1];
        }
        if(m1!=m2)printf("Problem with match (front) %ld %ld\n",m1,m2);

        ls_fit(s_xyz, e_xyz, npair, &rmsd_fit, fitted_xyz, R, orgi);
        nmatch++;
        rmsd_match[nmatch]=rmsd_fit;
    }
   
    for(k=1; k<=small; k++) { /*all small chain shift from j=1 to j=large */
        for(n1=1; n1<=3; n1++)
            e_xyz[k][n1]=tmpe_xyz[k][n1];
    }
    npair=small;
    for(k=0; k<= large-small ; k++) {
        m=0;
        for(k1=1+k; k1<=npair+k && k1<=large ; k1++) {
            m++;
            for(n1=1; n1<=3; n1++)
                s_xyz[m][n1]=tmps_xyz[k1][n1];
        }
        ls_fit(s_xyz, e_xyz, npair, &rmsd_fit, fitted_xyz, R, orgi);
        nmatch++;
        rmsd_match[nmatch]=rmsd_fit;
    }

    for(k=1; k<=nterm ; k++) { /* shift in the end */
        m1=0;
        for(k1=1; k1<=small-k ; k1++) {
            m1++;
            for(n1=1; n1<=3; n1++)
                e_xyz[m1][n1]=tmpe_xyz[k1][n1];
        }
        npair=m1;
        m2=0;
        for(k1=large-npair+1; k1<=large ; k1++) {
            m2++;
            for(n1=1; n1<=3; n1++)
                s_xyz[m2][n1]=tmps_xyz[k1][n1];
        }			
        if(m1!=m2)printf("Problem with match (end) %ld %ld\n",m1,m2);
        
        ls_fit(s_xyz, e_xyz, npair, &rmsd_fit, fitted_xyz, R, orgi);
        nmatch++;
        rmsd_match[nmatch]=rmsd_fit;
    }

    svalue=1000.0;
    for(k=1; k<=nmatch ; k++) {
        if(rmsd_match[k]<=svalue)
            svalue=rmsd_match[k];
/*            
        printf("k, i,j,small,large,rmsd %2d  %2d %2d %4d %4d %7.2f %7.2f\n", 
                k, i, j, small,large, rmsd_match[k], svalue);
*/            
    }
    *rmsd=svalue;
    
    free_dmatrix(e_xyz, 0, larg, 0, 4);
    free_dmatrix(s_xyz, 0, larg, 0, 4);
    free_dmatrix(tmpe_xyz , 0, larg, 0, 4);
    free_dmatrix(tmps_xyz , 0, larg, 0, 4);
    free_dmatrix(fitted_xyz,0, larg, 0, 4);
    free_dvector(rmsd_match,0, larg );
    free_dmatrix(R , 1, 4, 1, 4);
}



void struct_3d_align_new(long i1,long j1,long larg,long **chain_idx,char **AtomName,
                         long **seidx,double **xyz,long *RY, long nmatch,
                         long **seq1_match, long **seq2_match, double *rmsd)
/* 3D struct-alignment after sequence alignment. */
{
    long i=0,j=0,k=0, n=0, npair=0,natom=0;
    double orgi[4], *rmsd_match, rmsd_fit=0;
    double **e_xyz, **s_xyz, **fitted_xyz, **tmpe_xyz,**tmps_xyz, **R;

    e_xyz=dmatrix(0, larg, 0, 4);
    s_xyz=dmatrix(0, larg, 0, 4);
    tmpe_xyz = dmatrix(0, larg, 0, 4);
    tmps_xyz = dmatrix(0, larg, 0, 4);
    fitted_xyz = dmatrix(0, larg, 0, 4);
    rmsd_match = dvector(0, larg );
    R = dmatrix(1, 4, 1, 4);
    
    get_es_xyz(i1,chain_idx, AtomName, seidx, xyz, RY, tmpe_xyz);
    get_es_xyz(j1,chain_idx, AtomName, seidx, xyz, RY, tmps_xyz);

    natom=0;
    npair=0;

    for(i=1; i<=nmatch ; i++) {
        for(j=seq1_match[i][1]; j<=seq1_match[i][2]; j++) {
            k=j+1;
			
            natom++;
            for(n=1; n<=3; n++)
                e_xyz[natom][n]=tmpe_xyz[k][n];
        }
        for(j=seq2_match[i][1]; j<=seq2_match[i][2]; j++) {
            k=j+1;
            npair++;
            for(n=1; n<=3; n++)
                s_xyz[npair][n]=tmps_xyz[k][n];
        }
    }
    if(npair != natom )
        printf("Problem with LS fitting! npair=%ld natom=%ld\n",natom, npair);
         
    ls_fit(s_xyz, e_xyz, npair, &rmsd_fit, fitted_xyz, R, orgi);
   
    *rmsd=rmsd_fit;
        /*
     printf("LS fitting! npair=%d natom=%d  nmatch=%d rmas=%7.2f\n",
     natom, npair, nmatch,*rmsd);
        */
    free_dmatrix(e_xyz, 0, larg, 0, 4);
    free_dmatrix(s_xyz, 0, larg, 0, 4);
    free_dmatrix(tmpe_xyz , 0, larg, 0, 4);
    free_dmatrix(tmps_xyz , 0, larg, 0, 4);
    free_dmatrix(fitted_xyz,0, larg, 0, 4);
    free_dvector(rmsd_match,0, larg );
    free_dmatrix(R , 1, 4, 1, 4);
}


void get_es_xyz(long i, long **chain_idx, char **AtomName, long **seidx, 
                double **xyz, long *RY, double **e_xyz)
{
    long natom=0, n=0, k=0, m=0;
    natom=0;

    for (k=chain_idx[i][1]; k<=chain_idx[i][2]; k++){
        
        if(RY[k]==-1){
            for(m=seidx[k][1]; m<=seidx[k][2]; m++){
                if (!strcmp(AtomName[m], " CA ")){
                    natom++;
                    
                    for(n=1; n<=3; n++)
                        e_xyz[natom][n]=xyz[m][n];
                }
            }
        }else if (RY[k] >=0 ){
            for(m=seidx[k][1]; m<=seidx[k][2]; m++){
                if (!strcmp(AtomName[m], "  P ")){
                    natom++;
                    for(n=1; n<=3; n++)
                        e_xyz[natom][n]=xyz[m][n];
                }
            }
        }
    }
}


void broken_chain(long id, long **chain_idx, long **seidx, double **xyz, 
                  long *RY, char **AtomName, long **broken_chain_idx,long *nbr)
/* find if chain is broken. broken_chain_idx[i]=residue_number
For Nucleic acid: use xyz of C3'-P (between i,i+1) or O3'-O5' (between i,i+1).
For protein:  use xyz of Ca-Ca (between i,i+1).
 */
{
    long p=0, p1=0,p2=0, c3p=0, o3p=0, o5p=0, c=0, ca1=0, ca2=0 ;
    long i=0, k=0, ib1=0,ie1=0, ib2=0,ie2=0, nbroken=0;
    double dxyz[4];
    
    nbroken=0;
    for (i = chain_idx[id][1]; i <= chain_idx[id][2]-1; i++) {
        ib1 = seidx[i][1];
        ie1 = seidx[i][2];
        ib2 = seidx[i+1][1];
        ie2 = seidx[i+1][2];
        if (RY[i] >= 0){
            c3p = find_1st_atom(" C3'", AtomName, ib1, ie1, "");
            p   = find_1st_atom(" P  ", AtomName, ib2, ie2, "");
            o3p = find_1st_atom(" O3'", AtomName, ib1, ie1, "");
            o5p = find_1st_atom(" O5'", AtomName, ib2, ie2, "");
            if(c3p && p){
                for(k=1; k<=3; k++){
                    dxyz[k]=xyz[c3p][k]-xyz[p][k];
                }
                if(veclen(dxyz) > 3.2){
                    nbroken++;
                    broken_chain_idx[id][nbroken]=i;
                }
                c3p = 0;
                p = 0;                
            }else if(o3p && o5p){ 
                o3p = find_1st_atom(" O3'", AtomName, ib1, ie1, "");
                o5p = find_1st_atom(" O5'", AtomName, ib2, ie2, "");
                if(o3p && o5p){
                    for(k=1; k<=3; k++){
                        dxyz[k]=xyz[o5p][k]-xyz[o3p][k];
                    }
                    if(veclen(dxyz) > 3.5){
                        nbroken++;
                        broken_chain_idx[id][nbroken]=i;
                    }
                    o3p = 0;
                    o5p = 0;                
                }
            }else {
                p1 = find_1st_atom(" P  ", AtomName, ib1, ie1, "");
                p2 = find_1st_atom(" P  ", AtomName, ib2, ie2, "");
                if(p1 && p2){
                    for(k=1; k<=3; k++){
                        dxyz[k]=xyz[p2][k]-xyz[p1][k];
                    }
                    if(veclen(dxyz) > 7.5){
                        nbroken++;
                        broken_chain_idx[id][nbroken]=i;
                    }
                    p1 = 0;
                    p2 = 0;                
                }
            }

            
        }else if(RY[i] <= -1){ /*protein*/
              c = find_1st_atom(" C  ", AtomName, ib1, ie1, "");
            ca1 = find_1st_atom(" CA ", AtomName, ib1, ie1, "");
            ca2 = find_1st_atom(" CA ", AtomName, ib2, ie2, "");
            if(c && ca2){
                for(k=1; k<=3; k++){
                    dxyz[k]=xyz[ca2][k]-xyz[c][k];
                }
                if(veclen(dxyz) > 3.0){
                    nbroken++;
                    broken_chain_idx[id][nbroken]=i;
                }
            }else if(ca1 && ca2){
                for(k=1; k<=3; k++){
                    dxyz[k]=xyz[ca2][k]-xyz[ca1][k];
                }
                if(veclen(dxyz) > 4.2){
                    nbroken++;
                    broken_chain_idx[id][nbroken]=i;
                }
            }
            
        }
        
    }
    nbr[id]=nbroken;
    
}


void delete_ligand(long num, char **AtomName, char **ResName, char *ChainID,
                   long *ResSeq, double **xyz, char **Miscs, long *natom)
/* get  rid of ligands  */
{
    long i=0, j=0, k=0, m=0, atoms=0, ib=0, ie=0;
    long num_residue=0, nchain=0, nprot=0, nacid=0, alt=0;
    long  *ry, **chain_idx, **seidx;
    char str_id0[20], str_id1[20];

    seidx = residue_idx(num, ResSeq, Miscs, ChainID, ResName, &num_residue);

    ry = lvector(0, num_residue);
    chain_idx = lmatrix(1,500 , 1, 2);  /* # of chains max = 500 */    

    get_chain_idx(num_residue, seidx, ChainID,  ry, &nchain, chain_idx);
    
    atoms = 0;
    for (i=1; i<=nchain; i++){
        nprot=0;
        nacid=0;
            /*
        ib=chain_idx[i][1];
        ie=chain_idx[i][2];
        printf(" %4d %4d chain_ID:  %c  from residue %4d to %4d\n",
               chain_idx[i][1],chain_idx[i][2],ChainID[ seidx[ib][1]],
               ResSeq[ seidx[ib][1] ], ResSeq[ seidx[ie][1] ]);
            */
        for (k=chain_idx[i][1]; k<=chain_idx[i][2]; k++){
            ib = seidx[k][1];
            ie = seidx[k][2];
            ry[k] = residue_ident(AtomName, xyz, ib, ie);
            
            if(ry[k] == -2) continue;
            if(ry[k] == -1) nprot++;
            if(ry[k] >=0 ) nacid++;
        }

        for (k=chain_idx[i][1]; k<=chain_idx[i][2]; k++){
            ib = seidx[k][1];
            ie = seidx[k][2];
            if(ry[k] == -2) continue;
            if(nprot>nacid){ /* it must be protein */
                if(ry[k] >=0)continue;
            }else  /* it must be nucleic acids*/
                if(ry[k] == -1)continue;
            for(j = ib; j <= ie; j++){
                sprintf(str_id0, "%s%s%c%ld",AtomName[j], ResName[j], 
                        ChainID[j], ResSeq[j]); 
                alt=0;
                for(m = j+1; m <=ie; m++){
                    sprintf(str_id1, "%s%s%c%ld",AtomName[m], ResName[m], 
                            ChainID[m], ResSeq[m]); 
                    if(!strcmp(str_id1, str_id0)){
                        alt++;
                        break;
                    }
                }
                if(alt>0) continue;

                atoms++;
                strcpy(AtomName[atoms], AtomName[j]);
                strcpy(ResName[atoms], ResName[j]);
                ChainID[atoms] = ChainID[j];
                ResSeq[atoms] = ResSeq[j];
                for(m = 0 ; m <=NMISC; m++)
                    Miscs[atoms][m] = Miscs[j][m];
                for(m = 1 ; m <=3; m++)
                    xyz[atoms][m] = xyz[j][m];
                
            }
            
        }
        
    }
    *natom = atoms;
    free_lvector(ry , 0, num_residue);
    free_lmatrix(chain_idx , 1,500 , 1, 2);  
    free_lmatrix(seidx, 0, num_residue, 0, 2);
}