/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
**                                                                                 **
**                                                                                 **
**                                                       **                        **
**                                                       **                        **
**                       *********                       **                        **
**        **********    **********    ***        ***     **                        **
**        **********    **            ****      ****     **      **********        **
**        **            **            **  **  **  **     **      **********        **
**        **            ********      **  **  **  **     **      **                **
**        ******         ********     **    **    **     **      **                **
**        ******               ***    **    **    **     **************            **
**        **     ***     **    **********         **     **************            **
**        **     ***    *****************        ********        **                **
**        ***********   *********     **         *********       **                **
**        ***********    **    **     **         ***    **       **********        **
**               **  **  **    ******            **     **       **********        **
**               **  **  **    ******            ********                          **
**               **    * **    **                *********                         **
**               **    * **    **                **      **                        **
**               **     ***    **********        **      **                        **
**               **     ***    **********        **********                        **
**                                               *********                         **
**                                                                                 **
**                                                                                 **
**                          E N S E M B L E     T O O L S                          **
**                                                                                 **
**                                   Version 2.1                                   **
**                                                                                 **
**                                                                                 **
**                   Principal Investigator: Julie D. Forman-Kay                   **
**                                                                                 **
**                           Author:  Mickaël Krzeminski                           **
**                                                                                 **
**                               Date: November 2012                               **
**                                                                                 **
**                                                                                 **
*************************************************************************************
**                                                                                 **
**                                                                                 **
**   Copyright (C) The Hospital for Sick Children, 2001                            **
**                                                                                 **
**   Distribution of substantively modified versions of this module is prohibited  **
**  without the explicit permission of the copyright holder.                       **
**                                                                                 **
**   Any use of this work or derivative works in whole or in part for any          **
**  commercial purpose or for monetary gain is prohibited.                         **
**                                                                                 **
**                                                                                 **
**                                   NO WARRANTY                                   **
**   This software package is provided 'as is' without warranty of any kind,       **
**  expressed or implied.                                                          **
**                                                                                 **
**                                                                                 **
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#ifndef __PDB_HEADER__
#define __PDB_HEADER__


void _cleanPDB(PDBListData *ENSData)
{
    counter count;

    if (ENSData == NULL) return;

    for (count=0; count<(*ENSData).nbConfs; count++)
    {
        free((*ENSData).atomCoords[count]);
        (*ENSData).atomCoords[count] = NULL;
    }
    free((*ENSData).atomCoords); (*ENSData).atomCoords = NULL;

    for (count=0; count<(*ENSData).nbAtoms; count++)
    {
        free((*ENSData).atomHeaders[count]);
        (*ENSData).atomHeaders[count] = NULL;
    }
    free((*ENSData).atomHeaders); (*ENSData).atomHeaders = NULL;
}


char Resicode(char *residueName)
{
    if (!strcmp(residueName, "ALA")) return 'A';
    if (!strcmp(residueName, "CYS")) return 'C';
    if (!strcmp(residueName, "ASP")) return 'D';
    if (!strcmp(residueName, "GLU")) return 'E';
    if (!strcmp(residueName, "PHE")) return 'F';
    if (!strcmp(residueName, "GLY")) return 'G';
    if (!strcmp(residueName, "HIS")) return 'H';
    if (!strcmp(residueName, "ILE")) return 'I';
    if (!strcmp(residueName, "LYS")) return 'K';
    if (!strcmp(residueName, "LEU")) return 'L';
    if (!strcmp(residueName, "MET")) return 'M';
    if (!strcmp(residueName, "ASN")) return 'N';
    if (!strcmp(residueName, "PRO")) return 'P';
    if (!strcmp(residueName, "GLN")) return 'Q';
    if (!strcmp(residueName, "ARG")) return 'R';
    if (!strcmp(residueName, "SER")) return 'S';
    if (!strcmp(residueName, "THR")) return 'T';
    if (!strcmp(residueName, "VAL")) return 'V';
    if (!strcmp(residueName, "TRP")) return 'W';
    if (!strcmp(residueName, "TYR")) return 'Y';

    return 'X';
}


char isPDB (char *PDBFileName)
{
    /*
       This function is a fast checker of PDB files. Do not use with certainty as it
       verifies only the first 20 lines.
    */
    FILE *PDBFile;
    counter nbPDBLines = 0;
    char headerType[PDB_HEADER_LEN+1];        // To detect ATOM or HETATM
    char pdbFileLine[PDB_LINE_LEN];

    if ((PDBFile = fopen(PDBFileName, "rt")) == NULL) return '\0';

    headerType[PDB_HEADER_LEN] = '\0';
    while (fgets(pdbFileLine , MAX_PATHWAY_LEN, PDBFile) != NULL)
    {
        if (strlen(pdbFileLine ) > 5)
        {
            strncpy(headerType, pdbFileLine, 6);
            if (!(strcmp(headerType, "ATOM  ") && strcmp(headerType, "HETATM") && \
                  strcmp(headerType, "MODEL ") && strcmp(headerType, "ANISOU") && \
                  strcmp(headerType, "HEADER") && strcmp(headerType, "OBSLTE") && \
                  strcmp(headerType, "TITLE ") && strcmp(headerType, "SPLIT ") && \
                  strcmp(headerType, "CAVEAT") && strcmp(headerType, "COMPND") && \
                  strcmp(headerType, "SOURCE") && strcmp(headerType, "KEYWDS") && \
                  strcmp(headerType, "EXPDTA") && strcmp(headerType, "NUMMDL") && \
                  strcmp(headerType, "MDLTYP") && strcmp(headerType, "AUTHOR") && \
                  strcmp(headerType, "REVDAT") && strcmp(headerType, "SPRSDE") && \
                  strcmp(headerType, "JRNL  ") && strcmp(headerType, "REMARK") && \
                  strcmp(headerType, "DBREF ") && strcmp(headerType, "DBREF1") && \
                  strcmp(headerType, "DBREF2") && strcmp(headerType, "SEQADV") && \
                  strcmp(headerType, "SEQRES") && strcmp(headerType, "MODRES") && \
                  strcmp(headerType, "HET   ") && strcmp(headerType, "HETNAM") && \
                  strcmp(headerType, "HETSYN") && strcmp(headerType, "FORMUL") && \
                  strcmp(headerType, "HELIX ") && strcmp(headerType, "SHEET ") && \
                  strcmp(headerType, "SSBOND") && strcmp(headerType, "LINK  ") && \
                  strcmp(headerType, "CISPEP") && strcmp(headerType, "SITE  ") && \
                  strcmp(headerType, "CRYST1") && strcmp(headerType, "ORIGXn") && \
                  strcmp(headerType, "SCALEn") && strcmp(headerType, "MTRIXn") && \
                  strcmp(headerType, "CONECT") && strcmp(headerType, "MASTER")))
            {
                ++nbPDBLines;
                if (nbPDBLines == 20)
                {
                    break;
                }
            }
            else
            {
                fclose(PDBFile);
                return 0;
            }
        }
    }

    fclose(PDBFile);
    return 1;
}


int getMasses(float **atomMasses, PDBListData *ENSData)
{
    /*
       Returns an array that contains the masses of atoms from the headers of a PDB file.
    */
    counter atomNb, atomTypePos;
    char atomType, atomTypeString[2] = " \0";

    if (((*atomMasses) = malloc((*ENSData).nbAtoms * FLOAT_SIZE)) == NULL) fatal_error("MEMORY", "instance atomMass variable");

    for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++)
    {
        // Retrieving the first letter of the atom name
        atomType = ' ';
        for (atomTypePos=0; atomTypePos<4; atomTypePos++)
        {
            if ((*ENSData).atomHeaders[atomNb][PDB_ATOM_NAME_INDEX+atomTypePos] != ' ')
            {
                atomType = (*ENSData).atomHeaders[atomNb][PDB_ATOM_NAME_INDEX+atomTypePos];
                break;
            }
        }

        // Attributing the atom type masses
        switch (atomType)
        {
            case 'C':
                (*atomMasses)[atomNb] = MASS_C;
                break;
            case 'O':
                (*atomMasses)[atomNb] = MASS_O;
                break;
            case 'N':
                (*atomMasses)[atomNb] = MASS_N;
                break;
            case 'S':
                (*atomMasses)[atomNb] = MASS_S;
                break;
            case 'P':
                (*atomMasses)[atomNb] = MASS_P;
                break;
            case 'H':
            case '1':
            case '2':
            case '3':
                (*atomMasses)[atomNb] = MASS_H;
                break;
            default:
                (*atomMasses)[atomNb] = -1.;
        }

        if ((*atomMasses)[atomNb] + 1. < fEpsilon)
        {
            atomTypeString[0] = atomType;
            error("Unknown atom type: %s", atomTypeString);
            return 1;
        }
        else if (atomType == ' ')
        {
            error("No atom type specified!", NULL);
            return 1;
        }
    }

    return 0;
}


int getCoords(PDBListData *ENSData, char *ensFileName, confList *confNbsList)
{
    /*
      This function generates a variable of coordinates of all
      conformers found in the ensFileName file. This latter can
      be a CSP file, a list file, or a number list file.
      If the confList variable is different from NULL, only structure
      numbers from this list are taken into account.
    */

    counter confNb, atomNb;
    FILE *ENSFile;
    unsigned short PDBError = 0;
    unsigned long fileSize;
    char fileType[4] = {"\0"};
    char atomHeader[28] = {"\0"};

    // If the number of conformers is non-null, no need to record them again
    if ((*ENSData).nbConfs) return 0;

    // Opening the file
    if ((ENSFile = fopen(ensFileName, "rt")) == NULL)
    {
        error("NO_FILE", ensFileName);
        return 1;
    }

    // Determining the size of the file
    fseek(ENSFile, 0, SEEK_END);
    fileSize = ftell(ENSFile);                // !! Not more than 4 Gb !!
    if (!fileSize)
    {
        fclose(ENSFile);
        error("EMPTY", ensFileName);
        return 1;
    }
    fseek(ENSFile, 0, SEEK_SET);

    // Reading in the file type
    fread(&fileType, CHAR_SIZE, 3, ENSFile);
    fclose(ENSFile);

    if (!strcmp(fileType, "Csp"))
    {
        // Declaration of variables
        int32_t atomCoordInt;
        unsigned int coordStart, confBlockSize;
        unsigned long expFileSize;            // To check whether the CSP file is relevant

        // Reopening the file in binary mode
        ENSFile = fopen(ensFileName, "rb");
        fseek(ENSFile, 3 * CHAR_SIZE, SEEK_SET);

        // Reading the number of conformers and atoms
        fread(&(*ENSData).nbConfs, INT32_SIZE, 1, ENSFile);
        fread(&(*ENSData).nbAtoms, INT32_SIZE, 1, ENSFile);
        
        coordStart = (3 + 27 * (*ENSData).nbAtoms) * CHAR_SIZE + 2 * INT32_SIZE;
        confBlockSize = 3 * (*ENSData).nbAtoms * INT32_SIZE;

        // Comparing the effective and expected sizes of the file
        expFileSize = (unsigned long)(coordStart + (*ENSData).nbConfs * confBlockSize);
        if (fileSize != expFileSize)
        {
            fclose(ENSFile);
            printf ("%s is a corrupted CSP file (Expected size = %lu; Real size = %lu)!\n", ensFileName, expFileSize, fileSize);
            return 1;
        }
        fseek(ENSFile, 3 * CHAR_SIZE + 2 * INT32_SIZE, SEEK_SET);


        // Recording the atom atomHeaders and the residue numbers
        if (((*ENSData).atomHeaders = malloc((*ENSData).nbAtoms * sizeof(char *))) == NULL)
        {
            error("MEMORY", "instance (*ENSData).atomHeaders variable");
            return 1;
        }

        for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++)
        {
            if (((*ENSData).atomHeaders[atomNb] = malloc(28 * CHAR_SIZE)) == NULL)
            {
                error("MEMORY", "instance (*ENSData).atomHeaders[] variable");
                return 1;
            }
            fread(&atomHeader, CHAR_SIZE, 27, ENSFile);
            strcpy((*ENSData).atomHeaders[atomNb], atomHeader);
        }

        // Recording the coordinates
        if (confNbsList==NULL || (*confNbsList).confNb==NULL)
        {
            if (((*ENSData).atomCoords = malloc((*ENSData).nbConfs * sizeof(coordinates *)))==NULL)
            {
                error("MEMORY", "instance atomCoords variable");
                return 1;
            }

            for (confNb=0; confNb<(*ENSData).nbConfs; confNb++)
            {
                if (((*ENSData).atomCoords[confNb] = malloc((*ENSData).nbAtoms * COORD_SIZE))==NULL)
                {
                    error("MEMORY", "instance atomCoords[] variable");
                    return 1;
                }

                for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++)
                {
                    fread(&atomCoordInt, INT32_SIZE, 1, ENSFile);
                    (*ENSData).atomCoords[confNb][atomNb].x = (float)atomCoordInt / COORD_FACT;
                    fread(&atomCoordInt, INT32_SIZE, 1, ENSFile);
                    (*ENSData).atomCoords[confNb][atomNb].y = (float)atomCoordInt / COORD_FACT;
                    fread(&atomCoordInt, INT32_SIZE, 1, ENSFile);
                    (*ENSData).atomCoords[confNb][atomNb].z = (float)atomCoordInt / COORD_FACT;
                }
            }
        }
        else
        {
            // Checking for consistency of specified numbers
            // As confNbList is supposed to be sorted, only the last element is checked
            if ((*confNbsList).confNb[(*confNbsList).nb-1] > (*ENSData).nbConfs)
            {
                fclose(ENSFile);
                error("The specified file does not contain enough conformers!\n", NULL);
                return 1;
            }

            (*ENSData).nbConfs = (*confNbsList).nb;
            if (((*ENSData).atomCoords = malloc((*ENSData).nbConfs * sizeof(coordinates *)))==NULL)
            {
                error("MEMORY", "instance atomCoords variable");
                return 1;
            }

            for (confNb=0; confNb<(*ENSData).nbConfs; confNb++)
            {
                if (((*ENSData).atomCoords[confNb] = malloc((*ENSData).nbAtoms * COORD_SIZE))==NULL)
                {
                    error("MEMORY", "instance atomCoords[] variable");
                    return 1;
                }

                fseek(ENSFile, coordStart+((*confNbsList).confNb[confNb]-1)*confBlockSize, SEEK_SET);
                for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++)
                {
                    fread(&atomCoordInt, INT32_SIZE, 1, ENSFile);
                    (*ENSData).atomCoords[confNb][atomNb].x = (float)atomCoordInt / COORD_FACT;
                    fread(&atomCoordInt, INT32_SIZE, 1, ENSFile);
                    (*ENSData).atomCoords[confNb][atomNb].y = (float)atomCoordInt / COORD_FACT;
                    fread(&atomCoordInt, INT32_SIZE, 1, ENSFile);
                    (*ENSData).atomCoords[confNb][atomNb].z = (float)atomCoordInt / COORD_FACT;
                }
            }
        }

        fclose(ENSFile);

        return 0;
    }
    else if (!strcmp(fileType, "Ens"))
    {
        // Declaration of variables
        counter atomSeek;

        char trashChar[10], keyName[9];
        char trajCheck[14] = {"\0"};
        char atomHeader[28] = {"\0"};
        char residueNb[5] = {"\0"}, atomName[5] = {"\0"};
        char *ENSDirName = strdup(ensFileName);
        char *paramFileName = NULL;
        char **residueNbs = NULL, **atomNames = NULL, trajectoryName[MAX_PATHWAY_LEN];
        char **atomHeaders = NULL;

        unsigned short nbCsps = 0, cspNb, inTrajNb;
        unsigned short keyLen;
        unsigned short trajLen = 13;            // To detect the "trajectoriesh" from the parameter file
        short nbModules;
        short trashShort;
        short trajectoryNameLen;
        short found;

        unsigned int coordStart=0, confBlockSize = 0;    // The position before coordinates in the CSP file, and the size of all coordinated for one conformer
        unsigned int **atomHeaderPositions;        // To compare the positions of headers regarding a model

        uint32_t nbConfs, nbAtoms, nbRecConfs = 0;    // The number of conformers in each CSP file
        uint32_t *confToGet;                // The conformers to get read from the List file
        int32_t atomCoordInt;

        unsigned long expFileSize, paramFileSize;    // To check whether the CSP file (size) is relevant
        unsigned long fileSeek;

        double effNrg;

        CSPFileInfo *trajectories = NULL;

        FILE *paramFile, *CSPFile = NULL;        // The parameter file that contains the name of all CSP files

        // Reopening the file in binary mode
        ENSFile = fopen(ensFileName, "rb");
        fseek(ENSFile, 3 * CHAR_SIZE, SEEK_SET);

        // Size of a double variable
        fread(&DOUBLE_SIZE, SHORT_SIZE, 1, ENSFile);
        if (DOUBLE_SIZE != sizeof(double))
        {
            fclose(ENSFile);
            printf ("%s is an Ensemble file type generated on a different machine type!\n", ensFileName);
            return 1;
        }

        // Number of conformers
        fread(&(*ENSData).nbConfs, INT32_SIZE, 1, ENSFile);
        fseek(ENSFile, (*ENSData).nbConfs*INT32_SIZE, SEEK_CUR);    // Skipping the conformer numbers

        // Effective and ENSEMBLE energies
        fread(&effNrg, DOUBLE_SIZE, 1, ENSFile);

        // Number of modules
        fread(&nbModules, SHORT_SIZE, 1, ENSFile);

        // Checking the consistency of the file size (this must correspond to the expected file size)
        expFileSize = 3 * CHAR_SIZE + (nbModules + 2) * SHORT_SIZE + ((*ENSData).nbConfs + 1) * INT32_SIZE + (2 * nbModules + 1) * DOUBLE_SIZE;
        if (fileSize != expFileSize)
        {
            fclose(ENSFile);
            printf ("%s is a corrupted Ensemble file type (Expected size = %lu; Real size = %lu)!\n", ensFileName, expFileSize, fileSize);
            return 1;
        }

        // Reading the conformer numbers to retrieve
        fseek(ENSFile, 3*CHAR_SIZE+SHORT_SIZE+INT32_SIZE, SEEK_SET);
        if ((confToGet = malloc((*ENSData).nbConfs * INT32_SIZE)) == NULL)
        {
            error("MEMORY", "instance confToGet variable");
            return 1;
        }
        for (confNb=0; confNb<(*ENSData).nbConfs; confNb++) fread(&confToGet[confNb], INT32_SIZE, 1, ENSFile);

        fclose(ENSFile);

        // Memory allocation for recording the coordinates
        if (((*ENSData).atomCoords = malloc((*ENSData).nbConfs * sizeof(coordinates *)))==NULL)
        {
            error("MEMORY", "instance atomCoords variable");
            return 1;
        }

        // Finding the parameter file, which contains the names of the CSP files
        dirname(ENSDirName);
        if (!(strcmp(ENSDirName, ensFileName)))        // If the name is directly given, it returns the name, which is not what is wanted
        {
            ENSDirName = realloc(ENSDirName, 2*CHAR_SIZE);
            ENSDirName[0] = '.';
            ENSDirName[1] = '\0';
        }

        paramFileName = malloc((strlen(ENSDirName) + 22) * CHAR_SIZE);
        strcpy(paramFileName, ENSDirName);
        strcpy(&paramFileName[strlen(ENSDirName)], "/../Save/ensemble.dat");    // strcat would define a new variable, not strcpy

        if ((paramFile = fopen(paramFileName, "rb")) == NULL)
        {
            error("NO_FILE", paramFileName);
            return 1;
        }
        fseek(paramFile, 0, SEEK_END);
        paramFileSize = ftell(paramFile);                // !! Not more than 4 Gb !!
        fseek(paramFile, 0, SEEK_SET);

        // Retrieving the name of the CSP files from the parameter file
        for (fileSeek=0; fileSeek<=paramFileSize-trajLen; fileSeek++)
        {
            fseek(paramFile, fileSeek, SEEK_SET);
            fread(&trajCheck, CHAR_SIZE, trajLen, paramFile);

            if (!strcmp(trajCheck, "trajectoriesh"))
            {
                fread(&nbCsps, SHORT_SIZE, 1, paramFile);            // Number of CSP files
                if ((trajectories = malloc(nbCsps * sizeof(CSPFileInfo))) == NULL)
                {
                    error("MEMORY", "instance trajectories variable");
                    return 1;
                }
                break;
            }
        }

        if (!nbCsps)
        {
            error("%s is a corrupted user-defined parameter file!\n", ensFileName);
            return 1;
        }

        if ((atomHeaderPositions = malloc(nbCsps * sizeof(unsigned int *))) == NULL)
        {
            error("MEMORY", "instance atomHeaderPositions variable");
            return 1;
        }

        for (cspNb=0; cspNb<nbCsps; cspNb++)
        {
            fread(&trajectoryNameLen, SHORT_SIZE, 1, paramFile);            // Trajectory name length
            fread(&trajectoryName, CHAR_SIZE, trajectoryNameLen, paramFile);    // Trajectory name

            fread(&trashChar, CHAR_SIZE, 1, paramFile);                // 'h' because each trajectory is a hash in PERL
            fread(&trashShort, SHORT_SIZE, 1, paramFile);                // The number of elements in the hash (5: origin, tmp_name, first, last, type)

            for (inTrajNb=0; inTrajNb<5; inTrajNb++)
            {
                fread(&keyLen, SHORT_SIZE, 1, paramFile);            // Length of the key of the hash
                fread(keyName, CHAR_SIZE, keyLen+1, paramFile);            // Name of the key
                keyName[keyLen] = '\0';

                if (!strcmp(keyName, "origin"))
                {
                    fread(&keyLen, SHORT_SIZE, 1, paramFile);                        // Length of the corresponding value
                    if ((trajectories[cspNb].fileName = malloc((keyLen+1) * CHAR_SIZE)) == NULL) fatal_error("MEMORY", "instance trajectories[].fileName");
                    fread(trajectories[cspNb].fileName, CHAR_SIZE, keyLen, paramFile);    // Recording the corresponding value
                    trajectories[cspNb].fileName[keyLen] = '\0';

                    // Opening and checking out the consistency of the CSP (origin) file
                    if ((CSPFile = checkCSPFile(trajectories[cspNb].fileName, &nbConfs, &nbAtoms)) == NULL) return 1;

                    if (!cspNb)        // Reading the first CSP file
                    {
                        (*ENSData).nbAtoms = nbAtoms;

                        coordStart    = (3 + 27 * (*ENSData).nbAtoms) * CHAR_SIZE + 2 * INT32_SIZE;
                        confBlockSize = 3 * (*ENSData).nbAtoms * INT32_SIZE;

                        if ((atomHeaderPositions[0] = malloc((*ENSData).nbAtoms * INT_SIZE)) == NULL)
                        {
                            error("MEMORY", "instance atomHeaderPositions[0] variable");
                            return 1;
                        }

                        if ((residueNbs = malloc((*ENSData).nbAtoms * sizeof(char *))) == NULL)
                        {
                            error("MEMORY", "instance residueNbs variable");
                            return 1;
                        }

                        if ((atomNames  = malloc((*ENSData).nbAtoms * sizeof(char *))) == NULL)
                        {
                            error("MEMORY", "instance atomNames variable");
                            return 1;
                        }

                        // Recording the atom atomHeaders and the residue numbers
                        fseek(CSPFile, (3 * CHAR_SIZE + 2 * INT32_SIZE), SEEK_SET);        // Positionning the file handler right before the headers
                        if ((atomHeaders = malloc((*ENSData).nbAtoms * sizeof(char *))) == NULL)
                        {
                            error("MEMORY", "instance atomHeaders variable");
                            return 1;
                        }

                        if (((*ENSData).atomHeaders = malloc((*ENSData).nbAtoms * sizeof(char *))) == NULL)
                        {
                            error("MEMORY", "instance (*ENSData).atomHeaders variable");
                            return 1;
                        }

                        for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++)
                        {
                            if ((atomHeaders[atomNb] = malloc(28 * CHAR_SIZE)) == NULL)
                            {
                                error("MEMORY", "instance atomHeaders[] variable");
                                return 1;
                            }

                            if (((*ENSData).atomHeaders[atomNb] = malloc(28 * CHAR_SIZE)) == NULL)
                            {
                                error("MEMORY", "instance (*ENSData).atomHeaders[] variable");
                                return 1;
                            }
                            fread((*ENSData).atomHeaders[atomNb], CHAR_SIZE, 27, CSPFile);
                            (*ENSData).atomHeaders[atomNb][27] = '\0';

                            if ((atomNames[atomNb]  = malloc(5 * CHAR_SIZE)) == NULL)
                            {
                                error("MEMORY", "instance atomNames[] variable");
                                return 1;
                            }
                            strncpy(atomNames[atomNb], &(*ENSData).atomHeaders[atomNb][12], 4);
                            atomNames[atomNb][4] = '\0';

                            if ((residueNbs[atomNb] = malloc(5 * CHAR_SIZE)) == NULL)
                            {
                                error("MEMORY", "instance residueNbs[] variable");
                                return 1;
                            }
                            strncpy(residueNbs[atomNb], &(*ENSData).atomHeaders[atomNb][22], 4);
                            residueNbs[atomNb][4] = '\0';

                            atomHeaderPositions[0][atomNb] = atomNb;
                        }
                    }
                    else
                    {
                        if (nbAtoms != (*ENSData).nbAtoms)
                        {
                            error("ATOM_RECORD", trajectories[cspNb].fileName);
                            return 1;
                        }

                        fseek(CSPFile, (3 * CHAR_SIZE + 2 * INT32_SIZE), SEEK_SET);

                        // Sorting the atom headers based on the model
                        if ((atomHeaderPositions[cspNb] = malloc((*ENSData).nbAtoms * INT_SIZE)) == NULL)
                        {
                            error("MEMORY", "instance atomHeaderPositions[] variable");
                            return 1;
                        }

                        for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++)
                        {
                            fread(&atomHeader, CHAR_SIZE, 27, CSPFile);
                            strncpy(atomName,  &atomHeader[12], 4);
                            strncpy(residueNb, &atomHeader[22], 4);

                            if (strcmp(atomName, atomNames[atomNb]) || strcmp(residueNb, residueNbs[atomNb]))
                            {
                                found = 0;
                                for (atomSeek=0; atomSeek<(*ENSData).nbAtoms; atomSeek++)
                                {
                                    if (!(strcmp(atomName, atomNames[atomSeek]) || strcmp(residueNb, residueNbs[atomSeek])))
                                    {
                                        atomHeaderPositions[cspNb][atomNb] = atomSeek;
                                        break;
                                    }
                                }
                                if (!found)
                                {
                                    printf ("\n\n%s not found!\n\n", atomHeader);
                                    error("ATOM_RECORD", trajectories[cspNb].fileName);
                                    return 1;
                                }
                            }
                            else
                            {
                                atomHeaderPositions[cspNb][atomNb] = atomNb;
                            }
                        }
                    }
                    fclose(CSPFile);
                }
                else if (!strcmp(keyName, "first"))
                {
                    fread(&trajectories[cspNb].first, INT32_SIZE, 1, paramFile);
                }
                else if (!strcmp(keyName, "last"))
                {
                    fread(&trajectories[cspNb].last, INT32_SIZE, 1, paramFile);        // Length of one key of the hash
                }
                else if (!strcmp(keyName, "tmp_name"))
                {
                    fread(&keyLen, SHORT_SIZE, 1, paramFile);                // Length of the corresponding value
                    if ((trajectories[cspNb].tmpFileName = malloc((keyLen+1) * CHAR_SIZE)) == NULL)
                    {
                        error("MEMORY", "instance trajectories[].fileName");
                        return 1;
                    }
                    fread(trajectories[cspNb].tmpFileName, CHAR_SIZE, keyLen, paramFile);    // Recording the corresponding value
                    trajectories[cspNb].tmpFileName[keyLen] = '\0';
                }
                else if (!strcmp(keyName, "type"))
                {
                    fread(&keyLen, SHORT_SIZE, 1, paramFile);                // Length of the corresponding value
                    if ((trajectories[cspNb].type = malloc((keyLen+1) * CHAR_SIZE)) == NULL)
                    {
                        error("MEMORY", "instance trajectories[].fileName");
                        return 1;
                    }
                    fread(trajectories[cspNb].type, CHAR_SIZE, keyLen, paramFile);        // Recording the corresponding value
                    trajectories[cspNb].type[keyLen] = '\0';
                }
                else
                {
                    error("CORR_FILE", paramFileName);
                    return 1;
                }
            }
        }

        // Closing the parameter file
        fclose(paramFile);

        // Recording the coordinates
        for (confNb=0; confNb<(*ENSData).nbConfs; confNb++)
        {
            found = 0;
            for (cspNb=0; cspNb<nbCsps; cspNb++)
            {
                if (confToGet[confNb] >= trajectories[cspNb].first && confToGet[confNb] <= trajectories[cspNb].last)
                {
                    if ((CSPFile = fopen(trajectories[cspNb].fileName, "rb")) == NULL)
                    {
                        warning("NO_FILE", trajectories[cspNb].fileName);
                        return 1;
                    }
                    ++found;
                    break;
                }
            }
            if (!found)
            {
                error("One structure cannot be found in the CSP files!", NULL);
                return 1;
            }

            if (((*ENSData).atomCoords[nbRecConfs] = malloc((*ENSData).nbAtoms * COORD_SIZE))==NULL)
            {
                error("MEMORY", "instance atomCoords[] variable");
                return 1;
            }

            // Retrieving coordinates of the structure
            unsigned long confCoordStart = coordStart + confBlockSize * (confToGet[confNb] - trajectories[cspNb].first);
            for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++)
            {
                // Positionning the file handler
                fseek(CSPFile, confCoordStart + 3 * INT32_SIZE * atomHeaderPositions[cspNb][atomNb], SEEK_SET);

                // Getting x, y and z coordintates
                fread(&atomCoordInt, INT32_SIZE, 1, CSPFile);
                (*ENSData).atomCoords[nbRecConfs][atomNb].x = (float)atomCoordInt / COORD_FACT;
                fread(&atomCoordInt, INT32_SIZE, 1, CSPFile);
                (*ENSData).atomCoords[nbRecConfs][atomNb].y = (float)atomCoordInt / COORD_FACT;
                fread(&atomCoordInt, INT32_SIZE, 1, CSPFile);
                (*ENSData).atomCoords[nbRecConfs][atomNb].z = (float)atomCoordInt / COORD_FACT;
            }
            ++nbRecConfs;
            fclose(CSPFile);
        }

        // Releasing memory
        free(confToGet);
        free(ENSDirName);
        free(paramFileName);

        for (cspNb=0; cspNb<nbCsps; cspNb++)
        {
            free(trajectories[cspNb].fileName);
            free(trajectories[cspNb].tmpFileName);
            free(trajectories[cspNb].type);
            free(atomHeaderPositions[cspNb]);
        }
        free(atomHeaderPositions);
        free(trajectories);

        for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++)
        {
            free(atomNames[atomNb]);
            free(residueNbs[atomNb]);
            free(atomHeaders[atomNb]);
        }
        free(atomNames);
        free(residueNbs);
        free(atomHeaders);

        return 0;
    }
    else if (isPDB(ensFileName))
    {
        char end[PDB_TERM_LEN+1];            // To detect END or TER
        char headerType[PDB_HEADER_LEN+1];        // To detect ATOM or HETATM
        char atomCoordString[PDB_COORD_LEN+1];        // Atom coordinates
        char pdbFileLine[PDB_LINE_LEN];
        unsigned int maxNbAtoms = 0;

        if ((ENSFile = fopen(ensFileName, "rt")) == NULL) return 1;

        if (((*ENSData).atomCoords = malloc(sizeof(coordinates *))) == NULL)
        {
            error("MEMORY", "reallocate (*ENSData).atomCoords variable");
            return 1;
        }
        end[PDB_TERM_LEN]              = '\0';
        headerType[PDB_HEADER_LEN]     = '\0';
        atomCoordString[PDB_COORD_LEN] = '\0';
        (*ENSData).nbConfs             = 1;
        (*ENSData).nbAtoms             = 0;
        (*ENSData).atomHeaders = NULL;        // Initialization of the atomHeaders variable
        (*ENSData).atomCoords[0] = NULL;

        while (fgets(pdbFileLine, PDB_LINE_LEN, ENSFile) != NULL)
        {
            if ((*ENSData).nbAtoms == maxNbAtoms)
            {
                maxNbAtoms += PDB_MAX_NB_ATOMS;
                if (((*ENSData).atomCoords[0] = realloc((*ENSData).atomCoords[0],  maxNbAtoms * COORD_SIZE)) == NULL)
                {
                    error("MEMORY", "reallocate (*ENSData).atomCoords[0] variable");
                    return 1;
                }

                if (((*ENSData).atomHeaders   = realloc((*ENSData).atomHeaders, maxNbAtoms * sizeof(char *))) == NULL)
                {
                    error("MEMORY", "reallocate (*ENSData).atomHeaders variable");
                    return 1;
                }
            }

            if (strlen(pdbFileLine) > 6)
            {
                strncpy(headerType, pdbFileLine, 6);
                if (!strcmp(headerType, "ATOM  "))
                {
                    // Recording the header
                    strncpy(atomHeader, pdbFileLine, 27);
                    if (((*ENSData).atomHeaders[(*ENSData).nbAtoms] = malloc(28 * CHAR_SIZE)) == NULL)
                    {
                        error("MEMORY", "instance (*ENSData).atomHeaders[] variable");
                        return 1;
                    }
                    strcpy((*ENSData).atomHeaders[(*ENSData).nbAtoms], atomHeader);

                    // Reading the coordinates
                    strncpy(atomCoordString, &pdbFileLine[PDB_X_INDEX], PDB_COORD_LEN);
                    (*ENSData).atomCoords[0][(*ENSData).nbAtoms].x = atof(atomCoordString);
                    strncpy(atomCoordString, &pdbFileLine[PDB_Y_INDEX], PDB_COORD_LEN);
                    (*ENSData).atomCoords[0][(*ENSData).nbAtoms].y = atof(atomCoordString);
                    strncpy(atomCoordString, &pdbFileLine[PDB_Z_INDEX], PDB_COORD_LEN);
                    (*ENSData).atomCoords[0][(*ENSData).nbAtoms].z = atof(atomCoordString);

                    // Incrementing the number of atoms
                    ++(*ENSData).nbAtoms;
                }
                else if (!(strcmp(headerType, "ENDMDL") && strcmp(headerType, "END   ")))
                {
                    break;
                }
                else if (strcmp(headerType, "MODEL ") && strcmp(headerType, "ANISOU") && \
                     strcmp(headerType, "HEADER") && strcmp(headerType, "OBSLTE") && \
                     strcmp(headerType, "TITLE ") && strcmp(headerType, "SPLIT ") && \
                     strcmp(headerType, "CAVEAT") && strcmp(headerType, "COMPND") && \
                     strcmp(headerType, "SOURCE") && strcmp(headerType, "KEYWDS") && \
                     strcmp(headerType, "EXPDTA") && strcmp(headerType, "NUMMDL") && \
                     strcmp(headerType, "MDLTYP") && strcmp(headerType, "AUTHOR") && \
                     strcmp(headerType, "REVDAT") && strcmp(headerType, "SPRSDE") && \
                     strcmp(headerType, "JRNL  ") && strcmp(headerType, "REMARK") && \
                     strcmp(headerType, "DBREF ") && strcmp(headerType, "DBREF1") && \
                     strcmp(headerType, "DBREF2") && strcmp(headerType, "SEQADV") && \
                     strcmp(headerType, "SEQRES") && strcmp(headerType, "MODRES") && \
                     strcmp(headerType, "HET   ") && strcmp(headerType, "HETNAM") && \
                     strcmp(headerType, "HETSYN") && strcmp(headerType, "FORMUL") && \
                     strcmp(headerType, "HELIX ") && strcmp(headerType, "SHEET ") && \
                     strcmp(headerType, "SSBOND") && strcmp(headerType, "LINK  ") && \
                     strcmp(headerType, "CISPEP") && strcmp(headerType, "SITE  ") && \
                     strcmp(headerType, "CRYST1") && strcmp(headerType, "ORIGX1") && \
                     strcmp(headerType, "ORIGX2") && strcmp(headerType, "ORIGX3") && \
                     strcmp(headerType, "SCALE1") && strcmp(headerType, "SCALE2") && \
                     strcmp(headerType, "SCALE3") && strcmp(headerType, "MTRIX1") && \
                     strcmp(headerType, "MTRIX2") && strcmp(headerType, "MTRIX3") && \
                     strcmp(headerType, "CONECT") && strcmp(headerType, "MASTER"))
                {
                    printf("%s keyword found!\n", headerType);
                    error("PDB_FILE", ensFileName);
                    return 1;
                }
            }
            else if (strlen(pdbFileLine) > 3)        // "END" or "TER" ?
            {
                strncpy(end, pdbFileLine, PDB_TERM_LEN);
                if (!(strcmp(end, "END") && strcmp(end, "TER")))
                {
                    break;
                }
            }
        }

        (*ENSData).atomCoords[0] = realloc((*ENSData).atomCoords[0], (*ENSData).nbAtoms * COORD_SIZE);
        (*ENSData).atomHeaders   = realloc((*ENSData).atomHeaders, (*ENSData).nbAtoms * sizeof(char *));
    }
    else            // Presumably a PDB list file
    {
        // Declaration of necessary variables
        char end[PDB_TERM_LEN+1];            // To detect END or TER
        char headerType[PDB_HEADER_LEN+1];        // To detect ATOM or HETATM
        char atomCoordString[PDB_COORD_LEN+1];        // Atom coordinates
        char PDBFileName[MAX_PATHWAY_LEN], pdbFileLine[PDB_LINE_LEN];

        unsigned int maxNbAtoms = 0, maxNbConfs = 0;

        FILE *PDBFile = NULL;

        //if (((*ENSData).atomCoords = malloc(COORD_SIZE))==NULL) fatal_error("MEMORY", "instance (*ENSData).atomCoords variable");

        end[PDB_TERM_LEN]              = '\0';
        headerType[PDB_HEADER_LEN]     = '\0';
        atomCoordString[PDB_COORD_LEN] = '\0';
        (*ENSData).nbConfs             = 0;
        (*ENSData).nbAtoms             = 0;

        if ((ENSFile = fopen(ensFileName, "rt")) == NULL) return 1;

        while (fgets(PDBFileName, MAX_PATHWAY_LEN, ENSFile) != NULL)
        {
            if (PDBFileName[0] != '!' && PDBFileName[0] != '#')
            {
                if (PDBFileName[strlen(PDBFileName)-1] == '\n') PDBFileName[strlen(PDBFileName)-1] = '\0';
                if ((PDBFile = fopen(PDBFileName, "rt")) != NULL)        // Checking whether the file exists
                {
                    if ((*ENSData).nbConfs == maxNbConfs)
                    {
                        maxNbConfs += PDB_LIST_MAX_NB;
                        if (((*ENSData).atomCoords  = realloc((*ENSData).atomCoords,  maxNbConfs * sizeof(coordinates *))) == NULL)
                        {
                            error("MEMORY", "reallocate (*ENSData).atomCoords variable");
                            return 1;
                        }
                    }

                    if ((*ENSData).nbConfs)
                    {
                        if (((*ENSData).atomCoords[(*ENSData).nbConfs] = malloc((*ENSData).nbAtoms * COORD_SIZE)) == NULL)
                        {
                            error("MEMORY", "instance (*ENSData).atomCoords[] variable");
                            return 1;
                        }
                    }
                    else
                    {
                        (*ENSData).atomCoords[0] = NULL;
                        (*ENSData).atomHeaders = NULL;        // Initialization of the atomHeaders variable
                    }

                    atomNb = 0;
                    while (fgets(pdbFileLine, PDB_LINE_LEN, PDBFile) != NULL)
                    {
                        if ((!(*ENSData).nbConfs) && ((*ENSData).nbAtoms == maxNbAtoms))
                        {
                            maxNbAtoms += PDB_MAX_NB_ATOMS;
                            if (((*ENSData).atomCoords[0] = realloc((*ENSData).atomCoords[0],  maxNbAtoms * COORD_SIZE)) == NULL)
                            {
                                error("MEMORY", "reallocate (*ENSData).atomCoords[0] variable");
                                return 1;
                            }
                            //if (((*ENSData).atomCoords[0] = malloc(maxNbAtoms * COORD_SIZE)) == NULL)  fatal_error("MEMORY", "reallocate (*ENSData).atomCoords[0] variable");
                            if (((*ENSData).atomHeaders   = realloc((*ENSData).atomHeaders, maxNbAtoms * sizeof(char *))) == NULL)
                            {
                                error("MEMORY", "reallocate (*ENSData).atomHeaders variable");
                                return 1;
                            }
                        }

                        if (strlen(pdbFileLine) > 6)
                        {
                            strncpy(headerType, pdbFileLine, 6);
                            if (!strcmp(headerType, "ATOM  "))
                            {
                                strncpy(atomHeader, pdbFileLine, 27);
                                if (!(*ENSData).nbConfs)
                                {
                                    // Recording the header
                                    if (((*ENSData).atomHeaders[atomNb] = malloc(28 * CHAR_SIZE)) == NULL)
                                    {
                                        error("MEMORY", "instance (*ENSData).atomHeaders[] variable");
                                        return 1;
                                    }
                                    strcpy((*ENSData).atomHeaders[atomNb], atomHeader);
                                }
                                else
                                {
                                    // Comparing the header with the model
                                    if (strcmp(atomHeader, (*ENSData).atomHeaders[atomNb]))
                                    {
                                        error("ATOM_RECORD", PDBFileName);
                                        break;
                                    }
                                }

                                // Reading the coordinates
                                strncpy(atomCoordString, &pdbFileLine[PDB_X_INDEX], PDB_COORD_LEN);
                                (*ENSData).atomCoords[(*ENSData).nbConfs][atomNb].x = atof(atomCoordString);
                                strncpy(atomCoordString, &pdbFileLine[PDB_Y_INDEX], PDB_COORD_LEN);
                                (*ENSData).atomCoords[(*ENSData).nbConfs][atomNb].y = atof(atomCoordString);
                                strncpy(atomCoordString, &pdbFileLine[PDB_Z_INDEX], PDB_COORD_LEN);
                                (*ENSData).atomCoords[(*ENSData).nbConfs][atomNb].z = atof(atomCoordString);

                                // Incrementing the number of atoms
                                ++atomNb;
                            }
                            else if (!(strcmp(headerType, "ENDMDL") && strcmp(headerType, "END   ")))
                            {
                                break;
                            }
                            else if (strcmp(headerType, "MODEL ") && strcmp(headerType, "ANISOU") && \
                                 strcmp(headerType, "HEADER") && strcmp(headerType, "OBSLTE") && \
                                 strcmp(headerType, "TITLE ") && strcmp(headerType, "SPLIT ") && \
                                 strcmp(headerType, "CAVEAT") && strcmp(headerType, "COMPND") && \
                                 strcmp(headerType, "SOURCE") && strcmp(headerType, "KEYWDS") && \
                                 strcmp(headerType, "EXPDTA") && strcmp(headerType, "NUMMDL") && \
                                 strcmp(headerType, "MDLTYP") && strcmp(headerType, "AUTHOR") && \
                                 strcmp(headerType, "REVDAT") && strcmp(headerType, "SPRSDE") && \
                                 strcmp(headerType, "JRNL  ") && strcmp(headerType, "REMARK") && \
                                 strcmp(headerType, "DBREF ") && strcmp(headerType, "DBREF1") && \
                                 strcmp(headerType, "DBREF2") && strcmp(headerType, "SEQADV") && \
                                 strcmp(headerType, "SEQRES") && strcmp(headerType, "MODRES") && \
                                 strcmp(headerType, "HET   ") && strcmp(headerType, "HETNAM") && \
                                 strcmp(headerType, "HETSYN") && strcmp(headerType, "FORMUL") && \
                                 strcmp(headerType, "HELIX ") && strcmp(headerType, "SHEET ") && \
                                 strcmp(headerType, "SSBOND") && strcmp(headerType, "LINK  ") && \
                                 strcmp(headerType, "CISPEP") && strcmp(headerType, "SITE  ") && \
                                 strcmp(headerType, "CRYST1") && strcmp(headerType, "ORIGXn") && \
                                 strcmp(headerType, "SCALEn") && strcmp(headerType, "MTRIXn") && \
                                 strcmp(headerType, "CONECT") && strcmp(headerType, "MASTER"))
                            {
                                error("PDB_FILE", headerType);
                                break;
                            }
                        }
                        else if (strlen(pdbFileLine) > 3)        // "END"
                        {
                            strncpy(end, pdbFileLine, PDB_TERM_LEN);
                            if (!(strcmp(end, "END") && strcmp(end, "TER")))
                            {
                                break;
                            }
                        }
                    }
                    fclose(PDBFile);

                    if (!(*ENSData).nbConfs)
                    {
                        (*ENSData).nbAtoms = atomNb;
                        (*ENSData).atomCoords[0] = realloc((*ENSData).atomCoords[0], (*ENSData).nbAtoms * COORD_SIZE);
                        (*ENSData).atomHeaders   = realloc((*ENSData).atomHeaders, (*ENSData).nbAtoms * sizeof(char *));
                    }

                    if (PDBError)
                    {
                        // Releasing memory
                        for (confNb=0; confNb<(*ENSData).nbConfs; confNb++) free((*ENSData).atomCoords[confNb]); free((*ENSData).atomCoords);
                        for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++) free((*ENSData).atomHeaders[atomNb]); free((*ENSData).atomHeaders);

                        return 1;
                    }

                    ++(*ENSData).nbConfs;
                }
                else
                {
                    printf ("\n===================================\n");
                    error("NO_FILE", PDBFileName);
                    return 1;
                }
            }
        }
        if (!feof(ENSFile))
        {
            warning("The end of the file %s has not been reached", ensFileName);
        }
    }

    fclose(ENSFile);
    return 0;
}


int getAtomCoords(PDBListData *ENSData, char *ensFileName, char *atomToKeep)
{
    /*
       This function generates retrieves the coordinates found in ensFileName
       and selects only the ones that correspond to the atomToKeep
    */

    counter count, confNb, atomNb;
    char *atomName = malloc(5 * CHAR_SIZE);
    char *atomToKeepCopy = strdup(atomToKeep);
    uint32_t *atomsToKeep, nbAtomsToKeep = 0;

    // Getting all coordinates
    if (ensFileName != NULL)
    {
        if (getCoords(ENSData, ensFileName, NULL))
        {
            free(atomName); atomName = NULL;
            free(atomToKeepCopy); atomToKeepCopy = NULL;
            return 1;
        }
    }

    // Determining the index of the atoms to keep in the headers
    deepTrim(atomToKeepCopy);
    if ((atomsToKeep = malloc((*ENSData).nbAtoms * INT32_SIZE)) == NULL)
    {
        free(atomName); atomName = NULL;
        free(atomToKeepCopy); atomToKeepCopy = NULL;
        error("MEMORY", "instance atomsToKeep variable");
        return 1;
    }

    atomName[4] = '\0';
    for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++)
    {
        strncpy(atomName, &(*ENSData).atomHeaders[atomNb][12], 4);
        deepTrim(atomName);

        if (!strcmp(atomName, atomToKeepCopy))
        {
            for (count=0; count<28; count++)
            {
                (*ENSData).atomHeaders[nbAtomsToKeep][count] = (*ENSData).atomHeaders[atomNb][count];
            }
            atomsToKeep[nbAtomsToKeep] = atomNb;
            ++nbAtomsToKeep;
        }
    }

    if ((atomsToKeep = realloc(atomsToKeep, nbAtomsToKeep * INT32_SIZE)) == NULL)
    {
        free(atomName); atomName = NULL;
        free(atomToKeepCopy); atomToKeepCopy = NULL;
        return 1;
    }

    // Reallocating memory
    for (atomNb=nbAtomsToKeep; atomNb<(*ENSData).nbAtoms; atomNb++) free((*ENSData).atomHeaders[atomNb]);
    if (((*ENSData).atomHeaders = realloc((*ENSData).atomHeaders, nbAtomsToKeep * sizeof(char *))) == NULL)
    {
        free(atomName); atomName = NULL;
        free(atomsToKeep); atomsToKeep = NULL;
        free(atomToKeepCopy); atomToKeepCopy = NULL;
        return 1;
    }

    // Modifying the number of atoms
    (*ENSData).nbAtoms = nbAtomsToKeep;

    // Keeping only coordinates of atomToKeep atom
    for (confNb=0; confNb<(*ENSData).nbConfs; confNb++)
    {
        for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++)
        {
            (*ENSData).atomCoords[confNb][atomNb] = (*ENSData).atomCoords[confNb][atomsToKeep[atomNb]];
        }

        if (((*ENSData).atomCoords[confNb] = realloc((*ENSData).atomCoords[confNb], (*ENSData).nbAtoms * COORD_SIZE)) == NULL)
        {
            free(atomName); atomName = NULL;
            free(atomsToKeep); atomsToKeep = NULL;
            free(atomToKeepCopy); atomToKeepCopy = NULL;
            return 1;
        }
    }

    // Releasing memory
    free(atomName); atomName = NULL;
    free(atomsToKeep); atomsToKeep = NULL;
    free(atomToKeepCopy); atomToKeepCopy = NULL;

    return 0;
}


int getConfsCoords(PDBListData *ENSData, char *ensFileName, uintlist *confList)
{
    /*
       Retrieve the coordinates from ensFileName for confomer numbers listed in confList.
       confList must be a variable of type uintlist and sorted
    */

    counter confNb;

    // Getting all coordinates
    if (getCoords(ENSData, ensFileName, NULL)) return 1;

    // Keeping only coordinates of confList conformers
    for (confNb=0; confNb<(*confList).nbElements; confNb++)
    {
        if ((*confList).list[confNb]>(*ENSData).nbConfs || (*confList).list[confNb] < 1)
        {
            error_msg("The list contains wrong conformer numbers!", ensFileName);
            return 1;
        }
        else
        {
            (*ENSData).atomCoords[confNb] = (*ENSData).atomCoords[(*confList).list[confNb]];
        }
    }

    // Modifying the number of conformers
    (*ENSData).nbConfs = (*confList).nbElements;

    // Reducing the size of the variable that contains all coordinates
    (*ENSData).atomCoords = realloc((*ENSData).atomCoords, (*ENSData).nbConfs * sizeof(coordinates *));

    return 0;
}


int getNbResidues(PDBListData *ENSData)
{
    counter atomNb;
    unsigned int resiSeek = 0;
    char currResidueNb[5] = {"\0"};
    char prevResidueNb[5] = {"\0"};

    for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++)
    {
        strncpy(currResidueNb, &(*ENSData).atomHeaders[atomNb][22], 4);
        if (strcmp(currResidueNb, prevResidueNb))
        {
            strcpy(prevResidueNb, currResidueNb);
            ++resiSeek;
        }
    }

    return resiSeek;
}


char resiCode(char *residueName)
{
    if (!strcmp(residueName, "ALA")) return 'A';
    if (!strcmp(residueName, "CYS")) return 'C';
    if (!strcmp(residueName, "ASP")) return 'D';
    if (!strcmp(residueName, "GLU")) return 'E';
    if (!strcmp(residueName, "PHE")) return 'F';
    if (!strcmp(residueName, "GLY")) return 'G';
    if (!strcmp(residueName, "HIS")) return 'H';
    if (!strcmp(residueName, "ILE")) return 'I';
    if (!strcmp(residueName, "LYS")) return 'K';
    if (!strcmp(residueName, "LEU")) return 'L';
    if (!strcmp(residueName, "MET")) return 'M';
    if (!strcmp(residueName, "ASN")) return 'N';
    if (!strcmp(residueName, "PRO")) return 'P';
    if (!strcmp(residueName, "GLN")) return 'Q';
    if (!strcmp(residueName, "ARG")) return 'R';
    if (!strcmp(residueName, "SER")) return 'S';
    if (!strcmp(residueName, "THR")) return 'T';
    if (!strcmp(residueName, "VAL")) return 'V';
    if (!strcmp(residueName, "TRP")) return 'W';
    if (!strcmp(residueName, "TYR")) return 'Y';

    return 'X';
}


void pdb2seq(char **sequence, PDBListData *ENSData)
{
    counter atomHeaderNb;

    char prevResiNb[5] = {'\0'}, resiNb[5] = {'\0'}, residueName[4] = {'\0'};

    unsigned int nbResidues = 0, maxSeqLen = 0;

    for (atomHeaderNb=0; atomHeaderNb<(*ENSData).nbAtoms; atomHeaderNb++)
    {
        strncpy(resiNb, &(*ENSData).atomHeaders[atomHeaderNb][22], 4);
        if (strcmp(resiNb, prevResiNb))
        {
            if (nbResidues == maxSeqLen)
            {
                maxSeqLen += PDB_MAX_NB_RESIDUES;
                (*sequence) = realloc((*sequence), (maxSeqLen + 1) * CHAR_SIZE);    // +1 for the NULL terminating character
            }

            // Retrieving the 3L residue name
            strncpy(residueName, &(*ENSData).atomHeaders[atomHeaderNb][17], 3);
            (*sequence)[nbResidues] = resiCode(residueName);

            // Considering now the new residue number as the previous residue
            strcpy(prevResiNb, resiNb);

            // incrementing the number of residues
            ++nbResidues;
        }
    }

    (*sequence) = realloc((*sequence), (nbResidues+1) * CHAR_SIZE);
    (*sequence)[nbResidues] = '\0';
}


int getResidueNumbers(residueList *residueNbs, PDBListData *ENSData)
{
    /*
       This function does not return any value. Instead, it does fill
       the residueNbs variable with the number of residues and the
       residue numbers. Then, it is easy to verify that everything went
       alright by checking the number of residues.
    */

    counter atomNb;
    char *residueNbStr = malloc(5 * CHAR_SIZE);
    char previousResidueNbStr[5] = {'\0'};

    if (((*residueNbs).nbs = malloc((*ENSData).nbAtoms * INT32_SIZE)) == NULL)
    {
        free(residueNbStr);
        error("MEMORY", "instance residueNbs variable");
        return 1;
    }

    (*residueNbs).nb = 0;
    for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++)
    {
        strncpy(residueNbStr, &(*ENSData).atomHeaders[atomNb][22], 4);
        residueNbStr[4] = '\0';
        deepTrim(residueNbStr);
        if (strcmp(residueNbStr, previousResidueNbStr))
        {
            strcpy(previousResidueNbStr, residueNbStr);
            (*residueNbs).nbs[(*residueNbs).nb] = atoi(residueNbStr);
            ++(*residueNbs).nb;
        }
    }

    // Releasing memory
    free(residueNbStr);

    if (((*residueNbs).nbs = realloc((*residueNbs).nbs, (*residueNbs).nb*INT32_SIZE)) == NULL) return 1;

    return 0;
}


int WritePDBFile(char *PDBFileName, PDBListData *ENSData, unsigned int confNb)
{
    counter atomNb;
    FILE *PDBFile;

    if (confNb>(*ENSData).nbConfs) return 1;
    if ((PDBFile = fopen(PDBFileName, "w")) == NULL) return 1;

    for (atomNb=0; atomNb<(*ENSData).nbAtoms; atomNb++)
    {
        fprintf(PDBFile, "%s   %8.3f%8.3f%8.3f  1.00  0.00\n", (*ENSData).atomHeaders[atomNb], \
                                                               (*ENSData).atomCoords[confNb][atomNb].x,\
                                                               (*ENSData).atomCoords[confNb][atomNb].y,\
                                                               (*ENSData).atomCoords[confNb][atomNb].z);
    }
    fprintf(PDBFile, "END\n");
    fclose(PDBFile);

    return 0;
}

#endif