src/common/target_struct.cc

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <assert.h>
#include <string>



#include <libxml/tree.h>
#include <libxml/parser.h>
#include <libxml/xpath.h>
#include <libxml/xpathInternals.h>

#include "open_prospect.h"
#include "sm_matrix.h"

//declare shared non-public functions
char *get_first_child_content( xmlNodePtr parent, char *str);


//a non-public function found in seq_io.cc
void append_to_str(char **str1, char *str2);

//a non-public function found in rotamer.cc
int rot_mask_active(RotamerMask *mask, int i);

/* read sequence file in FASTA or other general formats            */

TargetStruct::TargetStruct(void) {
        this->len = 0;
        this->residue_offset = 0;
        this->twobody = NULL;
        this->hbond = NULL;
        this->bsheet = NULL;
        this->distconst = NULL;
        this->distconst_count = 0;
        this->residue = NULL;  
        this->structure = NULL;
}

TargetStruct::~TargetStruct() {
        FreeMem();
}

void TargetStruct::FreeMem(void) {      
        if (this->residue != NULL)
                delete this->residue;
        this->residue = NULL;    
        if (this->structure != NULL)
                delete this->structure;
        this->structure = NULL;    
        
        if(this->twobody != NULL) {
                free(this->twobody[0]);
                free(this->twobody);
        }
        twobody = NULL;
        if (this->distconst != NULL) {
                free(this->distconst);
        }
        distconst = NULL;
        if (this->hbond != NULL) {
                free( this->hbond[0] );
                free( this->hbond );
        }
        hbond = NULL;
        if (this->bsheet != NULL) {
                free( this->bsheet[0] );
                free( this->bsheet );
        }
        bsheet = NULL;
        len = 0;
}


TargetStruct::TargetStruct(const TargetStruct &rhs) {
        len = 0;
        twobody = NULL;
        hbond = NULL;
        distconst = NULL;
        bsheet = NULL;
        distconst_count = 0;
        residue = NULL;  
        residue_offset = rhs.residue_offset;
        
        if ( rhs.len > 0 && rhs.residue != NULL) {
                len = rhs.len;
                residue = new TargetResidue[len];       
                for (int i = 0; i < len; i++) {
                        residue[i] = rhs.residue[i];            
                }
                if ( rhs.twobody ) {
                        twobody = (pValType **)malloc(sizeof(pValType *) * len);
                        twobody[0] = (pValType *)malloc(sizeof(pValType *) * len * len );
                        for (int i = 0; i < len; i++) {
                                twobody[i] = twobody[0] + i * len;
                                for (int j = 0; j < len; j++) 
                                        twobody[i][j] = rhs.twobody[i][j];
                        }                       
                }
                if ( rhs.hbond ) {
                        hbond = (pValType **)malloc(sizeof(pValType *) * len);
                        hbond[0] = (pValType *)malloc(sizeof(pValType *) * len * len );
                        for (int i = 0; i < len; i++) {
                                hbond[i] = hbond[0] + i * len;
                                for (int j = 0; j < len; j++) 
                                        hbond[i][j] = rhs.hbond[i][j];
                        }                       
                }       
                if ( rhs.bsheet ) {
                        bsheet = (char **)malloc(sizeof(char *) * len);
                        bsheet[0] = (char *)malloc(sizeof(char *) * len * len );
                        for (int i = 0; i < len; i++) {
                                bsheet[i] = bsheet[0] + i * len;
                                for (int j = 0; j < len; j++) 
                                        bsheet[i][j] = rhs.bsheet[i][j];
                        }                       
                }                       
        }       
        if ( rhs.distconst && rhs.distconst_count > 0) {
                distconst_count = rhs.distconst_count;
                distconst = (TargetDistConstraint *)malloc(sizeof(TargetDistConstraint) * distconst_count);
                for (int i = 0; i < distconst_count; i++) {
                        distconst[i] = rhs.distconst[i];
                }
        }       
}

void TargetStruct::SetLen(int new_len) {
        if ( new_len != len ) {
                FreeMem();
                len = new_len;
                residue = new TargetResidue[len];
                structure = new ResidueAtoms[ len ];
        }
};

int TargetStruct::Setup( ProspectParam *param ) {
        int i, j;
        pValType sum;
        
        if (this->residue == NULL)
                return 1;
        
        for (i = 0; i < this->len; i++) {
                sum = 0;
                for (j = 0; j < 20; j++)
                        sum += this->residue[i].FREQ[j];
                if (sum < 0.000001) {
                        this->residue[i].FREQ[ AA1toAANum( this->residue[i].RS ) ] = 1;
                }
        }       
        this->twobody = (pValType **)malloc( sizeof(pValType *) * this->len );
        this->twobody[0] = (pValType *)malloc( sizeof(pValType) * this->len * this->len );
        for (i = 1; i < this->len; i++) {
                this->twobody[i] = this->twobody[0] + i*this->len;
        }  
        for (i = 0; i < this->len; i++) {
                for (j = 0; j <= i; j++) {
                        this->twobody[i][j] = calc_twobody_val(param, this,i,j);
                        this->twobody[j][i] = this->twobody[i][j];
                }
        }       
        if ( structure ) {
                for (i = 0; i < this->len; i++) {
                        for (int j = 0; j < this->structure[i].atom_count; j++) {
                                /*
                                 for (int k = 0; k < param->atom_data_count; k++) {
                                         if ( this->residue[i].structure.atom_name[j] == PDBAtomtoNum( param->atom_data[k].name ) ) {
                                                 this->residue[i].structure.radius[j] = param->atom_data[k].radius;
                                         }
                                 }
                                 */
                                this->structure[i].radius[j] =
                                param->atom_data->FindRadius( this->structure[i].atom_name[j] );
                        }
                }
        }
        return 0;
}

int TargetStruct::SetSequence(char *seq) {
        SetLen(  strlen( seq ) );
        for (int i = 0; i < this->len; i++) {
                this->residue[i].RS = toupper(seq[i]);
                for (int j = 0; j < 21; j++)
                        this->residue[i].FREQ[j] = 0;
                this->residue[i].FREQ[ AA1toAANum(this->residue[i].RS) ] = 1;
        }
        this->twobody = NULL;  
        this->distconst_count = 0;
        return 0;
}


void TargetStruct::SetChain(char new_chain) {
        for ( int i = 0; i < len; i++) {
                residue[i].chain = new_chain;
                if ( structure ) {
                        structure[i].chain = new_chain;
                }
        }       
}


/****************************************************************
*****File I/O Stuff**********************************************
****************************************************************/

void TargetStruct::ReadSingleDelStr( char *data ) {
        char *buffer = strdup(data) , RS, *cur_ptr;
        float tmp_float;
        int i = 0;
        cur_ptr = strtok( data, "\n\r" );
        do {
                sscanf( cur_ptr, "%c %f", &RS, &tmp_float );
                if ( i < len )
                        residue[i].single_del = tmp_float;
                i++;
        } while ( (cur_ptr = strtok( NULL, "\n\r") ) );
        free(buffer);
}

void TargetStruct::ReadSingleInsStr( char *data ) {
        char *buffer = strdup(data) , RS, *cur_ptr;
        float tmp_float;
        int i = 0;
        cur_ptr = strtok( data, "\n\r" );
        do {
                sscanf( cur_ptr, "%c %f", &RS, &tmp_float );
                if ( i < len )
                        residue[i].single_ins = tmp_float;
                i++;
        } while ( (cur_ptr = strtok( NULL, "\n\r") ) );
        free(buffer);
}


int TargetStruct::LoadXML(char *xml_path) {     
        char *read_buffer;
        xmlDocPtr doc;
        xmlNodePtr root;
        
        //  fprintf(stderr, "Loading %s\n", template_path);
        doc = xmlParseFile(xml_path);
        if (doc == NULL)
                return 1;
        
        root = xmlDocGetRootElement(doc);
        if (root == NULL)
                return 1;
        //Load in the residue information from the file         
        read_buffer = get_first_child_content(root, "ssp");
        if ( read_buffer ) {
                LoadSSPStr( read_buffer );
                free(read_buffer);
        }
        read_buffer = get_first_child_content(root, "singledel");
        if ( read_buffer ) {
                ReadSingleDelStr( read_buffer );
                free(read_buffer);
        }               
        read_buffer = get_first_child_content(root, "singleins");
        if ( read_buffer ) {
                ReadSingleInsStr( read_buffer );
                free(read_buffer);
        }       
        return 0;       
}


int TargetStruct::LoadSSP(char *ssp_path) {
        char buffer[2000];
        FILE *ssp_file;
        ssp_file = fopen(ssp_path, "r");
        if (ssp_file == NULL) {
                return 1;
        }       
        std::string my_string;
        while (fgets(buffer, sizeof(buffer), ssp_file) ) {
                my_string += buffer;            
        }
        fclose(ssp_file);       
        return LoadSSPStr( (char *)my_string.c_str() ); 
}


int TargetStruct::LoadSSPStr(char *ssp_text) {
        int i,j,num;
        char dummyc;
        char *buffer= strdup( ssp_text );
        TargetResidue  *tmp_residue;
        float sum;
        float ss_prob[20];

        this->distconst = NULL;
        tmp_residue = new TargetResidue[5000];
        
        i=0;
        
        char *cur_ptr = strtok( buffer, "\n\r" );
        while ( cur_ptr ) {
                j = sscanf(cur_ptr, "%d %c %c %f %f %f %c %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f",
                                   &num,
                                   &tmp_residue[i].RS,
                                   &tmp_residue[i].SS,
                                   &tmp_residue[i].SS_PROB[LOOP],
                                   &tmp_residue[i].SS_PROB[SHEET],
                                   &tmp_residue[i].SS_PROB[HELIX],
                                   &dummyc,
                                   &ss_prob[0],
                                   &ss_prob[1],
                                   &ss_prob[2],
                                   &ss_prob[3],
                                   &ss_prob[4],
                                   &ss_prob[5],
                                   &ss_prob[6],
                                   &ss_prob[7],
                                   &ss_prob[8],
                                   &ss_prob[9],
                                   &ss_prob[10],
                                   &ss_prob[11],
                                   &ss_prob[12],
                                   &ss_prob[13],
                                   &ss_prob[14],
                                   &ss_prob[15],
                                   &ss_prob[16],
                                   &ss_prob[17],
                                   &ss_prob[18],
                                   &ss_prob[19] );
                if (j != 27) {
                        return 1;
                }               
                sum = 0.0;
                for (j = 0; j < 20; j++) {
                        tmp_residue[i].FREQ[j] = ss_prob[j];
                        sum += ss_prob[j];
                }
                
                tmp_residue[i].FREQ[20] = 0;
                if (sum > 0.0) {
                        for (j=0;j<20;j++) {
                                tmp_residue[i].FREQ[j] = tmp_residue[i].FREQ[j]/sum;
                                assert(!isnan(tmp_residue[i].FREQ[j]));
                        }
                } else {
                        //bug: need to set some sort of default value
                } 
                tmp_residue[i].orig_num = i;
                i++;            
                cur_ptr = strtok( NULL, "\n\r" );
        }
        if ( i == 0 ) {
                delete tmp_residue;
                return 1;
        }
        SetLen( i );
        for (i = 0; i < len; i++) {
                this->residue[i] = tmp_residue[i];
        }
        delete tmp_residue;
        this->distconst_count = 0;
        return 0;
}


int TargetStruct::LoadConstraint(char *constraint_file) {
        FILE *file;
        char buffer[2048];
        int const_count = 0, const_size = 10, tmp;
        TargetDistConstraint *constraint = (TargetDistConstraint *)malloc(sizeof(TargetDistConstraint) * const_size);
        
        file = fopen(constraint_file, "r");
        
        while (fgets(buffer, sizeof(buffer), file)) {
                if (const_count >= const_size) {
                        const_size += 10;
                        constraint = (TargetDistConstraint *)realloc(constraint, sizeof(TargetDistConstraint) * const_size);
                }
                sscanf(buffer, "%d %d %f %f", &constraint[const_count].aa[0],
                           &constraint[const_count].aa[1],
                           &constraint[const_count].dist, &constraint[const_count].var);
                //make sure the second residue is after the first
                if ( constraint[const_count].aa[0] > constraint[const_count].aa[1] ) {
                        tmp = constraint[const_count].aa[0];
                        constraint[const_count].aa[0] = constraint[const_count].aa[1];
                        constraint[const_count].aa[1] = tmp;
                }
                const_count++;
        }
        fclose(file);
        
        this->distconst = constraint;
        this->distconst_count = const_count;
        return 0;
}





int TargetStruct::LoadPsiPred(char *psipred_path) {
        int i,j,num;
        char dummyc;
        FILE *psipred_file;
        i=0;
        TargetResidue  *tmp_residue;
        
        FreeMem();

        psipred_file = fopen(psipred_path, "r");
        if (psipred_file == NULL) {
                return 1;
        }
        
        this->distconst = NULL;
        tmp_residue = new TargetResidue[5000];
        while (!feof(psipred_file)) {
                fscanf(psipred_file,"%d%c%c%c%c%f%f%f",&num,&dummyc,
                           &tmp_residue[i].RS,
                           &dummyc,
                           &tmp_residue[i].SS,
                           &tmp_residue[i].SS_PROB[LOOP],
                           &tmp_residue[i].SS_PROB[HELIX],
                           &tmp_residue[i].SS_PROB[SHEET]);             
                for (j=0; j < 21;j++) {
                        tmp_residue[i].FREQ[j] = 0;
                }
                tmp_residue[i].FREQ[ AA1toAANum( tmp_residue[i].RS ) ] = 1;
                tmp_residue[i].orig_num = i;
                i++;
        }
        fclose(psipred_file);
        SetLen( num );
        for (i = 0; i < num; i++) {
                this->residue[i] = tmp_residue[i];
        }
        delete tmp_residue;
        this->distconst_count = 0;
        return 0;
}




int TargetStruct::LoadFasta(char *SeqPath) { 
        FILE *fd;
        short i, j = 0;
        char s[1024];
        int QuerySeqSize = 0;
        
        FreeMem();

        if (!(fd = fopen (SeqPath, "r"))) { 
                fprintf (stderr, "INPUT %s OPEN ERROR for %s:%d QUIT\n", SeqPath, __FILE__, __LINE__ );
                return 1; 
        }
        
        while (fgets (s, sizeof(s), fd)) {
                if (strncmp (s, "REM ", 4) != 0 && strncmp (s, ">", 1) != 0)  { 
                        for (i = 0; s[i] != '\0'; i++)  { 
                                if (s[i] != '\n' && s[i] != '\t' && s[i] != ' '  && s[i] != 13 && s[i] != 26) { 
                                        QuerySeqSize++;
                                }
                        }
                }
        }
        rewind(fd);
        
        j = 0;
        SetLen( QuerySeqSize );
        while (fgets (s, sizeof(s), fd)) { 
                if (strncmp (s, "REM ", 4) != 0 && strncmp (s, ">", 1) != 0)  { 
                        for (i = 0; s[i] != '\0'; i++) { 
                                if (s[i] != '\n' && s[i] != '\t' && s[i] != ' '  && s[i] != 13 && s[i] != 26) { 
                                        /* printf ("j=%d: %c -- %d\n", j, s[i], s[i]); */
                                        this->residue[j].RS = toupper(s[i]); j++; 
                                        for (int k = 0; k < 21; k++)
                                                this->residue[j].FREQ[k] = 0;
                                        this->residue[j].FREQ[ AA1toAANum(this->residue[j].RS) ] = 1;
                                }
                        }
                }
        }
        fclose (fd);
        this->twobody = NULL;  
        this->distconst_count = 0;
        return 0;
}




/* read secondary structure information about the seq in PHD format  */

int TargetStruct::LoadPHD( char *filename) { 
        FILE *fd;
        short i, k = -1;
        char     s[1024], *ptr;
        short QuerySeqSize = 0;
        
        
        
        
        if (!(fd = fopen (filename, "r"))) { 
                fprintf (stderr, "INPUT %s OPEN ERROR: %s:%d\n", filename, __FILE__, __LINE__); 
                return 1;
        }
        
        QuerySeqSize = 0;
        while (fgets (s, sizeof(s), fd)) { 
                if ( strstr(s, "AA") ) {
                        ptr = strstr(s, "|");
                        if (ptr ==NULL) {
                                fprintf(stderr, "Invalid PHD file\n");
                                return -1;
                        }
                        ptr++;
                        for (i = 0; i < 60; i++) {
                                if ( isalpha( ptr[i] ) ) 
                                        QuerySeqSize++;
                        }
                }
        }
        
        rewind(fd);
        
        //allocate the memory
        
        SetLen( QuerySeqSize );
        
        //read in the sequence and the structure  
        while (fgets (s, sizeof(s), fd)) { 
                if ( strstr(s, "AA") ) {
                        k++;
                        ptr = strstr(s, "|");
                        if (ptr == NULL) {
                                return -1;
                        }
                        ptr++;
                        for (i = 0; i < 60; i++) {
                                if ( isalpha(ptr[i]) ) { 
                                        this->residue[60*k+i].RS = toupper( ptr[i] );
                                }
                        }
                } else if (strstr (s, "prH") ) { 
                        ptr = strstr(s, "|");
                        ptr++;
                        for (i = 0; i < 60; i++) { 
                                if (isdigit(ptr[i]) && 60*k+i < QuerySeqSize) 
                                        this->residue[60*k+i].SS_PROB[HELIX] = ((pValType)(ptr[i] - '0')) * .1; 
                        }
                } else if (strstr (s, "prE") ) { 
                        ptr = strstr(s, "|");
                        ptr++;
                        for (i = 0; i < 60; i++) { 
                                if (isdigit(ptr[i]) && 60*k+i < QuerySeqSize)
                                        this->residue[60*k+i].SS_PROB[SHEET] = ((pValType) (ptr[i] - '0')) * .1; 
                        }
                } else if (strstr (s, "prL" ) ) { 
                        ptr = strstr(s, "|");
                        ptr++;
                        for (i = 0; i < 60; i++) { 
                                if (isdigit(ptr[i]) && 60*k+i < QuerySeqSize)
                                        this->residue[60*k+i].SS_PROB[LOOP]= ((pValType) (ptr[i] - '0')) * .1; 
                        }     
                } 
                /*
                 else if (strstr(s, "Rel" ) ) { 
                         ptr = strstr(s, "|");
                         ptr++;
                         for (i = 0; i < 60; i++) { 
                                 if (isdigit(ptr[i]) && 60*k+i < QuerySeqSize)
                                         StructureInfo[StructTypeCount][60*k+i] = (short) (ptr[i] - '0'); 
                         }
                 }
                 */
        }
        fclose (fd); 
        
        this->distconst_count = 0;
        
        return 0;
}



int TargetStruct::LoadChk(char *ChkPath) { 
        long i, j, len;
        char *qseq;
        double tmp_dbl[20];
        FILE *file = fopen(ChkPath, "rb");
        
        fread( &len, sizeof(long), 1, file);

        SetLen(len);
        
        qseq = (char *)malloc(sizeof(char) * len);
        fread(qseq, 1, len, file);
        
        for (i = 0; i < len; i++) 
                this->residue[i].RS = qseq[i];
        
        for (i = 0; i < len; i++) {
                
                fread(tmp_dbl, sizeof(double), 20, file);
                for (j = 0; j < 20; j++)
                        this->residue[i].FREQ[j] = tmp_dbl[j];
                this->residue[i].FREQ[20] = 0;
        }
        fclose(file);
        
        this->distconst_count = 0;
        
        return 0;
}


/*
 int ReadChkMem(char *chk_data, char **Seq, double ***freq, short *QuerySeqSize) { 
         long i, j, len;
         char *qseq;
         char *seq_start;
         double *freq_start;
         
         len = ((long *)chk_data)[0];
         *QuerySeqSize = len;
         *Seq = (char *)malloc( len + 1);
         seq_start = chk_data + sizeof(long);
         freq_start = (double *)(chk_data + sizeof(long) + len);
         
         for (i = 0; i < len; i++) {
                 (*Seq)[i] = seq_start[i];
         }
         (*Seq)[len] = 0;
         
         *freq = (double **)malloc(sizeof(double *) * (len) );
         (*freq)[0] = (double *)malloc(sizeof(double) * (len) * 20);
         for (i = 0; i < len; i++) {
                 (*freq)[i] = (*freq)[0] + (i * 20);
                 for (j = 0; j < 20; j++) 
                         (*freq)[i][j] = freq_start[(i * 20) + j];
         }
         return 0;    
 }
 */

int TargetStruct::LoadFreq(char *FreqPath) {
        long i,j;
        int temp_int;
        char s[2048];
        float temp_freq[20];
        FILE *file = NULL;
        
        if ( FreqPath == NULL ){
                fprintf(stderr,"NULL pointer FreqPath %s %d  \n",__FILE__,__LINE__);
                return -1;
        }
        if ( (file = fopen(FreqPath, "r")) == NULL  ){
                fprintf(stderr,"NULL pointer opening %s %s %d  \n",FreqPath,__FILE__,__LINE__);
                return -1;
        }
        
        do {
                if (!fgets(s, sizeof(s), file)) {
                        fclose(file);
                        return -1;
                }
        } while ( 1 != sscanf(s, "%d", &temp_int));
        if (temp_int == 0) {
                fprintf(stderr,"Invalid Freqfile\n");
                fclose(file);
                return -1;
        }
        
        SetLen( temp_int );
        //read in the sequence
        fgets(s, sizeof(s), file);
        for (i = 0; i < this->len; i++) 
                this->residue[i].RS = s[i];
        
        //start reading in the frequencies
        i = 0;
        while (fgets (s, sizeof(s), file) && i < this->len) { 
                sscanf( s, "%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f",
                                &temp_freq[0], &temp_freq[1], &temp_freq[2], &temp_freq[3], &temp_freq[4], 
                                &temp_freq[5], &temp_freq[6], &temp_freq[7], &temp_freq[8], &temp_freq[9],
                                &temp_freq[10], &temp_freq[11], &temp_freq[12], &temp_freq[13], &temp_freq[14],
                                &temp_freq[15], &temp_freq[16], &temp_freq[17], &temp_freq[18], &temp_freq[19] );
                for (j = 0; j < 20; j++) {
                        this->residue[i].FREQ[j] = temp_freq[j];
                }
                this->residue[i].FREQ[20] = 0; //the X residue
                i++;
        }
        
        this->distconst_count = 0;
        
        
        return 0;
} 



int TargetStruct::Copy(TargetStruct *rhs, int start, int stop) {        
        if ( stop < 0)
                stop = rhs->len;
        if ( stop > rhs->len)
                stop = rhs->len;
        if (start < 0)
                start = 0;      
        
        SetLen( stop - start );
        for (int i = start; i < stop; i++) {
                this->residue[i-start] = rhs->residue[i];
        }
        for (int i = start; i < stop; i++) {
                this->structure[i-start] = rhs->structure[i];
        }       
        return 0;
}



int TargetStruct::Append(TargetStruct *rhs, int start, int stop) {      
        if ( stop < 0)
                stop = rhs->len;
        if ( stop > rhs->len)
                stop = rhs->len;
        if (start < 0)
                start = 0;      
        int orig_len = len;     
        TargetStruct tmp_targ;
        tmp_targ.Copy( this );  
        SetLen( orig_len + stop - start );      
        for ( int i = 0; i < orig_len; i++ ) {
                this->residue[i] = tmp_targ.residue[i];
                this->structure[i] = tmp_targ.structure[i];
        }       
        for (int i = start; i < stop; i++) {
                this->residue[orig_len + i-start] = rhs->residue[i];
                this->structure[orig_len + i-start] = rhs->structure[i];
        }       
        return 0;
}




int TargetStruct::LoadPDB(char *filename, char chain) {
        FILE *file;
        if (NULL == (file = fopen( filename, "r"))) {
                return 1;
        }       
        
        FreeMem();

        int tmp_len = 0;
        char  cur_residue_name[20];
        char buffer[2000];
        strcpy(cur_residue_name, "");   
        while (fgets(buffer, sizeof(buffer), file)) {
                char tmp_str[7];
                if ( ( !strncmp(buffer, "ATOM", 4) || !strncmp(buffer, "HETATM", 4))  && (chain == 0 || toupper(chain) == toupper(buffer[21]) ) ) {
                        strncpy( tmp_str, &(buffer[17]), 3 );
                        tmp_str[3] = 0;
                        char cur_RS = AA3toAA1( tmp_str );
                        if ( cur_RS != 'X' ) {                  
                                strncpy( tmp_str, &(buffer[22]), 5 );
                                tmp_str[5] = 0;
                                if (strcmp( cur_residue_name, tmp_str ) ) {
                                        tmp_len++;
                                        strcpy(cur_residue_name, tmp_str);
                                }
                        }
                }
                if ( !strncmp(buffer, "ENDMDL", 6) ) {
                        break;
                }
        }
        rewind(file);
        
        SetLen( tmp_len );
        
        for (int i = 0; i < len; i++) {
                this->structure[i].atom_count = 0;
                this->residue[i].RS = -1;
        }
        
        int res_num = -1;
        strcpy(cur_residue_name, "");   
        char tmp_x_str[30], tmp_y_str[30], tmp_z_str[30];
        strcpy(tmp_x_str, "");
        strcpy(tmp_y_str, "");
        strcpy(tmp_z_str, "");
        while (fgets(buffer, sizeof(buffer), file) ) {
                char tmp_str[50];
                if ( ( !strncmp(buffer, "ATOM", 4) || !strncmp(buffer, "HETATM", 4))  && (chain == 0 || toupper(chain) == toupper(buffer[21]) ) ) {
                        strncpy( tmp_str, &(buffer[17]), 3 );
                        tmp_str[3] = 0;
                        char cur_RS = AA3toAA1( tmp_str );
                        if ( cur_RS != 'X' ) {                  
                                strncpy( tmp_str, &(buffer[22]), 5 );
                                tmp_str[5] = 0;
                                if (strcmp( cur_residue_name, tmp_str ) ) {
                                        res_num++;
                                        strcpy( cur_residue_name, tmp_str );
                                }
                                if (res_num >= 0 && this->structure[res_num].atom_count < kMaxResAtoms) {
                                        this->structure[ res_num ].ReadAtomLine( buffer );                                      
                                        this->residue[res_num].RS = cur_RS;
                                }
                        }
                }
                if ( !strncmp(buffer, "ENDMDL", 6) ) {
                        break;
                }
        }
        fclose(file);
        
        //from here on out, the processing is done from the ATOM records that where read in
        //setup_residue_information( the_data );                                
        return 0;
}




/*
 int Prospect_Target_StructureSet(ProspectParam *param, 
                                                                  TemplateStruct *template_data, TargetStruct *target_data,
                                                                  AlignmentStruct *align_data, int *rotamer_selection ) {
         
         return 1;
 }
 */

int TargetStruct::CopyBackBone( TemplateStruct *template_data,
                                                                AlignmentStruct *align_data,
                                                                int start, int stop) {
        if ( stop < 0 || stop > this->len )
                stop = this->len;
        if (start < 0)
                start = 0;
        for (int i = start; i < stop; i++) {
                int j = align_data->GetTargetAlign(i);
                if (j != -1) {
                        this->structure[i] = template_data->structure[j];
                        this->residue[i].psi = template_data->residue[j].psi;
                        this->residue[i].phi = template_data->residue[j].phi;                   
                } else {
                        this->structure[i].atom_count = 0;
                }
        }
        return 0;
}




int TargetStruct::SideChain_Init(ProspectParam *param) {
        for (int i = 0; i < this->len; i++) {
                if ( this->structure[i].atom_count >= 4 ) {
                        int res_num = AA1toAANum( this->residue[i].RS);
                        if (res_num >= 0 && res_num < 20) {
                                ResidueInfoElement *res_info = &param->residue_info[ res_num ];
                                
                                assert( param->residue_info[ res_num ].structure.atom_count < kMaxResAtoms );
                                int k =  this->structure[i].atom_count;
                                assert( k < kMaxResAtoms );
                                
                                for (int j = 0; j < res_info->structure.atom_count; j++) {
                                        //strcpy( this->residue[i].structure.name[k], res_info->structure.name[j] );
                                        /*
                                         for (int l = 0; l < param->atom_data_count; l++) {
                                                 if ( this->residue[i].structure.atom_name[k] == PDBAtomtoNum( param->atom_data[l].name) ) {
                                                         this->residue[i].structure.radius[k] = param->atom_data[l].radius;
                                                 }
                                         }
                                         */
                                        this->structure[i].radius[k] = param->atom_data->FindRadius( this->structure[i].atom_name[k] );
                                        k++;
                                }
                                this->structure[i].atom_count = k;                              
                        }
                }
        }       
        return 0;
}


pValType calc_LRET(vec3 a, vec3 b, pValType ra, pValType rb) {
        pValType dist[3], d;
        pValType rad_dist = ra + rb;    
        dist[0] = a[0] - b[0];
        dist[1] = a[1] - b[1];
        dist[2] = a[2] - b[2];
        dist[0] = dist[0] * dist[0];
        dist[1] = dist[1] * dist[1];
        dist[2] = dist[2] * dist[2];
        d = dist[0] + dist[1] + dist[2];
        if ( d > rad_dist * rad_dist)
                return 0.0;
        d = sqrt(d);
        //assert( d != 0.0 );
        if ( d <= rad_dist*0.8254 )
                return 10;
        return 57.273 * (1 - d/rad_dist);       
}

pValType TargetStruct::Calc_LRET(ProspectParam *param) {
        pValType energy = 0.0;
        
        for (int i = 0; i < this->len; i++) {
                for (int ia = 0; ia < this->structure[i].atom_count; ia++) {
                        for (int ja = ia+1; ja < this->structure[i].atom_count; ja++) {
                                energy += calc_LRET( this->structure[i].atom[ia], this->structure[i].atom[ja],
                                                                         this->structure[i].radius[ia], this->structure[i].radius[ja]
                                                                         );
                        }
                        for (int j = i+1; j < this->len; j++) {
                                for (int ja = 0; ja < this->structure[j].atom_count; ja++) {
                                        energy += calc_LRET( this->structure[i].atom[ia], this->structure[j].atom[ja],
                                                                                 this->structure[i].radius[ia], this->structure[j].radius[ja]
                                                                                 );
                                }
                        }
                }                                       
        }
        return energy;  
}





int TargetStruct::SideChain_Pack(ProspectParam *param) {
        RotamerMask rotamer_mask;
        
        if (param->v_level > 0) 
                fprintf(stderr, "SideChain: Calculating energy\n");
        Prospect_Target_RotamerMask_Init(param, this, &rotamer_mask);
        if (param->v_level > 0) 
                fprintf(stderr, "SideChain: DEE\n");
        Prospect_RotamerMask_DEE(&rotamer_mask);
        if (param->v_level > 0) 
                fprintf(stderr, "SideChain: DEE_Split\n");
        Prospect_RotamerMask_DEE_Split(&rotamer_mask);
        
        
        //      Prospect_RotamerMask_DEE2(&rotamer_mask);
        if (param->v_level > 1) {
                fprintf(stderr, "Before tree decomp:\n");
                for (int i = 0; i < this->len; i++) {
                        fprintf(stderr, "Residue %d '%c' (%d of %ld rotamers)\n", 
                                        i, this->residue[i].RS, 
                                        rot_mask_active(&rotamer_mask, i), rotamer_mask.count[i]);
                }
        }
        if (param->v_level > 0) 
                fprintf(stderr, "SideChain: Tree decomposition\n");
        Prospect_RotamerMask_TreeDecomp(&rotamer_mask);
        
        
        //given the current elimination
        
        for (int i = 0; i < this->len; i++) {
                if ( this->structure[i].atom_count >= 4) {
                        int count = 0;
                        int a = 0;
                        for (int rot_a = 0; rot_a < rotamer_mask.count[i]; rot_a++) {
                                if ( rotamer_mask.state[i][rot_a] ) {
                                        if (count == 0)
                                                a = rot_a;
                                        count++;
                                }
                        }
                        if ( this->residue[i].RS == 'P' || this->residue[i].RS == 'G' || count == 0 ) {
                                Prospect_Target_Rotamer(param, &this->structure[i], this->residue[i].RS, 0);                            
                        } else {
                                Prospect_Target_Rotamer(param, &this->structure[i], this->residue[i].RS, rotamer_mask.rot_num[i][a]);
                        }
                        /*
                         fprintf(stderr, "Residue %d '%c' rot: %d (%d of %d rotamers) single: %f\n", 
                                         i, this->residue[i].RS, rotamer_mask.rot_num[i][a], count, rotamer_mask.count[i],
                                         rotamer_mask.single_score[a]);
                         
                         */
                }
        }               
        Prospect_Target_RotamerMask_Free(&rotamer_mask);        
        return 0;
}



int TargetStruct::Count_Atoms(void) {
        if ( structure == NULL )
                return 0;
        int count = 0;
        for (int i = 0; i < this->len; i++) {
                count += this->structure[i].atom_count;
        }
        return count;   
}

void TargetStruct::ReNumberAtoms(int start) {
        if ( structure == NULL )
                return;
        int cur_atom = start;   
        for ( int i = 0; i < len; i++) {
                for (int j = 0; j < structure[i].atom_count; j++) {
                        structure[i].atom_num[j] = cur_atom;
                        cur_atom++;
                }
        }
}


int TargetStruct::Count_SetResidues(void) {
        if ( structure == NULL )
                return 0;
        int count = 0;
        for (int i = 0; i < this->len; i++) {
                if ( this->structure[i].atom_count >= 4 ) 
                        count++;
        }
        return count;   
}


int TargetStruct::Count_UnSetResidues() {
        if ( structure == NULL )
                return len;
        int count = 0;
        for (int i = 0; i < this->len; i++) {
                if ( this->structure[i].atom_count < 4 ) 
                        count++;
        }
        return count;   
}


//BUG::need some way to describe residues that are not defined.
vec3 * TargetStruct::GetCaLocArray(void) {
        vec3 *data_out = (vec3 *)malloc(sizeof(vec3) * len );
        for ( int i = 0; i < len; i++ ) {
                if ( structure[i].bb_num[ anCa ] != -1 ) {
                        structure[i].GetBBLoc( anCa, data_out[i] );
                } else {
                        data_out[i][0] = 0.0 / 0.0;
                        data_out[i][1] = 0.0 / 0.0;
                        data_out[i][2] = 0.0 / 0.0;
                }
        }
        return data_out;        
}



int TargetStruct::FindResIndexFromResNum( int res_num, char res_mod ) {
        if ( structure == NULL )
                return -1;
        for ( int i = 0; i < len; i++ ) {
                if ( res_num == structure[i].res_num && res_mod == structure[i].res_mod ) {
                        return i;
                }
        }
        return -1;
}


/**********
Hydrogen bonding stuff
***********/


#define kDist_nh 1.0


void TargetStruct::Calc_hloc(int residue, vec3 h_loc) {
        float oc_dist, nca_dist, nh_dist, cah_dist;
        long k;
        
#if 1
        
        ResidueAtoms *residue_struct = &this->structure[ residue ];
        
        oc_dist = sqrt( pow( residue_struct->atom[ (int)residue_struct->bb_num[ anC ] ][0] - residue_struct->atom[ (int)residue_struct->bb_num[ anO ] ][0], 2) +
                                        pow( residue_struct->atom[ (int)residue_struct->bb_num[ anC ] ][1] - residue_struct->atom[ (int)residue_struct->bb_num[ anO ] ][1], 2) +
                                        pow( residue_struct->atom[ (int)residue_struct->bb_num[ anC ] ][2] - residue_struct->atom[ (int)residue_struct->bb_num[ anO ] ][2], 2) );
        
        for (k = 0; k < 3; k++) {
                h_loc[k] = residue_struct->atom[ (int)residue_struct->bb_num[ anN ] ][k] + 
                (residue_struct->atom[ (int)residue_struct->bb_num[ anC ] ][k] - residue_struct->atom[ (int)residue_struct->bb_num[ anO ] ][k]) / oc_dist;
        }
        
        
        cah_dist = sqrt( pow( residue_struct->atom[ (int)residue_struct->bb_num[ anCa ] ][0] - h_loc[0], 2) +
                                         pow( residue_struct->atom[ (int)residue_struct->bb_num[ anCa ] ][1] - h_loc[1], 2) +
                                         pow( residue_struct->atom[ (int)residue_struct->bb_num[ anCa ] ][2] - h_loc[2], 2) );
        
        nca_dist = sqrt( pow( residue_struct->atom[ (int)residue_struct->bb_num[ anN ] ][0] - residue_struct->atom[ (int)residue_struct->bb_num[ anCa ] ][0], 2) +
                                         pow( residue_struct->atom[ (int)residue_struct->bb_num[ anN ] ][1] - residue_struct->atom[ (int)residue_struct->bb_num[ anCa ] ][1], 2) +
                                         pow( residue_struct->atom[ (int)residue_struct->bb_num[ anN ] ][2] - residue_struct->atom[ (int)residue_struct->bb_num[ anCa ] ][2], 2) );
        
        
        nh_dist = sqrt( pow( residue_struct->atom[ (int)residue_struct->bb_num[ anN ] ][0] - h_loc[0], 2) +
                                        pow( residue_struct->atom[ (int)residue_struct->bb_num[ anN ] ][1] - h_loc[1], 2) +
                                        pow( residue_struct->atom[ (int)residue_struct->bb_num[ anN ] ][2] - h_loc[2], 2) );
        
        
        if ( 0 < (nca_dist*nca_dist + nh_dist * nh_dist - cah_dist * cah_dist) / (2*nca_dist*nh_dist) ) {       
                for (k = 0; k < 3; k++) {
                        h_loc[k] = residue_struct->atom[ (int)residue_struct->bb_num[ anN ] ][k] - 
                        (residue_struct->atom[ (int)residue_struct->bb_num[ anC ] ][k] - residue_struct->atom[ (int)residue_struct->bb_num[ anO ] ][k]) / oc_dist;
                }
        }
        
#else
        float nc_dist;
        
        nc_dist = sqrt( pow( pdb->atoms[ residue.n ].loc[0] - pdb->atoms[ residue.c ].loc[0], 2) +
                                        pow( pdb->atoms[ residue.n ].loc[1] - pdb->atoms[ residue.c ].loc[1], 2) +
                                        pow( pdb->atoms[ residue.n ].loc[2] - pdb->atoms[ residue.c ].loc[2], 2) );
        
        nca_dist = sqrt( pow( pdb->atoms[ residue.n ].loc[0] - pdb->atoms[ residue.ca ].loc[0], 2) +
                                         pow( pdb->atoms[ residue.n ].loc[1] - pdb->atoms[ residue.ca ].loc[1], 2) +
                                         pow( pdb->atoms[ residue.n ].loc[2] - pdb->atoms[ residue.ca ].loc[2], 2) );
        
        for (k = 0; k < 3; k++) {
                h_loc[k] = pdb->atoms[ residue.n ].loc[k] - 
                ( ( pdb->atoms[ residue.c ].loc[k] - pdb->atoms[ residue.n ].loc[k] ) / nc_dist +
                  ( pdb->atoms[ residue.ca ].loc[k] - pdb->atoms[ residue.n ].loc[k] ) / nca_dist );
        }
        
        /*
         nh_dist = sqrt( pow( h_loc[0] - pdb->atoms[ residue.n ].loc[0], 2) +
                                         pow( h_loc[1] - pdb->atoms[ residue.n ].loc[1], 2) +
                                         pow( h_loc[2] - pdb->atoms[ residue.n ].loc[2], 2) );
                                        
         for (k = 0; k < 3; k++) {
                 h_loc[k] = pdb->atoms[ residue.n ].loc[k] + 
                 kDist_nh * ( h_loc[k] - pdb->atoms[ residue.n ].loc[k] ) / nh_dist;
                 
         }
         */
        
#endif
        
}


#define q1     0.42
#define q2     0.20
#define f_val  332

float TargetStruct::Calc_hval(int i, int j) {
        vec3 h_loc;
        float on_dist, ch_dist, oh_dist, cn_dist, nh_dist;
        //float oc_dist, d_nc_dist, d_ch_dist, d_nh_dist;
        
        
        Calc_hloc( j, h_loc);
        /*
         oc_dist = sqrt( pow( pdb->atoms[ pdb->residues[j].c ].loc[0] - pdb->atoms[ pdb->residues[j].o ].loc[0], 2) +
                                         pow( pdb->atoms[ pdb->residues[j].c ].loc[1] - pdb->atoms[ pdb->residues[j].o ].loc[1], 2) +
                                         pow( pdb->atoms[ pdb->residues[j].c ].loc[2] - pdb->atoms[ pdb->residues[j].o ].loc[2], 2) );
         
         for (long k = 0; k < 3; k++) {
                 h_loc[k] = pdb->atoms[ pdb->residues[j].n ].loc[k] + 
                 (pdb->atoms[ pdb->residues[j].c ].loc[0] - pdb->atoms[ pdb->residues[j].o].loc[0]) / oc_dist;
                 
         }
         */
        
        ResidueAtoms *residue_i = &this->structure[i], 
                *residue_j = &this->structure[j];
        
        on_dist = sqrt( pow( residue_i->atom[ (int)residue_i->bb_num[ anO ]][0] - residue_j->atom[ (int)residue_j->bb_num[ anN ]][0], 2) +
                                        pow( residue_i->atom[ (int)residue_i->bb_num[ anO ]][1] - residue_j->atom[ (int)residue_j->bb_num[ anN ]][1], 2) +
                                        pow( residue_i->atom[ (int)residue_i->bb_num[ anO ]][2] - residue_j->atom[ (int)residue_j->bb_num[ anN ]][2], 2) );
        
        ch_dist = sqrt( pow( residue_i->atom[ (int)residue_i->bb_num[ anC ]][0] - h_loc[0], 2) +
                                        pow( residue_i->atom[ (int)residue_i->bb_num[ anC ]][1] - h_loc[1], 2) +
                                        pow( residue_i->atom[ (int)residue_i->bb_num[ anC ]][2] - h_loc[2], 2) );
        
        oh_dist = sqrt( pow( residue_i->atom[ (int)residue_i->bb_num[ anO ]][0] - h_loc[0], 2) +
                                        pow( residue_i->atom[ (int)residue_i->bb_num[ anO ]][1] - h_loc[1], 2) +
                                        pow( residue_i->atom[ (int)residue_i->bb_num[ anO ]][2] - h_loc[2], 2) );
        
        
        cn_dist = sqrt( pow( residue_i->atom[ (int)residue_i->bb_num[ anC ]][0] - residue_j->atom[ (int)residue_j->bb_num[ anN ]][0], 2) +
                                        pow( residue_i->atom[ (int)residue_i->bb_num[ anC ]][1] - residue_j->atom[ (int)residue_j->bb_num[ anN ]][1], 2) +
                                        pow( residue_i->atom[ (int)residue_i->bb_num[ anC ]][2] - residue_j->atom[ (int)residue_j->bb_num[ anN ]][2], 2) );
                                
        nh_dist = sqrt( pow( residue_i->atom[ (int)residue_i->bb_num[ anN ]][0] - h_loc[0], 2) +
                                        pow( residue_i->atom[ (int)residue_i->bb_num[ anN ]][1] - h_loc[1], 2) +
                                        pow( residue_i->atom[ (int)residue_i->bb_num[ anN ]][2] - h_loc[2], 2) );
        
        return q1 * q2 * ( 1 / on_dist + 1 / ch_dist - 1 / oh_dist - 1 / cn_dist) * f_val;
}




void TargetStruct::Setup_Hbonds( void ) {
        long i,j;
        float ca_ca_dist;
        //      float nc_dist;
        
        this->hbond = (float  **)malloc(sizeof(float*) * this->len );
        this->hbond[0] = (float  *)malloc(sizeof(float) * this->len * this->len );
        for (i = 1; i < this->len; i++) {
                this->hbond[i] = this->hbond[0] + this->len * i;        
        }
        
        for (i = 0; i < this->len; i++) {
                for (j = 0; j < this->len; j++) {
                        if ( abs(j-i) >= 3 ) {
                                ResidueAtoms *residue_i = &this->structure[i], 
                                *residue_j = &this->structure[j];
                                
                                ca_ca_dist = sqrt( pow( residue_i->atom[ (int)residue_i->bb_num[ anCa ] ][0] - residue_j->atom[ (int)residue_j->bb_num[ anCa ] ][0], 2) +
                                                                   pow( residue_i->atom[ (int)residue_i->bb_num[ anCa ] ][1] - residue_j->atom[ (int)residue_j->bb_num[ anCa ] ][1], 2) +
                                                                   pow( residue_i->atom[ (int)residue_i->bb_num[ anCa ] ][2] - residue_j->atom[ (int)residue_j->bb_num[ anCa ] ][2], 2) );
                                if ( ca_ca_dist > 9) {
                                        this->hbond[ i ][ j ] = 10.0;
                                } else {
                                        this->hbond[ i ][ j ] = Calc_hval(i, j);
                                }
                        } else {
                                this->hbond[ i ][ j ] = 10.0;
                        }
                }
        }       
}


int TargetStruct::DSSP(void) {
        if ( hbond == NULL)
                Setup_Hbonds();
        
        
        char *dssp_array = (char *)malloc(len + 7);
        memset( dssp_array, 0, len + 7 );
        
#define dsHCutoff            -0.5
#define dsParallelBridge     0x01
#define dsAntiParallelBridge 0x02
#define ds3Helix             0x04
#define ds4Helix             0x08
#define ds5Helix             0x10
#define ds3Turn              0x20
#define ds4Turn              0x40
#define ds5Turn              0x80
        
        bsheet = (char **)malloc( sizeof( char * ) * len );
        bsheet[0] = (char *)malloc( len * len );
        for (int i = 0; i < len; i++) {
                bsheet[i] = bsheet[0] + len * i;
                for (int j = 0; j < len; j++) {
                        bsheet[i][j] = 0;
                }
        }
        
        //parallel bridge scan
        for (int i = 1; i < len-1; i++) {
                for (int j = 0; j < len; j++) {
                        if ( hbond[i-1][j] < dsHCutoff && hbond[j][i+1] < dsHCutoff) {
                                dssp_array[i] |= dsParallelBridge;
                                dssp_array[j] |= dsParallelBridge;
                                bsheet[i][j] = 1;
                                bsheet[j][i] = 1;
                        }
                }
        }       
        
        //anti-parallel bridge scan
        for (int i = 0; i < len; i++) {
                for (int j = 0; j < len; j++) {
                        if ( hbond[i][j] < dsHCutoff && hbond[j][i] < dsHCutoff) {
                                dssp_array[i] |= dsAntiParallelBridge;
                                dssp_array[j] |= dsAntiParallelBridge;
                                bsheet[i][j] = 1;
                                bsheet[j][i] = 1;
                        }
                }
        }
        
        for (int i = 1; i < len-1; i++) {
                for (int j = 1; j < len-1; j++) {
                        if ( hbond[i-1][j+1] < dsHCutoff && hbond[j-1][i+1] < dsHCutoff) {
                                dssp_array[i] |= dsAntiParallelBridge;
                                dssp_array[j] |= dsAntiParallelBridge;
                                bsheet[i][j] = 1;
                                bsheet[j][i] = 1;
                        }
                }
        }       
        //3 turn search
        for (int i = 0; i < len-3; i++) {
                if (hbond[i][i+3] < dsHCutoff) {
                        dssp_array[i] |= ds3Turn;
                }
        }
        //4 turn search
        for (int i = 0; i < len-4; i++) {
                if (hbond[i][i+4] < dsHCutoff) {
                        dssp_array[i] |= ds4Turn;
                }
        }
        //5 turn search
        for (int i = 0; i < len-5; i++) {
                if (hbond[i][i+5] < dsHCutoff) {
                        dssp_array[i] |= ds5Turn;
                }
        }
        //3 helix search
        for (int i = 1; i < len ; i++) {
                if ( dssp_array[i-1] & ds3Turn && dssp_array[i] & ds3Turn ) {
                        dssp_array[i] |= ds3Helix;
                        dssp_array[i+1] |= ds3Helix;
                        dssp_array[i+2] |= ds3Helix;
                }
        }
        //4 helix search
        for (int i = 1; i < len ; i++) {
                if ( dssp_array[i-1] & ds4Turn && dssp_array[i] & ds4Turn ) {
                        dssp_array[i] |= ds4Helix;
                        dssp_array[i+1] |= ds4Helix;
                        dssp_array[i+2] |= ds4Helix;
                        dssp_array[i+3] |= ds4Helix;
                }
        }
        //5 helix search
        for (int i = 1; i < len ; i++) {
                if ( dssp_array[i-1] & ds5Turn && dssp_array[i] & ds5Turn ) {
                        dssp_array[i] |= ds5Helix;
                        dssp_array[i+1] |= ds5Helix;
                        dssp_array[i+2] |= ds5Helix;
                        dssp_array[i+3] |= ds5Helix;
                        dssp_array[i+4] |= ds5Helix;
                        dssp_array[i+5] |= ds5Helix;
                }
        }
        for (int i = 0; i < len; i++) {
                if ( dssp_array[i] & dsParallelBridge || dssp_array[i] & dsAntiParallelBridge ) {
                        residue[i].SS = 'E';
                } else if ( dssp_array[i] & ds3Helix || dssp_array[i] & ds4Helix || dssp_array[i] & ds5Helix ) {
                        residue[i].SS = 'H';
                } else {
                        residue[i].SS = 'C';
                }               
        }
        
        free(dssp_array);
        return 0;
}




#define D_CB_CA 1.538
#define A_CB_CA_N 113.5
#define A_CB_CA_C 108

//hypothetical Cb calculation for glycine
int ResidueAtoms_CalcCb(ResidueAtoms *res) {
        /* Given coordinates of C, N, CA, and O atoms, calculate
        CB coordinates.
        Assume O and CB are at different sides of C-N-CA plane.
        */
        double n2, c2, ca2; /* |n|^2, |c|^2, |ca|^2 */
        double d_ca_n, d_ca_c; /* distance of CA-N and CA-C */
        double ntemp, ctemp; /*intermediate variables*/
        double xz, xz0; /*more intermediate variables*/
        double yz, yz0; /*more intermediate variables*/
        double zz, zz0; /*more intermediate variables*/
        double x1,y1,z1,x2,y2,z2,z,d1,d2; /*two sets of solutions of CB*/
        double tx=19.491, ty=5.806, tz=22.004;
        int Gly_flag = 1;
        
        
        vec3 n, o, c, ca, cb;
        
        if (  res->bb_num[ anN ] == -1 ||  res->bb_num[ anO ] == -1 ||  res->bb_num[ anC ] == -1 ||  res->bb_num[ anCa ] == -1 ) {
                return 1;
        }
        
        n[0] = res->atom[ res->bb_num[ anN ] ][0];
        n[1] = res->atom[ res->bb_num[ anN ] ][1];
        n[2] = res->atom[ res->bb_num[ anN ] ][2];
        o[0] = res->atom[ res->bb_num[ anO ] ][0];
        o[1] = res->atom[ res->bb_num[ anO ] ][1];
        o[2] = res->atom[ res->bb_num[ anO ] ][2];
        c[0] = res->atom[ res->bb_num[ anC ] ][0];
        c[1] = res->atom[ res->bb_num[ anC ] ][1];
        c[2] = res->atom[ res->bb_num[ anC ] ][2];
        ca[0] = res->atom[ res->bb_num[ anCa ] ][0];
        ca[1] = res->atom[ res->bb_num[ anCa ] ][1];
        ca[2] = res->atom[ res->bb_num[ anCa ] ][2];
        
        
        n2 = n[0]*n[0]+n[1]*n[1]+n[2]*n[2];
        c2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
        ca2 = ca[0]*ca[0]+ca[1]*ca[1]+ca[2]*ca[2];
        
        d_ca_n= (ca[0]-n[0])*(ca[0]-n[0])+(ca[1]-n[1])*(ca[1]-n[1])+(ca[2]-n[2])*(ca[2]-n[2]);
        d_ca_n= sqrt(d_ca_n);
        d_ca_c= (ca[0]-c[0])*(ca[0]-c[0])+(ca[1]-c[1])*(ca[1]-c[1])+(ca[2]-c[2])*(ca[2]-c[2]);
        d_ca_c= sqrt(d_ca_c);
        
        ntemp = d_ca_n*d_ca_n-2*d_ca_n*D_CB_CA*cos(M_PI*A_CB_CA_N/180.0);
        ntemp = (ntemp+ca2-n2)/2;
        
        ctemp = d_ca_c*d_ca_c -2*d_ca_c*D_CB_CA*cos(M_PI*A_CB_CA_C/180.0);
        ctemp = (ctemp+ca2-c2)/2;
        
        yz = (ca[0]-c[0])*(ca[1]-n[1])-(ca[0]-n[0])*(ca[1]-c[1]);
        if(yz==0) {
                Gly_flag = 0;
                //      printf("Unsolvable GLY, yz=%f\n",yz);
        }
        yz0 = (ntemp*(ca[0]-c[0])-ctemp*(ca[0]-n[0]))/yz;
        yz = ((n[2]-ca[2])*(ca[0]-c[0])-(c[2]-ca[2])*(ca[0]-n[0]))/yz;
        
        xz0 = (ntemp-yz0*(ca[1]-n[1]))/(ca[0]-n[0]);
        xz = ((n[2]-ca[2])-yz*(ca[1]-n[1]))/(ca[0]-n[0]);
        if(ca[0]==n[0]) {
                //      printf("Unsolvable GLY, ca[0]=n[0]=%f\n",n[0]);
                Gly_flag = 0;
        }
        
        z = 1+xz*xz+yz*yz;
        zz = -(xz*(xz0-ca[0])+yz*(yz0-ca[1])-ca[2])/z;
        zz0 = (xz0-ca[0])*(xz0-ca[0])+(yz0-ca[1])*(yz0-ca[1])+ca[2]*ca[2]-D_CB_CA*D_CB_CA;
        
        zz0 = (xz*(xz0-ca[0])+yz*(yz0-ca[1])-ca[2])*(xz*(xz0-ca[0])+yz*(yz0-ca[1])-ca[2])-z*zz0;
        
        if(zz0<0) {
                Gly_flag = 0;
                //      printf("Unsolvable GLY, zz0=%f\n",zz0);
        }
        zz0 = sqrt(zz0)/z;
        
        /*Two sets of solutions */
        z1 = zz+zz0;
        z2 = zz-zz0;
        x1 = xz0+xz*z1;
        x2 = xz0+xz*z2;
        y1 = yz0+yz*z1;
        y2 = yz0+yz*z2;
        
        /* Compare the two for the distance of CB-O, choose the longer one */
        
        d1 = (x1-o[0])*(x1-o[0])+(y1-o[1])*(y1-o[1])+(z1-o[2])*(z1-o[2]);
        d2 = (x2-o[0])*(x2-o[0])+(y2-o[1])*(y2-o[1])+(z2-o[2])*(z2-o[2]);
        
        if(d1>d2){
                cb[0]=x1;
                cb[1]=y1;
                cb[2]=z1;
        } else {
                cb[0]=x2;
                cb[1]=y2;
                cb[2]=z2;
        }       
        if (Gly_flag == 0) {
                x1=(n[0]+c[0])/2;
                y1=(n[1]+c[1])/2;
                z1=(n[2]+c[2])/2;
                x1=ca[0]-x1;
                y1=ca[1]-y1;
                z1=ca[2]-z1;
                zz=sqrt(x1*x1+y1*y1+z1*z1);
                x1=x1*D_CB_CA/zz;
                y1=y1*D_CB_CA/zz;
                z1=z1*D_CB_CA/zz;
                cb[0]=ca[0]+x1;
                cb[1]=ca[1]+y1;
                cb[2]=ca[2]+z1;
        }       
        if ( res->bb_num[ anCb ] == -1 ) {
                res->bb_num[ anCb ] = res->atom_count;
                res->atom_count++;
        }
        res->atom[ res->bb_num[ anCb ] ][0] = cb[0];
        res->atom[ res->bb_num[ anCb ] ][1] = cb[1];
        res->atom[ res->bb_num[ anCb ] ][2] = cb[2];
        return 0;
}





char * TargetStruct::WritePDBText(int start, int stop) {        
        if ( stop < 0 || stop > this->len)
                stop = this->len;
        if (start < 0)
                start = 0;
        
        char *outtext = NULL;
        char buffer[ 4000 ];    
        int atom_count = 1;
        for (int i = start; i < stop; i++) {
        //      if (this->residue[i].structure_state != rsUndefined) {

//              char tmp[5];
//                      AA1toAA3(this->residue[i].RS, tmp);     
                for (int j = 0; j < this->structure[i].atom_count; j++) {
                
                        /*
                        if ( this->structure[i].atom_num[j] >= 0 )
                                        atom_count = this->structure[i].atom_num[j];
                                sprintf(buffer, "ATOM  %5d  %-3s %+3s %c%5d   %8.3f%8.3f%8.3f  1.00%6.2f\n", 
                                                atom_count, 
                                                pdb_atomname[ this->structure[i].atom_name[j] ], 
                                                tmp, 
                                                this->residue[i].chain, 
                                                i + residue_offset,
                                                this->structure[i].atom[j][0], 
                                                this->structure[i].atom[j][1], 
                                                this->structure[i].atom[j][2],
                                                this->structure[i].temp_factor[j] );
                                atom_count++;
*/
                        char *tmp_text = this->structure[i].WriteAtomLine( j, this->residue[i].RS, i, this->residue[i].chain );
                        if ( tmp_text ) {
                                append_to_str(&outtext, tmp_text);                      
                                free( tmp_text );
                        }
                }
        }
        return outtext; 
}


int TargetStruct::WritePDB(char *outpath, int start, int stop) {
        FILE *file = NULL;
        
        if (outpath == NULL) {
                file = stdout;
        } else {
                file=fopen(outpath, "w");
        }
        if (file == NULL) {
                return 1;
        }       
        char *text = WritePDBText(start, stop);
        fprintf(file, "%s", text);
        free(text);
        if(file != stdout) {
                fclose(file);
        }
        return 0;
}



void TargetStruct::CleanUndefinedResidues(void) {
        
        if ( structure == NULL )
                return;
        
        int new_len = 0;
        for ( int i = 0; i < len; i++ ) {
                if ( structure[i].bb_num[ anCa ] != -1 && structure[i].bb_num[ anN ] != -1 &&
                         structure[i].bb_num[ anN ] != -1  && structure[i].bb_num[ anO ] != -1 ) {
                        new_len++;                      
                }
        }       
        if ( new_len != len ) {         
                ResidueAtoms  *old_structure = new ResidueAtoms[ len ];
                TargetResidue *old_residue = new TargetResidue[ len ];
                for ( int i = 0; i < len; i++ ) {
                        old_structure[i] = structure[i];
                        old_residue[i]   = residue[i];  
                }
                int old_len = len;
                SetLen(new_len);
                int cur_res = 0;
                for ( int i = 0; i < old_len; i++ ) {
                        if ( old_structure[i].bb_num[ anCa ] != -1 && old_structure[i].bb_num[ anN ] != -1 &&
                                 old_structure[i].bb_num[ anN ] != -1  && old_structure[i].bb_num[ anO ] != -1 ) {
                                structure[cur_res] = old_structure[ i ];
                                residue[cur_res]   = old_residue[ i ];                          
                                cur_res++;                      
                        }
                }               
        }       
}



char *TargetStruct::GetSeq(void) {
        char *out = (char *)malloc(sizeof(char) * (len+1));
        for (int i = 0; i < len; i++) {
                out[i] = residue[i].RS;
        }
        out[ len ] = 0;
        return out;
}


float choufasman_a[20] = {
        142,
        98,
        101,
    67,
        70,
        151 ,
    111 ,
        57 ,
        100 ,
        108  ,
        121 ,
        114  ,
        145 ,
        113  ,
        57   ,
        77 ,
        83 ,
        108  ,
        69  ,
        106     
};

float choufasman_b[20] = {
        83,
        93,
        54,
        89,
        119,
        037,
        110,
        75,
        87,
        160,
        130,
        74,
        105,
        138,
        55,
        75,
        119,
        137,
        147,
        170
};

float choufasman_t[20] = {
        66,
        95,
        146,
        156,
        119,
        74,
        98,
        156,
        95,
        47,
        59,
        101,
        60,
        60,
        152,
        143,
        96,
        96,
        114,
        50      
};

float choufasman_f[20][4] = {
        {0.06  ,  0.076,   0.035,   0.058},
    {0.070,   0.106,   0.099,   0.085},
    {0.147,   0.110,   0.179,  0.081},
        {0.161,   0.083,   0.191,   0.091},
        {0.149,   0.050,   0.117,   0.128},
        {0.056,   0.060,   0.077,   0.064},
    {0.074,   0.098,   0.037,   0.098},
        {0.102,   0.085,   0.190,   0.152},
        {0.140,   0.047,   0.093,   0.054},
        {0.043,   0.034,   0.013,   0.056},
        {0.061,   0.025,   0.036 ,  0.070},
        {0.055,   0.115,   0.072,   0.095},
        {0.068,   0.082,   0.014,   0.055},
        {0.059,   0.041,   0.065,   0.065},
        {0.102,   0.301,   0.034,   0.068},
        {0.120,   0.139,   0.125,   0.106},
        {0.086,   0.108,   0.065,   0.079},
        {0.077,   0.013,   0.064,   0.167},
        {0.082,   0.065,   0.114,   0.125},
        {0.062,   0.048,   0.028,   0.053}
};

void TargetStruct::ChouFasman(void) {
        
        for ( int i = 0; i < len-6; i++) {
                int a_hit_count = 0;
                for ( int j = 0; j < 6; j++ ) {
                        if ( choufasman_a[ AA1toAANum( residue[i+j].RS ) ] > 100 ) {
                                a_hit_count++;
                        }
                }
                if ( a_hit_count >= 4 ) {
                        //possible hit, extend in both directions
                }
                
                
                
                
        }
}


void TargetStruct::MoveMeanToOrigin(void) {
        if ( structure == NULL )
                return;
        vec3 mean = {0,0,0};
        int atom_count = 0;
        for (int i = 0; i < len; i++) {
                ResidueAtoms *res = &structure[i];
                for ( int j = 0; j < res->atom_count; j++) {
                        mean[0] += res->atom[j][0];
                        mean[1] += res->atom[j][1];
                        mean[2] += res->atom[j][2];
                        atom_count++;
                }
        }
        mean[0] /= (double)atom_count;
        mean[1] /= (double)atom_count;
        mean[2] /= (double)atom_count;
        for (int i = 0; i < len; i++) {
                ResidueAtoms *res = &structure[i];
                for ( int j = 0; j < res->atom_count; j++) {
                        res->atom[j][0] -= mean[0];
                        res->atom[j][1] -= mean[1];
                        res->atom[j][2] -= mean[2];
                }
        }
}



void TargetStruct::Translate(vec3 trans) {
        if (structure == NULL) return;
        for (int i = 0; i < len; i++) {
                smAdd(trans, structure[i].atom, structure->atom_count);
        }
}

void TargetStruct::Rotate (mat3 rotate){
        if (structure == NULL) return;
        for (int i = 0; i < len; i++) {
                smMultMat(rotate, structure[i].atom, structure->atom_count);
        }       
}

//BUG:: if a residue doesn't have a defined Ca atom, this could get ugly...
int TargetStruct::GetCaCount(void) {
        if ( structure == NULL ) return 0;
        return len;
}

vec3 *TargetStruct::GetCaArray(void) {
        if ( structure == NULL ) return NULL;
        vec3* array = (vec3 *)malloc(sizeof(vec3) * len );
        for (int i = 0; i < len; i++) {
                smCopy( structure[i].atom[ structure[i].bb_num[ anCa ] ], array[i] );
        }
        return array;
}

---