/*
FILE:     struct_align.C
*/
/*
VERSION:  1.600
*/
/*
DATE:     2/3/2005
*/
/*
  Comments and Questions to: sw-help@rcsb.rutgers.edu
*/
/*
COPYRIGHT 1999-2005 Rutgers - The State University of New Jersey

This software is provided WITHOUT WARRANTY OF MERCHANTABILITY OR
FITNESS FOR A PARTICULAR PURPOSE OR ANY OTHER WARRANTY, EXPRESS OR
IMPLIED.  RUTGERS MAKE NO REPRESENTATION OR WARRANTY THAT THE
SOFTWARE WILL NOT INFRINGE ANY PATENT, COPYRIGHT OR OTHER
PROPRIETARY RIGHT.

The user of this software shall indemnify, hold harmless and defend
Rutgers, its governors, trustees, officers, employees, students,
agents and the authors against any and all claims, suits,
losses, liabilities, damages, costs, fees, and expenses including
reasonable attorneys' fees resulting from or arising out of the
use of this software.  This indemnification shall include, but is
not limited to, any and all claims alleging products liability.
*/
/*
               PDB SOFTWARE LICENSE AGREEMENT

BY CLICKING THE ACCEPTANCE BUTTON OR INSTALLING OR USING 
THIS "SOFTWARE, THE INDIVIDUAL OR ENTITY LICENSING THE  
SOFTWARE ("LICENSEE") IS CONSENTING TO BE BOUND BY AND IS 
BECOMING A PARTY TO THIS AGREEMENT.  IF LICENSEE DOES NOT 
AGREE TO ALL OF THE TERMS OF THIS AGREEMENT
THE LICENSEE MUST NOT INSTALL OR USE THE SOFTWARE.

1. LICENSE AGREEMENT

This is a license between you ("Licensee") and the Protein Data Bank (PDB) 
at Rutgers, The State University of New Jersey (hereafter referred to 
as "RUTGERS").   The software is owned by RUTGERS and protected by 
copyright laws, and some elements are protected by laws governing 
trademarks, trade dress and trade secrets, and may be protected by 
patent laws. 

2. LICENSE GRANT

RUTGERS grants you, and you hereby accept, non-exclusive, royalty-free 
perpetual license to install, use, modify, prepare derivative works, 
incorporate into other computer software, and distribute in binary 
and source code format, or any derivative work thereof, together with 
any associated media, printed materials, and on-line or electronic 
documentation (if any) provided by RUTGERS (collectively, the "SOFTWARE"), 
subject to the following terms and conditions: (i) any distribution 
of the SOFTWARE shall bind the receiver to the terms and conditions 
of this Agreement; (ii) any distribution of the SOFTWARE in modified 
form shall clearly state that the SOFTWARE has been modified from 
the version originally obtained from RUTGERS.  

2. COPYRIGHT; RETENTION OF RIGHTS.  

The above license grant is conditioned on the following: (i) you must 
reproduce all copyright notices and other proprietary notices on any 
copies of the SOFTWARE and you must not remove such notices; (ii) in 
the event you compile the SOFTWARE, you will include the copyright 
notice with the binary in such a manner as to allow it to be easily 
viewable; (iii) if you incorporate the SOFTWARE into other code, you 
must provide notice that the code contains the SOFTWARE and include 
a copy of the copyright notices and other proprietary notices.  All 
copies of the SOFTWARE shall be subject to the terms of this Agreement.  

3. NO MAINTENANCE OR SUPPORT; TREATMENT OF ENHANCEMENTS 

RUTGERS is under no obligation whatsoever to: (i) provide maintenance 
or support for the SOFTWARE; or (ii) to notify you of bug fixes, patches, 
or upgrades to the features, functionality or performance of the 
SOFTWARE ("Enhancements") (if any), whether developed by RUTGERS 
or third parties.  If, in its sole discretion, RUTGERS makes an 
Enhancement available to you and RUTGERS does not separately enter 
into a written license agreement with you relating to such bug fix, 
patch or upgrade, then it shall be deemed incorporated into the SOFTWARE 
and subject to this Agreement. You are under no obligation whatsoever 
to provide any Enhancements to RUTGERS or the public that you may 
develop over time; however, if you choose to provide your Enhancements 
to RUTGERS, or if you choose to otherwise publish or distribute your 
Enhancements, in source code form without contemporaneously requiring 
end users or RUTGERS to enter into a separate written license agreement 
for such Enhancements, then you hereby grant RUTGERS a non-exclusive,
royalty-free perpetual license to install, use, modify, prepare
derivative works, incorporate into the SOFTWARE or other computer
software, distribute, and sublicense your Enhancements or derivative
works thereof, in binary and source code form.

4. FEES.  There is no license fee for the SOFTWARE.  If Licensee
wishes to receive the SOFTWARE on media, there may be a small charge
for the media and for shipping and handling.  Licensee is
responsible for any and all taxes.

5. TERMINATION.  Without prejudice to any other rights, Licensor
may terminate this Agreement if Licensee breaches any of its terms
and conditions.  Upon termination, Licensee shall destroy all
copies of the SOFTWARE.

6. PROPRIETARY RIGHTS.  Title, ownership rights, and intellectual
property rights in the Product shall remain with RUTGERS.  Licensee 
acknowledges such ownership and intellectual property rights and will 
not take any action to jeopardize, limit or interfere in any manner 
with RUTGERS' ownership of or rights with respect to the SOFTWARE.  
The SOFTWARE is protected by copyright and other intellectual 
property laws and by international treaties.  Title and related 
rights in the content accessed through the SOFTWARE is the property 
of the applicable content owner and is protected by applicable law.  
The license granted under this Agreement gives Licensee no rights to such
content.

7. DISCLAIMER OF WARRANTY.  THE SOFTWARE IS PROVIDED FREE OF 
CHARGE, AND, THEREFORE, ON AN "AS IS" BASIS, WITHOUT WARRANTY OF 
ANY KIND, INCLUDING WITHOUT LIMITATION THE WARRANTIES THAT IT 
IS FREE OF DEFECTS, MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE 
OR NON-INFRINGING.  THE ENTIRE RISK AS TO THE QUALITY AND 
PERFORMANCE OF THE SOFTWARE IS BORNE BY LICENSEE.  SHOULD THE 
SOFTWARE PROVE DEFECTIVE IN ANY RESPECT, THE LICENSEE AND NOT 
LICENSOR ASSUMES THE ENTIRE COST OF ANY SERVICE AND REPAIR.  
THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL PART OF 
THIS AGREEMENT.  NO USE OF THE PRODUCT IS AUTHORIZED HEREUNDER 
EXCEPT UNDER THIS DISCLAIMER.

8. LIMITATION OF LIABILITY.  TO THE MAXIMUM EXTENT PERMITTED BY
APPLICABLE LAW,  IN NO EVENT WILL LICENSOR BE LIABLE FOR ANY 
INDIRECT, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING 
OUT OF THE USE OF OR INABILITY TO USE THE SOFTWARE, INCLUDING, 
WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK 
STOPPAGE, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL 
OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF ADVISED OF THE
POSSIBILITY THEREOF. 
*/
/* the input coordinate file either is PDF format or CIF */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <unistd.h>
#include <sys/types.h>

#include "nrutil.h"
#include "util.h"
#include "match.h"
#include "sg_nsym.h"
#include "atom_mass.h"
#include "version.h"

char METHOD[10];  //X-RAY or NMR

#define  NMISC 29
#define MAXL 512
void write_entity_poly_text(FILE *fadit_in,  FILE *fadit_out);
void parse_adit_text_to_mmcif(FILE *fout, char *inpfile);
void get_pdb_one_model(char *inpfile,  char *outfile);
float single_residue_mass(char *residue_type, char residue_name);
float residue_mass(char *seq, char *chain_id, char *residue_type);
void write_matth_coeff_to_data_template( char *inpfile);
int get_symmetry_num(char *sg_name_inp, char *sg_name_out);
float get_cell_volume(float *cell);
void values_in_bracket(FILE *inp, char *first_item, char *value, int *key);
void output_info();


char *getenv();

struct tmp{
    char space_group[20];
    char a[20];
    char b[20];
    char c[20];
    char alpha[20];
    char beta[20];
    char gamma[20];
} cell;

FILE *fout;

int main(int argc, char **argv)
{
    char inpfile[MAXL],  inpfile_cif[MAXL];
    int pid, i;
    
    
    if(argc == 1 || (argc ==2 && !strcmp_case(argv[1], "-h"))){
        usage();
    }
    pid=getpid();
    
    for(i=0; i<argc-1; i++){
      if(!strcmp_case(argv[i], "-NMR")){
          strcpy(METHOD,"NMR");
          break;
      }
    }

    for(i=0; i<argc-1; i++){
          
        if(!strcmp_case(argv[i], "-PDB") || !strcmp_case(argv[i], "-IPDB")){
            
            
            if(is_pdb(argv[i+1])<20){
                printf("Warning! The file (%s) is not a protein PDB file.\n",argv[i+1]);
            }
            
            if(!strcmp_case(METHOD, "NMR")){
                sprintf(inpfile, "nmr-%s",argv[i+1]);
                get_pdb_one_model(argv[i+1], inpfile);
                process_file(inpfile, "PDB");
                delete_file(inpfile);
            }else{
                process_file(argv[i+1],"PDB" );
            }
            
            output_info();
            
        
        }else if(!strcmp_case(argv[i], "-CIF")|| !strcmp_case(argv[i], "-ICIF")){

            sprintf(inpfile_cif, "%s%d.pdb",argv[i+1], pid);
            cif2pdb(argv[i+1], inpfile_cif);
            process_file(inpfile_cif, "CIF");
            delete_file(inpfile_cif);
            output_info();
    
        }else if(!strcmp_case(argv[i], "-EXT")){
            strcpy(inpfile,argv[i+1]);
            get_mmcif(inpfile);
        
        }else if (!strcmp_case(argv[i], "-SOL")){
            strcpy(inpfile,argv[i+1]);
            write_matth_coeff_to_data_template(inpfile);
            printf("\nThe new file (%s) has been generated\n",inpfile);
            
        }
        
        
    }
//    delete_file("entity-poly.mmcif");
//    delete_file("cgi_value");  // keep this when using  web page.
    
    
}
void output_info()
{
    printf("\n\nOutput file names:  data_template.text and log_script.inp\n\n");
            
    printf("To get the complete structure information, \n");
    printf("enter necessary information to data_template.text and log_script.inp.\n");
    printf("EITHER  type the command 'extract  -EXT log_script.inp',\n");
    printf("OR  use pdb_extract (type pdb_extract -h for help).\n");
}

void get_pdb_one_model(char *inpfile,  char *outfile)
/* extract the nmr coordinate. If found the best, write the best to
   parfile_nmr.pdb, if not , write the first model to parfile_nmr.pdb
*/
{
    long  j,  model;
    char str[256],substr[80];
    FILE *fp, *fout;
    
    
    if((fp = fopen(inpfile, "r"))==NULL)
        open_file_error(inpfile, "get_pdb_one_model");
    fout = fopen(outfile, "w");
    

    model = -999;
    while (fgets(str, sizeof str, fp) != NULL) {
        upperstr(str);
        if(!strncmp(str, "REMARK", 6) &&
           strstr(str, "BEST REPRESENTATIVE CONFORMER IN THIS ENSEMBLE :") ){
            strcpy(substr, strstr(str,":")+1);
            if(sscanf(substr, "%ld", &j)==1) {
                model=j;
            }else
                model = -999;
        }
        if(!strncmp(str, "ATOM", 4)  || !strncmp(str, "HETA", 4) )
            break;
    }

    rewind(fp);
    if(model == -999){
        fprintf(fout,"REMARK   NMR structrue:    model = 1\n"); 
        while (fgets(str, sizeof str, fp) != NULL) {
            upperstr(str);
            if(!strncmp(str, "ATOM", 4)  || !strncmp(str, "HETA", 4) ){
                fprintf(fout,"%s", str);
                while (fgets(str, sizeof str, fp) != NULL) {
                    if(!strncmp(str, "MODEL", 5) || !strncmp(str, "ENDMDL", 6)
                       || !strncmp(str, "END", 3)){
                        fprintf(fout,"END\n");
                        fclose(fout);
                        return;
                    }
                    fprintf(fout,"%s", str);
                }
            }
        }
    }else{
        j=-999;
        while (fgets(str, sizeof str, fp) != NULL) {
            if(!strncmp(str, "MODEL", 5) && sscanf(str, "%*s %ld", &j)==1){
                if(j == model) break;
            }
        }
        if(j != -999){ /*found model*/
            fprintf(fout,"REMARK  NMR structrue: (the best model = %ld ).\n",j); 
                        
            while (fgets(str, sizeof str, fp) != NULL) {
                if(!strncmp(str, "MODEL", 5) || !strncmp(str, "ENDMDL", 6)
                   || !strncmp(str, "END", 3)){
                    fprintf(fout,"END\n");
                    fclose(fout);
                    return;
                }
                fprintf(fout,"%s", str);
            }
        }else{ /* not found model*/
            rewind(fp);
            while (fgets(str, sizeof str, fp) != NULL) {
                fprintf(fout,"%s", str);
            }
        }
                
    }
    
}
    

void usage(void)
{
    printf("\nUsage: extract [option] input_file\n");
    printf("-EXPT  followed by experiment method (NMR or XRAY (defult)).\n");
    printf("-PDB  followed by a coordinate file in PDB format.\n");
    printf("-CIF  followed by a coordinate file in mmCIF format.\n");
    printf("-EXT  followed by the file log_script.inp (generated).\n\n");
    
    printf("EXAMPLE for X-RAY structures:\n");
    printf("  Type 'extract  -PDB  pdb_file_name' \n");
    printf("EXAMPLE for NMR structures:\n");
    printf("  Type 'extract  -PDB  pdb_file_name -EXPT NMR' \n\n");
    
    printf("Enter necessary information to data_template.text and log_script.inp\n");
    printf("\nEITHER\n");
    printf("use 'extract  -EXT log_script.inp ' to get the complete mmCIF output.\n");
    printf("OR\n");
    printf("use pdb_extract (type pdb_extract -h for help) to get the complete mmCIF output.\n\n");
    
    exit(0);
}


   
void process_file(char *pdbfile, const char *format)
/* do all sorts of calculations */
{
    char outfile[MAXL],ALT_LIST[MAXL],chain_nam1,chain_nam2, tmp[MAXL];
    char *ChainID, *bseq, **AtomName, **ResName, **Miscs;
    char **chain_res, **chain_res_new, *matched,matched_comp[1000];
    char *temp=NULL, adit_text[MAXL],adit_inp[MAXL], str[MAXL],text_path[MAXL];
    
    long i=0,j=0,j1=0, k=0, k1=0, k2=0, num=0,natom=0, num_residue=0;
    long large=0, len1=0, len2=0, ntot=0, half_resd=0, nchain=0, npoly=0;
    long *ResSeq, *RY, **seidx, **chain_idx, **broken_chain_idx,*nbroken;
    long **seq1_match, **seq2_match, nmatch;
    double **xyz, rmsd,rmsd_new, percent;
    FILE  *fcgi=NULL, *mmcif=NULL, *fadit_inp_in=NULL, *fadit_inp_out=NULL;    
    FILE  *fadit_text_in=NULL, *fadit_text_out=NULL;
    
//    fp = fopen(pdbfile, "r");
    
    sprintf(outfile, "%s.out", pdbfile);
    fout=fopen(outfile, "w");
    fcgi=fopen("cgi_value", "w"); /* export a file for cgi (web page)*/
    fadit_text_out=fopen("data_template.text", "w");
    fadit_inp_out=fopen("log_script.inp", "w");
    mmcif=fopen("entity-poly.mmcif", "w");

    fprintf(mmcif,"_entity_poly.entity_id\n");
    fprintf(mmcif,"_entity_poly.type\n");
    fprintf(mmcif,"_entity_poly.pdbx_seq_one_letter_code\n");
    fprintf(mmcif,"_entity_poly.pdbx_strand_id\n");
    
/* read in the PDB file and get the number of atoms*/
    num = number_of_atoms(pdbfile);
    if(num<=2){
        printf("Error! please check the input option or the input file\n");
        return;
    }
  
    AtomName = cmatrix(0, num, 0, 4);
    ResName = cmatrix(0, num, 0, 3);
    ChainID = cvector(0, num);
    ResSeq = lvector(0, num);
    xyz = dmatrix(0, num, 0, 3);
    Miscs = cmatrix(0, num, 0, NMISC);

    seq1_match = lmatrix(0,40, 0,3);
    seq2_match = lmatrix(0,40, 0,3);
    
        
    printf("\nPDB file name: %s\n\n",  pdbfile);
    fprintf(fout,"PDB data file name: %s\n",  pdbfile);

    read_pdb_head(pdbfile);

    strcpy(adit_text, "");
    strcpy(adit_inp, "");


    
    (!strcmp(CCP4,"yes")) ? temp = getenv("CCP4"):temp = getenv("RCSBROOT");
        
    if(temp == NULL) {
        if(!strcmp(CCP4,"yes"))
            printf("Wrong! You did not set CCP4 environment for the shell.\n");
        else
            printf("Wrong! You did not set RCSBROOT environment for the shell.\n");
        printf("Please look at the README instruction in the package.\n");
        return;
    }else{  
        strcpy(text_path,"");
        strcpy(text_path,temp);
        
        if(!strcmp(CCP4,"yes")){
            strcat(text_path,CCP4_2_DATA);
        }else{
            strcat(text_path,EXTRACT_2_DATA);
        }
        strcat(adit_text,text_path);
        strcat(adit_inp, text_path);

        
        if(!strcmp_case(METHOD, "NMR")){
            strcat(adit_text,"nmr_data_template.txt");
            strcat(adit_inp, "nmr_log_script.inp");
        
            if((fadit_inp_in = fopen(adit_inp,"r"))==NULL){
                printf("WARNING! Can not open the file (%s).\n",adit_inp);
            }
            if((fadit_text_in = fopen(adit_text,"r"))==NULL){
                printf("Please use the file 'data_template.text' for pdb_extract.\n");
            }
        }else{
            strcat(adit_text,"data_template.text");
            strcat(adit_inp, "log_script.inp");
 
            if((fadit_inp_in = fopen(adit_inp,"r"))==NULL){
                printf("WARNING! Can not open the file 'log_script.inp'.\n");
            }
            if((fadit_text_in=fopen(adit_text,"r"))==NULL){
                printf("Please use the file 'data_template.text' for pdb_extract.\n");
            }
            
        }
        
    }
    
    if(fadit_text_in!=NULL) write_entity_poly_text(fadit_text_in,  fadit_text_out);
    
    if(fadit_inp_in != NULL) { 
        while(fgets(str, sizeof str, fadit_inp_in)!=NULL){
            fprintf(fadit_inp_out, "%s", str);
        }
        fclose(fadit_inp_out);
    }
   
   
    fprintf(fcgi, "a= \" %s\": b = \" %s\": c = \" %s\":  alpha = \" %s\":  beta = \" %s\":  gamma =\" %s\": space_group = \" %s\"\n",
            cell.a, cell.b, cell.c, cell.alpha, cell.beta, cell.gamma, cell.space_group); 

 
    num=read_pdb(pdbfile,AtomName, ResName, ChainID, ResSeq, xyz, Miscs,
                 ALT_LIST);

    delete_ligand(num, AtomName,ResName,ChainID, ResSeq,xyz, Miscs, &natom);
    if(natom<=0){
        printf("\n%s is neither Protein nor nucleic acids!\n",  pdbfile);
        fprintf(fout,"\n%s is neither Protein nor nucleic acids!\n",  pdbfile);
        return;
    }
    printf("PDB data file name: %s (%ld)\n",  pdbfile, num);
   
/* get the numbering information of each residue.
   seidx[i][j]; i = 1-num_residue  j=1,2
*/
    seidx = residue_idx(natom, ResSeq, Miscs, ChainID, ResName, &num_residue);
    bseq = cvector(0, num_residue);
    RY = lvector(0, num_residue);
    get_one_lett_seq(num_residue, seidx, AtomName, ResName, ChainID,
                     ResSeq, Miscs,xyz, bseq, RY);  
     
    chain_idx = lmatrix(0,500 , 0, 2);  /* # of chains max = 500 */    
    get_chain_idx(num_residue, seidx, ChainID, RY, &nchain, chain_idx);

    chain_res = cmatrix(0, nchain+1, 0, num_residue);
    chain_res_new = cmatrix(0, nchain+1, 0, num_residue+4);
    half_resd = 1+num_residue/2;
    broken_chain_idx = lmatrix(0,nchain+1, 0, half_resd);
    nbroken = lvector(0,nchain+1);
    matched = cvector(0, 2*nchain) ;


    for (i=1; i<=nchain; i++){ 
        k=chain_idx[i][1];
        chain_nam1 = ChainID[ seidx[k][1] ];
        k=chain_idx[i][2] - chain_idx[i][1] +1;
        broken_chain(i, chain_idx, seidx, xyz, RY, AtomName, broken_chain_idx,
                     nbroken);
        
        printf("%ld, Chain-ID  %c: Residue Number %4ld\n",i, chain_nam1, k);
        fprintf(fout,"\n%ld, Chain-ID  %c: Residue Number %4ld\n",i,chain_nam1,k);

        j=0;
        j1=0;
        strcpy(chain_res[i], "");
        for (k=chain_idx[i][1]; k<=chain_idx[i][2]; k++){
            chain_res[i][j]=bseq[k];
            chain_res_new[i][j1]=bseq[k];
            
            if(bseq[k]=='X' || bseq[k]=='?'){
                fprintf(fout,"(%s)",ResName[ seidx[k][1] ]);
            }else{
                fprintf(fout,"%c",bseq[k]);
            }
            for (k1=1; k1<=nbroken[i]; k1++){
                if(k==broken_chain_idx[i][k1]){
                    fprintf(fout,"????");
                    chain_res_new[i][++j1]='?';
                    chain_res_new[i][++j1]='?';
                    chain_res_new[i][++j1]='?';
                    chain_res_new[i][++j1]='?';
                }                
            }
            
            if((j1+1)%70==0) fprintf(fout,"\n");
            j++;
            j1++;
        } 
        chain_res[i][j]='\0';
        chain_res_new[i][j1]='\0';
        fprintf(fout,"\n");
        
        for (k=1; k<=nbroken[i]; k++){
            k1=broken_chain_idx[i][k];
            printf("Chain (%ld) %c: broken at residue %c: (sequence %ld)\n"
                   ,i,  chain_nam1, bseq[k1], ResSeq[seidx[k1][1]]);
            fprintf(fout,"Chain (%ld) %c: broken at residue %c: (sequence %ld)\n"
                   ,i,chain_nam1, bseq[k1], ResSeq[seidx[k1][1]]);
        }
        
    }


    npoly=0;
    ntot=0;
    strcpy(matched,"");
    strcpy(matched_comp,"");
    for (i=1; i<=nchain; i++){
        k1=chain_idx[i][1];
        chain_nam1 = ChainID[ seidx[k1][1] ];
        if(strchr(matched_comp, chain_nam1)) continue; /* matched with this chain*/
        if(RY[k1]== -2)continue; /*rid of ligand*/
        if(chain_idx[i][2]-chain_idx[i][1]<1)continue;/* too few (only one)*/
        if(chain_idx[i][2]-chain_idx[i][1]==nbroken[i])continue; /*all broken*/

        k=0;
        if(i!=nchain){
            for (j=i+1; j<=nchain; j++){
                k2=chain_idx[j][1];
                if(RY[k1] != RY[k2]) continue;
            
                chain_nam2 = ChainID[ seidx[k2][1] ];
            
                if(RY[k1] >=0 && RY[k2]>=0) { /* Nucleic acids */
                    fprintf(fout,"Nucleic Acids %c %c\n",chain_nam1, chain_nam2);
                    
                
//                    fasta_seq_align(i, j, chain_res, &score);
                    
                    seq_alignment(i, j, chain_res, &percent,
                                  &nmatch, seq1_match, seq2_match);

                    if(percent>98.0){
                        matched_comp[ntot]=chain_nam2;
                        ntot++;
                        matched[k]=chain_nam2;
                        k++;
                        matched[k]=',';
                        k++;
                        fprintf(fout,"%s RNA/DNA: Similarity(%c:%c:)= %6.2f; nmatch= %2ld rmsd = ?\n",
                                pdbfile, chain_nam1, chain_nam2, percent, nmatch);
                    
                    }

                }else if(RY[k1]== -1 && RY[k2]== -1) {  /*protein*/
                    seq_alignment(i, j, chain_res, &percent,
                                  &nmatch, seq1_match, seq2_match);
                    if(percent>50){ 

                        len1=strlen(chain_res[i])+2;
                        len2=strlen(chain_res[j])+2;

                        large = (len1 >= len2) ?  len1 : len2;
                       
                        struct_3d_align_new(i, j, large,chain_idx, AtomName, seidx,
                                            xyz, RY, nmatch, seq1_match, seq2_match,
                                            &rmsd_new);
                        if(len1<=len2){
                            struct_3d_align(i, j, large,chain_idx, AtomName, seidx,
                                            xyz, RY, &rmsd);
                        }else {                    
                           struct_3d_align(j, i, large,chain_idx, AtomName, seidx,
                                        xyz, RY, &rmsd);
                        }
                    }else{
                        rmsd=99;
                        rmsd_new=99;
                    }
                    if(percent>50.0){
                        matched[k]=chain_nam2;
                        matched_comp[ntot]=chain_nam2;
                        ntot++;
                        k++;
                        matched[k]=',';
                        k++;
                        fprintf(fout,"%s PROTEIN:  Similarity(%c:%c:)= %6.2f; rmsd= %6.2lf rmsd_new= %6.2lf nmatch=%2ld\n",
                                pdbfile,chain_nam1, chain_nam2, percent, rmsd, rmsd_new,nmatch);
                    }
                }
            }
        }
        
        matched[k]='\0';
        if(i==nchain) strcpy(matched,"");

        
        npoly++;
        fprintf(fadit_text_out,"<molecule_entity_id=\"%ld\" >\n",npoly);
        fprintf(mmcif," %ld ",npoly);
        fprintf(fcgi,"ENTITY_ID = \" %ld\" :",npoly);
        

        if(RY[k1] >=0){            
            fprintf(mmcif, " POLYMER ");            
            fprintf(fadit_text_out, "<molecule_entity_type=\"POLYMER\" >\n");   
            fprintf(fcgi, " ENTITY_TYPE=\"POLYMER\"  :");            
        }else if(RY[k1] ==-1){
            fprintf(mmcif, " polypeptide(L) ");            
            fprintf(fadit_text_out, "<molecule_entity_type=\"polypeptide(L)\" >\n");
            fprintf(fcgi, " ENTITY_TYPE=\"polypeptide(L)\" :");
        }
        fprintf(fadit_text_out, "<molecule_one_letter_sequence=\" \n");
        fprintf(fcgi, " ONE_LETTER_SEQUNCE=\" ");
        
        j1=0;

        fprintf(mmcif,"\n;");            
        for (k=chain_idx[i][1]; k<=chain_idx[i][2]; k++){
            if(bseq[k]=='X' || bseq[k]=='?'){
                fprintf(mmcif,"(%s)", ResName[ seidx[k][1] ]);            
                fprintf(fadit_text_out,"(%s)", ResName[ seidx[k][1] ]);
                fprintf(fcgi,"(%s)",ResName[ seidx[k][1] ]);
            }else{
                fprintf(mmcif, "%c", bseq[k]);            
                fprintf(fadit_text_out,"%c",bseq[k]);
                fprintf(fcgi,"%c",bseq[k]);
            }
            for (k1=1; k1<=nbroken[i]; k1++){
                if(k==broken_chain_idx[i][k1]){
                    fprintf(mmcif,"????");
                    fprintf(fadit_text_out,"????");
                    fprintf(fcgi,"????");
                }                
            }
            j1++;
            if((j1)%70==0) {
                fprintf(mmcif,"\n");
                fprintf(fadit_text_out,"\n");
            }
            
        }
        fprintf(mmcif,"\n;\n");
        fprintf(fadit_text_out,"\" >\n");
        fprintf(fcgi,"\" :");
        
        sprintf(tmp, "%c,%s",chain_nam1,matched);
        tmp[strlen(matched)+1]='\0';
        if(tmp[0]==' ')tmp[0]='?';
        
        fprintf(mmcif, " %s\n",tmp);
        fprintf(fadit_text_out, "< molecule_chain_id=\"%s\" >\n\n",tmp);
        fprintf(fcgi, " CHAIN_ID=\"%s\" ;",tmp);

    }

    fprintf(fadit_text_out,"<molecule_entity_id=\" \" >\n");
    fprintf(fcgi,"ENTITY_ID = \" \" :");
    
    fprintf(fadit_text_out, "<molecule_entity_type=\" \" >\n");            
    fprintf(fcgi, " ENTITY_TYPE=\" \"  :");            

    fprintf(fadit_text_out, "<molecule_one_letter_sequence=\"  \" >\n");
    fprintf(fcgi, " ONE_LETTER_SEQUNCE=\"  \" :");

    fprintf(fadit_text_out, "<molecule_chain_id=\" \" >\n\n");
    fprintf(fcgi, " CHAIN_ID=\"  \" ;");

    while(fadit_text_in && fgets(str, sizeof str, fadit_text_in)!=NULL){
        fprintf(fadit_text_out, "%s", str);
    }

/*    
    fprintf(fadit_text_out, "NEVER CHANGE ANYTHING BELOW!! OTHERWISE, THINGS WILL MESS UP!\n");
    fprintf(fadit_text_out, "THE FORMAT OF THE APPENDANT BELOW IS (%s)\n", format);
    fprintf(fadit_text_out, "==============================================================\n");
// note: if the above lines changed, the route (write_text_appendant ) has to change
    
    while(fgets(str, sizeof str, fp)!=NULL){
        fprintf(fadit_text_out, "%s", str);
    }
*/
    
    
    fclose(fout);
    fclose(fadit_text_out);
    fclose(mmcif);
    write_matth_coeff_to_data_template((char *) "data_template.text");
  
    free_cmatrix(AtomName, 0, num, 0, 4);
    free_cmatrix(ResName, 0, num, 0, 3);
    free_cvector(ChainID, 0, num);
    free_lvector(ResSeq, 0, num);
    free_dmatrix(xyz, 0, num, 0, 3);
    free_cmatrix(Miscs, 0, num, 0, NMISC);
    free_lmatrix(seidx, 0, num_residue, 0, 2);
    free_lmatrix(seq1_match , 0,40, 0,3);
    free_lmatrix(seq2_match , 0,40, 0,3);
    
    free_cvector(bseq, 0, num_residue);
    free_lvector(RY, 0, num_residue);
    free_lmatrix(chain_idx, 0,500 , 0, 2);
    free_cmatrix(chain_res , 0, nchain+1, 0, num_residue);
    free_cmatrix(chain_res_new , 0, nchain+1, 0, num_residue+4);

    free_lmatrix(broken_chain_idx , 0, nchain+1,0, half_resd);
    free_lvector(nbroken ,0,nchain+1);
    free_cvector(matched , 0, 2*nchain) ;
    delete_file(outfile);
    
    return;
    
    
}

void write_matth_coeff_to_data_template(char *inpfile)
{
    char  value[50000],seq[50000],  item[60], str[100];
    char seq_type[60], chain_id[512];
    char sg_name_H_M[20],sg_name_table[20];
    int  key=0,  n1=0, n2=0, n3=0, nsym;
    float mass=0, cell[7], cell_volume=0, matthew_coeff=0, solvent_const=0;
    FILE  *fp=NULL, *ftmp;

        
    if((fp=fopen(inpfile, "r"))==NULL){
        return;
    }
    
    ftmp = tmpfile();
    while(fgets(str, sizeof str, fp) !=NULL){
        fprintf(ftmp, "%s", str);
    }
    
    rewind(ftmp);
    while(1) {
        values_in_bracket(ftmp,  item, value,  &key);
        if(key>0)break;
        
        if(!strcmp_case(item, "space_group")){
            strcpy(sg_name_H_M, value);
            
        }else if(!strcmp_case(item, "unit_cell_a") ||
                 !strcmp_case(item, "a")){
            cell[0] = atof (value);
            
        }else if(!strcmp_case(item, "unit_cell_b")||
                 !strcmp_case(item, "b")){
            cell[1] = atof (value);
           
        }else if(!strcmp_case(item, "unit_cell_c")||
                 !strcmp_case(item, "c")){
            cell[2] = atof (value);
            
        }else if(!strcmp_case(item, "unit_cell_alpha")||
                 !strcmp_case(item, "alpha")){
            cell[3] = atof (value);
           
        }else if(!strcmp_case(item, "unit_cell_beta")||
                 !strcmp_case(item, "beta")){
            cell[4] = atof (value);
            
        }else if(!strcmp_case(item, "unit_cell_gamma")||
                 !strcmp_case(item, "gamma")){
            cell[5] = atof (value);
        }

        
        if(strstr_case(item, "entity_type")){
            strcpy(seq_type, value);
            n1=1;
        }else if(strstr_case(item, "one_letter")){
            strcpy(seq, value);
            n2=1;
        }else if(strstr_case(item, "chain_id")){
            strcpy(chain_id, value);
            n3=1;
        }
        if(n1>0 && n2>0 && n3>0 && strlen(seq)>0 && strlen(chain_id)>0){
            if(strstr_case(seq_type,"POLYMER")){
                mass = mass + residue_mass(seq, chain_id, (char *)"nucleic_acids");
                
            }else if(strstr_case(seq_type,"polypeptide")){
                mass = mass + residue_mass(seq, chain_id, (char *)"amino_acids");
            }
            n1=0; n2=0; n3=0;
        }
    }
    
    cell_volume = get_cell_volume(cell);
    nsym = get_symmetry_num(sg_name_H_M, sg_name_table);
    mass = nsym*mass;
    
    if(cell_volume >1 && mass>0){
        matthew_coeff = cell_volume/mass;
        solvent_const =(1.0-1.0/(matthew_coeff * 0.602*1.34))*100.0;
    }

//    printf("(%10s) %2d %7.2f %7.2f ::%10.1f %10.1f\n",sg_name_H_M, nsym,
//                   matthew_coeff,solvent_const, mass, cell_volume);

    if(matthew_coeff<10 && matthew_coeff>0.5 && solvent_const<95 && solvent_const>5){
        fp=fopen("data_template.text", "w");
        rewind(ftmp);
        while(fgets(str, sizeof str, ftmp) !=NULL){
            
            if(pattern5(str, "<", "crystals_solvent_content_", "=", ">" , "")){
                fprintf(fp, "<crystals_solvent_content_1 = \"%.1f \">\n",solvent_const);
            }else if (pattern5(str, "<", "crystals_matthews_coefficient_", "=", ">" , "")){
                fprintf(fp, "<crystals_matthews_coefficient_1 = \"%.1f \">\n",matthew_coeff);
            }else{
                fprintf(fp,"%s", str);
            }
        }
    }
    fclose(ftmp);
}

float get_cell_volume(float *cell)
{
    float  alpha, beta,gamma ,cell_volume;
    
    alpha = 3.14159 * cell[3]/180;
    beta  = 3.14159 * cell[4]/180;
    gamma = 3.14159 * cell[5]/180;
    
    cell_volume = cell[0] * cell[1] * cell[2] *
        sqrt(1.- cos(alpha)*cos(alpha) - cos(beta)*cos(beta) - cos(gamma)*cos(gamma)
             + 2.0 * cos(alpha) * cos(beta) * cos(gamma));
    return cell_volume;
}

float residue_mass(char *seq, char *chain_id, char *residue_type)
{
    int i, j, k, nres=0, n_chain=1, n_modify=0, len_id, len_seq,nb=0;
    float mass = 0, mass_tot=0;
    
    
    len_id = strlen(chain_id);
    len_seq = strlen(seq);
    if(len_id<=0 || len_seq<=0) return 0;

    for(i=0; i< len_id; i++){
        if(chain_id[i] == ',')n_chain++;
    }
    
    for(i=0; i< len_seq; i++){
        if(seq[i] == '?' ){
            nb++;
            continue;
        }
        
        if(seq[i] == ' ' || seq[i] == '\n')continue;
        if(seq[i] == '('){
            j=i;
            k=0;
            for(i=j; i< len_seq; i++){
                if(seq[i] == ' ' || seq[i] == '\n')continue;
                k++;
                if(seq[i] == ')'){
                    n_modify++;
                    if(k>5)printf("Warning! Characters (=%d) in () > 3\n", k);
                    
                    break;
                }
            }
            if(k>0) continue;      
        }
        nres++;
//        printf("%c", seq[i]);
        if(!strcmp_case(residue_type, "nucleic_acids")){
            mass = mass + single_residue_mass((char *)"nucleic_acids", seq[i]);
        }else if (!strcmp_case(residue_type, "amino_acids")){
            mass = mass + single_residue_mass((char *)"amino_acids", seq[i]);
        } 
    }
    
//    printf("\n modify= res=%d %d  %d  %.2f\n", n_chain, n_modify, nres, mass );

    
    
    if(!strcmp_case(residue_type, "nucleic_acids")){
        mass_tot = (mass + n_modify*300)* n_chain;
        return mass_tot ;
    } else if (!strcmp_case(residue_type, "amino_acids")){
        mass_tot = (mass + n_modify*110)* n_chain;
        return mass_tot ;
    } else{
        mass_tot = (mass + n_modify*110)* n_chain;
        return mass_tot ;
    }
    
    
}

float single_residue_mass(char *residue_type, char residue_name)
{
    int i;

    if(strstr_case(residue_type, "nucleic_acids")) {
        for (i = 0; i < NUM_MASS_TABLE_NUCLEIC_ACIDS; i++) {
            if (mass_table_nucleic_acids[i].Residue_1[0] == residue_name){
                return mass_table_nucleic_acids[i].Residue_Mass_H;
            }
        }
        return 290;
    }else if(strstr_case(residue_type, "amino_acids")) {
        for (i = 0; i < NUM_MASS_TABLE_AMINO_ACIDS; i++) {
            if (mass_table_amino_acids[i].Residue_1[0] == residue_name){
                return mass_table_amino_acids[i].Residue_Mass_H;
            }
        }
        return 110;
           
    }else{
        return 110;
    }
}

int get_symmetry_num(char *sg_name_inp, char *sg_name_out)
{
       int i, n=0;
       char sg_name_inp_tmp[20],sg_name_table_tmp[20];

       strcpy(sg_name_inp_tmp, sg_name_inp);
       rid_of_space(sg_name_inp_tmp);
           
       for (i = 0; i < NUM_SPACE_GROUP; i++) {
           strcpy(sg_name_table_tmp, space_group_lib[i].space_group_name);
           rid_of_space(sg_name_table_tmp);
           
           if (!strcmp_case(sg_name_inp_tmp, sg_name_table_tmp)){
               strcpy(sg_name_out, space_group_lib[i].space_group_name);
               n = space_group_lib[i].nsym;
               break;
           }
       }
       if(n==0){ //space group name wrong!
           strcpy(sg_name_out, sg_name_inp);
           n=1;
       }
       return n;
}


void write_entity_poly_text(FILE *fadit_in,  FILE *fadit_out)
{
    char str[256];

    if(!strcmp_case(METHOD, "NMR")){
        while(fgets(str, sizeof str, fadit_in)!=NULL){
            fprintf(fadit_out, "%s", str);
	    //            if(strstr_space(str, "cases of residues modeled as Ala or Gly due to missing side-chain density")){
            if(strstr_space(str, "Please manually add the sequence")){
                fprintf(fadit_out,"\n");
                break;
            } 
        }
    }else{ //Xray
        while(fgets(str, sizeof str, fadit_in)!=NULL){
            fprintf(fadit_out, "%s", str);
            if(strstr_space(str, "Enter Crystallographic Data")) {
                
                fprintf(fadit_out,"\n<space_group = \"%s\"> (use International Table conventions)\n",
                        cell.space_group);
                fprintf(fadit_out,"<space_group_number = \"? \">\n\n");
                fprintf(fadit_out,"<unit_cell_a     = \"%s\" >\n", cell.a);
                fprintf(fadit_out,"<unit_cell_b     = \"%s\" >\n", cell.b);
                fprintf(fadit_out,"<unit_cell_c     = \"%s\" >\n", cell.c);
                fprintf(fadit_out,"<unit_cell_alpha = \"%s\" >\n", cell.alpha);
                fprintf(fadit_out,"<unit_cell_beta  = \"%s\" >\n", cell.beta);
                fprintf(fadit_out,"<unit_cell_gamma = \"%s\" >\n", cell.gamma);

                break;
            }
        }
        
        while(fgets(str, sizeof str, fadit_in)!=NULL){
            fprintf(fadit_out, "%s", str);
            if(strstr_space(str, "cases of residues modeled as Ala or Gly due to missing side-chain density")){
                fprintf(fadit_out,"\n");
                break;
            } 
        }
        
    }
    
    
}



void read_pdb_head(char *pdbfile)
{
    char str[100];
    FILE *fp;
	
    if((fp = fopen(pdbfile, "r"))==NULL){
        open_file_error(pdbfile, "read_pdb_head");
    }
    
    strcpy(cell.space_group, "?");
    strcpy(cell.a, "?");
    strcpy(cell.b, "?");
    strcpy(cell.c, "?");
    strcpy(cell.alpha, "?");
    strcpy(cell.beta, "?");
    strcpy(cell.gamma, "?");
    while (fgets(str, sizeof str, fp) != NULL) {
        upperstr(str);
        if (!strncmp(str, "CRYST1", 6)){
            strncpy(cell.a,str+7, 9);
            cell.a[10]='\0';
            
            strncpy(cell.b,str+16, 9);
            cell.b[10]='\0';
            
            strncpy(cell.c,str+25, 9);
            cell.c[10]='\0';
            
            strncpy(cell.alpha,str+34, 7);
            cell.alpha[8]='\0';
            
            strncpy(cell.beta,str+41, 7);
            cell.beta[8]='\0';
            
            strncpy(cell.gamma,str+48, 7);
            cell.gamma[8]='\0';
            
            strncpy(cell.space_group,str+55, 10);
            cell.space_group[11] = '\0';
            break;
        }else if (!strncmp(str, "ATOM", 4)){
            break;
        }

    }
    fclose(fp);
 
}

void get_mmcif(char *inpfile)
/* parse the information from log_script.inp and use the pdb_extract program
   to get a mmcif file.
*/
{
    char lett[10000];  /* characters in <> */
    char value[10000]; /* value in " " of <> */
    char item[80];  /* characters for the first item in <> */
    
    char sf_dt[5],sf_prog[20], sf_cryst[5],sf_diff[5], sf_file[MAXL];
    char ref_sf_dt[5],ref_sf_prog[20],ref_sf_file[MAXL],log_sf_file[MAXL];
    char scale_p[20], scale_log[MAXL],scale_cif[MAXL];
    char index_p[20],index_log[MAXL],index_cif[MAXL];
    char mr_e[10],mr_p[20],mr_log[MAXL];
    char phase_e[10],phase_p[20],phase_log[MAXL],phase_pdb[MAXL],phase_cif[MAXL];
    char dm_d[20], dm_log[MAXL], dm_pdb[MAXL], dm_cif[MAXL];
    char ref_r[20], ref_log[MAXL], ref_pdb[MAXL], ref_cif[MAXL];
    
    char data_template[MAXL], data_template_tmp[MAXL];
    char tmp[2000], sf_output[MAXL], stat_output[MAXL];
    char command_str[2000], sf_command[2000],ref_sf_command[2000], stat_command[2000];
    char *temp,  pdb_extract[MAXL], pdb_extract_sf[MAXL];
    int  entity = 0;
    long key=0;
    FILE *finp, *fadit;

    
    if((finp=fopen(inpfile,"r"))==NULL){
        open_file_error(inpfile, "get_mmcif");
    }
    
    (!strcmp(CCP4,"yes"))? temp = getenv("CCP4"): temp = getenv("RCSBROOT");
    if(temp != NULL) {
        strcpy(pdb_extract, "");
        strcpy(pdb_extract_sf, "");
        strcat(pdb_extract, temp);
        strcat(pdb_extract, "/bin/pdb_extract ");
        strcat(pdb_extract_sf, temp);
        strcat(pdb_extract_sf, "/bin/pdb_extract_sf ");
    }else{
        if(!strcmp(CCP4,"yes"))
            printf("Wrong! You did not set CCP4 environment for the shell.\n");
        else
            printf("Wrong! You did not set RCSBROOT environment for the shell.\n");
        printf("Please look at the README instruction in the package.\n");
        return;
    }
    


    strcpy(ref_sf_command, "");
    strcpy(sf_command, "");
    strcpy(stat_command, "");
    
    strcpy(phase_cif,""); strcpy(phase_pdb,""); strcpy(phase_log,"");
    strcpy(dm_cif,""); strcpy(dm_pdb,""); strcpy(dm_log,"");
    strcpy(ref_cif,""); strcpy(ref_pdb,"");  strcpy(ref_log,"");
    strcpy(scale_cif,""); strcpy(scale_log,"");strcpy(log_sf_file,"");
    strcpy(index_cif,""); strcpy(index_log,"");
    
    
    strcpy(sf_output,"pdb_extract_sf.mmcif"); /*output default*/
    strcpy(stat_output,"pdb_extract.mmcif");  /*output default*/
    
    while(1) { /* parsing values  in < >*/
        element_in_bracket(finp,  item, value, lett, &key);
        if(key>0)break;
        if( strlen(value)<=0 || strlen(lett)<=0) continue;
        
        if(!strcmp_case(item,"scale_data_type") || /*SF for scale*/
           !strcmp_case(item,"data_type")){ 
            strcpy(sf_dt, value);
        }else if(!strcmp_case(item,"scale_program_name") ||
                 !strcmp_case(item,"program_name") ){
            strcpy(sf_prog, value);
        }else if(strstr_case(item,"crystal_number")){
            strcpy(sf_cryst, value);
        }else if(strstr_case(item,"diffract_number")){
            strcpy(sf_diff, value);
        }else if(strstr_case(item,"scale_data_file_name") ||
                 !strncmp_case(item,"data_file_name", 14) ){
            strcpy(sf_file, value);
        }else if(strstr_case(item,"scale_log_file_name")||
                 !strncmp_case(item,"log_file_name", 13) ){
            strcat(log_sf_file, value);
            strcat(log_sf_file, " ");

        }else if(strstr_case(item,"reflection_data_type")|| /*SF for refine*/
                 !strcmp_case(item,"refine_data_type")){ 
            strcpy(ref_sf_dt, value);
        }else if(strstr_case(item,"reflection_data_format")||
                 !strcmp_case(item,"refine_program_name")){
            strcpy(ref_sf_prog, value);
        }else if(strstr_case(item,"reflection_data_file_name")||
                 !strcmp_case(item,"refine_data_file_name")){
            strcpy(ref_sf_file, value);

            
        }else if(strstr_case(item,"data_indexing_software")){/*Data Indexing */
            strcpy(index_p, value);
        }else if(strstr_case(item,"data_indexing_LOG_file_name")){
            strcpy(index_log, value);
        }else if(strstr_case(item,"data_indexing_CIF_file_name")){
            strcpy(index_cif, value);

        }else if(strstr_case(item,"data_scaling_software")){/*Data Scaling */
            strcpy(scale_p, value);
        }else if(strstr_case(item,"data_scaling_LOG_file_name")){
            strcpy(scale_log, value);
        }else if(strstr_case(item,"data_scaling_CIF_file_name")){
            strcpy(scale_cif, value);


        }else if(!strncmp_case(item,"mr_software",10)){  /*MR*/
            strcpy(mr_e, "MR");
            strcpy(mr_p, value);
        }else if(!strncmp_case(item,"mr_log_file_log_",14)){
            strcat(mr_log, value);
            strcat(mr_log, " ");

            
        }else if(!strncmp_case(item,"phasing_method",12)){  /*phasing*/
            strcpy(phase_e, value);
        }else if(!strncmp_case(item,"phasing_software",12)){
            strcpy(phase_p, value);
        }else if(!strncmp_case(item,"phasing_log_file_cif_",21)){
            strcat(phase_cif, value);            
            strcat(phase_cif, " ");            
        }else if(!strncmp_case(item,"phasing_log_file_pdb_",21)){
            strcat(phase_pdb, value);
            strcat(phase_pdb, " ");
        }else if(!strncmp_case(item,"phasing_log_file_log_",21)){
            strcat(phase_log, value);
            strcat(phase_log, " ");
            
        }else if(!strncmp_case(item,"dm_software",10)){  /*dm*/
            strcpy(dm_d, value);
        }else if(!strncmp_case(item,"dm_log_file_cif_",16)){
            strcat(dm_cif, value);            
            strcat(dm_cif, " ");            
        }else if(!strncmp_case(item,"dm_log_file_pdb_",16)){
            strcat(dm_pdb, value);
            strcat(dm_pdb, " ");
        }else if(!strncmp_case(item,"dm_log_file_log_",16)){
            strcat(dm_log, value);
            strcat(dm_log, " ");
            
        }else if(!strncmp_case(item,"refine_software",12)){ /*refinement*/
            strcpy(ref_r, value);
        }else if(!strncmp_case(item,"refine_log_file_cif_",20)){
            strcat(ref_cif, value);            
            strcat(ref_cif, " ");            
        }else if(!strncmp_case(item,"refine_log_file_pdb_",20)){
            strcat(ref_pdb, value);
            strcat(ref_pdb, " ");
        }else if(!strncmp_case(item,"refine_log_file_log_",20) ){
            strcat(ref_log, value);
            strcat(ref_log, " ");
        
        }else if(!strcmp_case(item,"data_template_file") ){  /*Data Template File*/
            strcpy(data_template, value);

        }else if(!strncmp(item,"sf_output",8) ){  /*output file*/
            strcpy(sf_output, value);
        }else if(!strncmp(item,"statistics_output",16) ){
            strcpy(stat_output, value);
        }
        
        
        if(strlen(sf_prog)>0 && strlen(sf_dt)>0 && strlen(sf_file)>0
           && strlen(sf_cryst)>0 && strlen(sf_diff)>0){

            sprintf(tmp, " -c %s -w %s -idat %s ", sf_cryst, sf_diff, sf_file);
            strcat(sf_command, tmp);

            strcpy(sf_cryst, "");
            strcpy(sf_diff, "");
            strcpy(sf_file, "");
        }

        if(strlen(ref_sf_prog)>0 && strlen(ref_sf_dt)>0 && strlen(ref_sf_file)>0){
            sprintf(tmp, "-rt %s -rp %s -idat %s ",
                    ref_sf_dt, ref_sf_prog, ref_sf_file);
            strcat(ref_sf_command, tmp);
            
            strcpy(ref_sf_dt,   "");
            strcpy(ref_sf_prog, "");
            strcpy(ref_sf_file, "");
        }
        
    }
    
   
/* ------------------------Use pdb_extract_sf command-----------------------*/
    strcpy(command_str, ""); //initialize for the next
    
    
    if(strlen(sf_command)>0 && strlen(ref_sf_command)<=0 ){
        
        sprintf(command_str, "%s -dt %s -dp %s %s -o %s ",
                pdb_extract_sf, sf_dt, sf_prog, sf_command, sf_output);
        
        printf("\n============Structure factors used for phasing===============\n");
        printf("The script used:\n%s\n",command_str);
        system(command_str);
        
    }else if( strlen(sf_command)<=0 && strlen(ref_sf_command)>0 ){        
        sprintf(command_str, "%s  %s -o %s ",
                pdb_extract_sf, ref_sf_command, sf_output);
       
        printf("\n============Structure factosr used for refinement=============\n");
        printf("The script used:\n%s\n",command_str);
        system(command_str);
        
    }else if( strlen(sf_command)>0 && strlen(ref_sf_command)>0 ){
        sprintf(command_str, "%s %s -dt %s -dp %s %s -o %s ",pdb_extract_sf, 
                ref_sf_command, sf_dt, sf_prog, sf_command, sf_output);
        printf("\n=========Structure factors used for phasing & refinement=======\n");
        printf("The script used:\n%s\n",command_str);
        system(command_str);
    }else{
        printf("No structure factors will be extracted\n");
    }

/* --------------------------Use pdb_extract command------------------------*/

    strcpy(stat_command, "");
    strcat(stat_command, pdb_extract);
    
// extract scale statistics for indexing
    if (strlen(index_p)>0 &&
        (strlen(index_log)>0 || strlen(index_cif) ) ){

        strcat(stat_command, " -i ");
        strcat(stat_command, index_p);

        if(strlen(index_log)>0){
            strcat(stat_command, " -ilog ");
            strcat(stat_command, index_log);
        }
        
        if(strlen(index_cif)>0){
            strcat(stat_command, " -icif ");
            strcat(stat_command, index_cif);
        }
    }

    
// extract scale statistics for scaling
    if (strlen(scale_p)>0 &&
        (strlen(scale_log)>0 || strlen(scale_cif) ) ){

        strcat(stat_command, " -s ");
        strcat(stat_command, scale_p);

        if(strlen(scale_log)>0){
            strcat(stat_command, " -ilog ");
            strcat(stat_command, scale_log);
        }
        
        if(strlen(scale_cif)>0){
            strcat(stat_command, " -icif ");
            strcat(stat_command, scale_cif);
        }
    }

// extract scale statistics for phasing
    if (strlen(log_sf_file)>0 && strlen(sf_prog)>0){ 
        sprintf(tmp, " -sp  %s  -ilog %s ", sf_prog, log_sf_file);
        strcat(stat_command, tmp);
    }

/* extract MR statistics*/
    if (strlen(mr_p)>0 && strlen(mr_log)>0 ){  
        sprintf(tmp, " -e  MR  -m   %s  -ilog  %s ", mr_p, mr_log);
        strcat(stat_command, tmp);
    } 
 
/* extract phasing statis*/
    if(strlen(phase_e)>0 && strlen(phase_p)>0 && 
       (strlen(phase_cif)>0 || strlen(phase_pdb)>0 || strlen(phase_log)>0)){
        sprintf(tmp, " -e  %s  -p   %s  ", phase_e , phase_p);
        strcat(stat_command, tmp);
          
        if(strlen(phase_cif)>0){
            strcat(stat_command, " -icif ");
            strcat(stat_command, phase_cif);
        }
        
        
        if(strlen(phase_log)>0){
            strcat(stat_command, " -ilog ");
            strcat(stat_command, phase_log);
        }
        
        if(strlen(phase_pdb)>0){
            strcat(stat_command, " -ipdb ");
            strcat(stat_command, phase_pdb);
        }
        
    }
   
/* extract dm statis*/
    if(strlen(dm_d)>0  &&  
       (strlen(dm_cif)>0 || strlen(dm_pdb)>0 || strlen(dm_log)>0)){
        strcat(stat_command, " -d ");
        strcat(stat_command, dm_d);

        if(strlen(dm_cif)>0){
            strcat(stat_command, " -icif ");
            strcat(stat_command, dm_cif);
        }
        
        if(strlen(dm_pdb)>0){
            strcat(stat_command, " -ipdb ");
            strcat(stat_command, dm_pdb);
        }
        
        if(strlen(dm_log)>0){
            strcat(stat_command, " -ilog ");
            strcat(stat_command, dm_log);
        }
    }
    
/* extract refine statis*/        
    if(strlen(ref_r)>0  &&  
       (strlen(ref_cif)>0 || strlen(ref_pdb)>0 ||  strlen(ref_log)>0)){

        
        strcat(stat_command, " -r ");
        strcat(stat_command, ref_r);

        if(strlen(ref_cif)>0){
            strcat(stat_command, " -icif ");
            strcat(stat_command, ref_cif);
        }
        
        if(strlen(ref_pdb)>0){
            strcat(stat_command, " -ipdb ");
            strcat(stat_command, ref_pdb);
        }
        
        if(strlen(ref_log)>0){
            strcat(stat_command, " -ilog ");
            strcat(stat_command, ref_log);
        }   
    }

/* extract data template file*/        
    if(strlen(data_template)>0){
        strcat(stat_command, " -ient ");
        strcat(stat_command, data_template);
    }else { //test if the log_script.inp has entity_poly information
        
        rewind(finp);
        while(1) { 
            element_in_bracket(finp,  item, value, lett, &key);
            if(key>0)break;
            if( strlen(value)<=0 || strlen(lett)<=0) continue;
            rid_of_front_end_space(value);
            if(!strcmp_case(item,"ONE_LETTER_SEQUNCE") && strlen(value)>0){  
                entity=1;
                break;
            }
        }

        if(entity==1){
            sprintf(data_template_tmp, "adit_text_%d", getpid());
            fadit =  fopen(data_template_tmp, "w"); 
            parse_adit_text_to_mmcif(fadit, inpfile); 
            fclose (fadit);
            strcat(stat_command, " -ient ");
            strcat(stat_command, data_template_tmp);
        }
       
    }
    
     
    strcat(stat_command, " -o ");
    strcat(stat_command, stat_output);

    printf("\n============Structure Solution Statistics===============\n");
    printf("The script used:\n%s\n", stat_command);
    system(stat_command); /* phasing + dm + refinement stat */
    if( strlen(data_template_tmp)>1) delete_file(data_template_tmp);
    fclose(finp);
    printf("\nOutput file names:  %s & %s\n\n", stat_output, sf_output);
}



void write_entity_poly(FILE *fp, char *ent_id,char *ent_type,
                       char *ent_code,char *ent_chain)
{
    long i=0, len=0;

    len=strlen(ent_code);
    fprintf(fp,"%s  ",ent_id);
    fprintf(fp,"%s \n",ent_type);
    fprintf(fp,";");
    for(i=0;i<len;i++){
        fprintf(fp,"%c",ent_code[i]);
        if((i+1)%70==0) fprintf(fp,"\n");
    }
    fprintf(fp,"\n;");
    fprintf(fp,"\n%s\n",ent_chain);
}

        
void element_in_bracket(FILE *inp,char *item, char *value, char *lett,long *key)
/* item: the identifier.  size: total length of lett.
   lett: the total characters in <>.  key: =0 not the end, =1 end of file
*/
{
    int i=0,j=0,n=1, m=0,k=0;
    int c;
    long len;
    char *tmp;
    
    *key=0;
/* skip every letter until < is met */
    for(j=0; (c=getc(inp)) != '<' ; j++){
        if(c==EOF)  {
            *key=1;
            return;
        }
    }
/*take off the space between < and the first letter */    
    for(j=0; (c=getc(inp)) == ' ' ; j++){
        if(c==EOF)  {
            *key=1;
            return;
        }
    }
    lett[0]=c; 
    for(j=0; (c=getc(inp)) != '>' ; j++){
        if(c==EOF)  {
            *key=1;
            return;
        }
        if(c!='\n'){   /* omit \n, put all the chars in one line.*/                
            lett[n]=c;
            n=n+1;
        }
    }
    lett[n] = '\0';
    if(pattern5(lett, "=", "\"", "\"", "", "")<=0)return;
//    if(!strchr(lett, '=') && !strchr(lett, '\"'))return;
    
    len=n;
    tmp=cvector(0,len);

    i=strcspn(lett,"="); /* get the first item */

    m=0;
    for (k=0;k<i; ++k){
        if(isspace(lett[k])) continue;
        
        item[m]=tolower(lett[k]);
        m++;
    }
    item[m]='\0';
        

    strcpy(tmp, strchr(lett,'=')+1);  /*get string after the first word */
    j=0;
    n=strlen(tmp);
    for (i=0;i<n; ++i){
        if(!isspace(tmp[i]) && (tmp[i]!='\"')){
            value[j]=tmp[i];
            j++;
        }
    }    
    value[j]='\0';
    free_cvector(tmp, 0, len);
    
}