开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 并行计算 > 查看源码
InputDataLoader.cpp
资源名称：fastgrid-4.2.1-141-src.tar.gz　[点击查看]
上传用户：chinafayin
上传日期：2022-04-05
资源大小：153k
文件大小：42k
源码类别：

并行计算

开发平台：

Visual C++

InputDataLoader.cpp：源码内容
  1. /*
  2.     FastGrid (formerly AutoGrid)
  4.     Copyright (C) 2009 The Scripps Research Institute. All rights reserved.
  5.     Copyright (C) 2009 Masaryk University. All rights reserved.
  7.     AutoGrid is a Trade Mark of The Scripps Research Institute.
  9.     This program is free software; you can redistribute it and/or
  10.     modify it under the terms of the GNU General Public License
  11.     as published by the Free Software Foundation; either version 2
  12.     of the License, or (at your option) any later version.
  14.     This program is distributed in the hope that it will be useful,
  15.     but WITHOUT ANY WARRANTY; without even the implied warranty of
  16.     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17.     GNU General Public License for more details.
  19.     You should have received a copy of the GNU General Public License
  20.     along with this program; if not, write to the Free Software
  21.     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  22. */
  24. #include <cstring>
  25. #include <cctype>
  26. #include <climits>
  27. #include "InputDataLoader.h"
  28. #include "Utils.h"
  29. #include "Exceptions.h"
  31. // GPF tokens
  32. enum GPFTokens
  33. {
  34.     GPF_NULL = 0,
  35.     GPF_COMMENT,
  36.     GPF_RECEPTOR,
  37.     GPF_GRIDFLD,
  38.     GPF_NPTS,
  39.     GPF_SPACING,
  40.     GPF_GRIDCENTER,
  41.     GPF_LIGAND_TYPES,
  42.     GPF_MAP,
  43.     GPF_NBP_COEFFS,
  44.     GPF_NBP_R_EPS,
  45.     GPF_ELECMAP,
  46.     GPF_DIEL,
  47.     GPF_FMAP,
  48.     GPF_DISORDER,
  49.     GPF_SMOOTH,
  50.     GPF_SOL_PAR,
  51.     GPF_CONSTANT,
  52.     GPF_COVALENTMAP,
  53.     GPF_RECEPTOR_TYPES,
  54.     GPF_PARAM_FILE,
  55.     GPF_DSOLVMAP,
  56.     GPF_QASP
  57. };
  59. InputDataLoader::InputDataLoader(LogFile *logFile): logFile(logFile)
  60. {
  61.     fldFilenameAVS[0] = 0;
  62.     floatingGridFilename[0] = 0;
  63.     receptorFilename[0] = 0;
  64.     xyzFilename[0] = 0;
  65.     parameterLibraryFilename[0] = 0;
  67.     memset(receptorTypes, 0, sizeof(receptorTypes));
  68.     numReceptorTypes = 0;
  69.     numGridPointsPerMap = INIT_NUM_GRID_PTS;
  70.     numReceptorAtoms = 0;
  71.     covalentPoint = 0;
  73.     // for NEW3 desolvation terms
  74.     solparQ = .01097;   // unweighted value restored 3:9:05
  75.     invDielCal = 0;
  76.     rSmooth = 0;
  77.     gridSpacing = 0.375;     // One quarter of a C-C bond length.
  78.     covHalfWidthSquaredInv = 1.0;
  79.     covBarrier = 1000.0;
  81.     distDepDiel = false;
  82.     disorderH = false;
  84.     receptorAtom = 0;
  85. }
  87. InputDataLoader::~InputDataLoader()
  88. {
  89.     if (receptorAtom)
  90.     {
  91.         delete [] charge;
  92.         delete [] vol;
  93.         delete [] solpar;
  94.         delete [] atomType;
  95.         delete [] hbond;
  96.         delete [] receptorAtom;
  97.     }
  98. }
  100. void InputDataLoader::load(const char *gridParameterFilename, GridMapList &gridmaps, ParameterLibrary &parameterLibrary)
  101. {
  102.     // LIGAND: maximum is MAX_MAPS
  103.     // each type is now at most two characters plus ''
  104.     // currently ligandAtomTypes is sparse... some types are not set
  105.     char ligandTypes[MAX_MAPS][3] = {0};
  107.     // number of different receptor atom types declared on receptorTypes line in GPF
  108.     int receptor_types_gpf_ct = 0;
  109.     int has_receptor_types_in_gpf = 0;
  111.     // array of numbers of each type
  112.     // NB: this is a sparse int array, some entries are 0
  113.     int receptor_atom_type_count[NUM_RECEPTOR_TYPES] = {0};
  115.     Vec3d gridExtent;
  116.     Vec3d gridCornerMax;
  117.     Vec3d cmax = -BIG;
  118.     Vec3d cmin = BIG;
  119.     Vec3d csum = 0;
  120.     Vec3d cmean;
  122.     // LINE_LEN
  123.     char line[LINE_LEN];
  124.     char GPFLine[LINE_LEN];
  125.     int length = LINE_LEN;
  127.     char atom_name[6];
  128.     char record[LINE_LEN];
  129.     char temp_char = ' ';
  130.     char token[LINE_LEN];
  131.     static const char xyz[] = "xyz"; // used to print headings
  132.     FILE *receptorFile, *xyz_fileptr = 0;
  134.     double q_tot = 0.0;
  135.     double diel;
  136.     double q_max = -BIG, q_min = BIG;
  137.     double ri;
  138.     double temp_vol, temp_solpar;
  140.     int GPF_keyword = -1;
  141.     int indcom = 0;
  142.     int infld;
  143.     int mapIndex = -1;
  145.     // Initializes the grid parameter file
  146.     FILE *GPF = stdin;
  147.     if (gridParameterFilename[0])
  148.     {
  149.         GPF = boincOpenFile(gridParameterFilename, "r");
  150.         if (!GPF)
  151.         {
  152.             logFile->printErrorFormatted(ERROR, "Sorry, I can't find or open Grid Parameter File "%s"", gridParameterFilename);
  153.             logFile->printErrorFormatted(ERROR, "Unsuccessful Completion.n");
  154.             throw ExitProgram(911);
  155.         }
  156.     }
  158.     // Read in the grid parameter file...
  159.     ParameterEntry thisparm;
  160.     ParameterEntry *foundParam;
  162.     while (fgets(GPFLine, LINE_LEN, GPF) != 0)
  163.     {
  164.         GPF_keyword = parseGPFLine(GPFLine);
  166.         // This first "switch" figures out how to echo the current GPF line.
  167.         switch (GPF_keyword)
  168.         {
  169.         case -1:
  170.             logFile->printFormatted("GPF> %s", GPFLine);
  171.             logFile->printError(WARNING, "Unrecognized keyword in grid parameter file.n");
  172.             continue;           // while fgets GPFLine...
  174.         case GPF_NULL:
  175.         case GPF_COMMENT:
  176.             logFile->printFormatted("GPF> %s", GPFLine);
  177.             break;
  179.         default:
  180.             logFile->printFormatted("GPF> %s", GPFLine);
  181.             indcom = strIndex(GPFLine, "#");
  182.             if (indcom != -1)
  183.                 GPFLine[indcom] = '';    // Truncate str. at the comment
  184.             break;
  186.         }                       // first switch
  188.         // This second switch interprets the current GPF line.
  189.         switch (GPF_keyword)
  190.         {
  191.         case GPF_NULL:
  192.         case GPF_COMMENT:
  193.             break;
  195.         case GPF_RECEPTOR:
  196.             {
  197.                 // read in the receptor gridParameterFilename
  198.                 sscanf(GPFLine, "%*s %s", receptorFilename);
  199.                 logFile->printFormatted("nReceptor Input File :t%snnReceptor Atom Type Assignments:nn", receptorFilename);
  201.                 // try to open receptor file
  202.                 if ((receptorFile = boincOpenFile(receptorFilename, "r")) == 0)
  203.                 {
  204.                     logFile->printErrorFormatted(ERROR, "can't find or open receptor PDBQT file "%s".n", receptorFilename);
  205.                     logFile->printError(FATAL_ERROR, "Unsuccessful completion.nn");
  206.                 }
  208.                 // get the upper bound of total number of atoms
  209.                 fseek(receptorFile, 0, SEEK_END);
  210.                 int sizeOfFile = ftell(receptorFile);
  211.                 fseek(receptorFile, 0, SEEK_SET);
  212.                 char *fileContent = new char[sizeOfFile];
  213.                 if (!fread(fileContent, sizeOfFile, 1, receptorFile))
  214.                     logFile->printError(FATAL_ERROR, "Cannot read the receptor file.");
  215.                 fseek(receptorFile, 0, SEEK_SET);
  217.                 int numAtomsMax = 1;
  218.                 for (int i = 0; i < sizeOfFile; i++)
  219.                     if (fileContent[i] == 'n')
  220.                         ++numAtomsMax;
  221.                 
  222.                 delete [] fileContent; 
  223.                 
  224.                 // reserve space for atoms
  225.                 charge = new double[numAtomsMax];
  226.                 vol = new double[numAtomsMax];
  227.                 solpar = new double[numAtomsMax];
  228.                 atomType = new int[numAtomsMax];
  229.                 hbond = new HBondType[numAtomsMax];
  230.                 receptorAtom = new Vec4d[numAtomsMax];
  232.                 // start to read in the lines of the receptor file
  233.                 int ia = 0;
  234.                 while ((fgets(line, length, receptorFile)) != 0)
  235.                 {
  236.                     sscanf(line, "%6s", record);
  237.                     if (strncmp(record, "ATOM", 4) == 0 || // Amino Acid or DNA/RNA atoms
  238.                         strncmp(record, "HETA", 4) == 0 || // Non-standard heteroatoms
  239.                         strncmp(record, "CHAR", 4) == 0)
  240.                     {
  241.                         // Partial Atomic Charge - not a PDB record
  243.                         strncpy(atom_name, &line[12], 4);
  244.                         /* atom_name is declared as an array of 6 characters, the PDB atom name is 4 characters (C indices 0, 1, 2 and 3) but let's ensure that the fifth character (C index 4)
  245.                            is a null character, which terminates the string. */
  246.                         atom_name[4] = '';
  248.                         // Output the serial number of this atom...
  249.                         logFile->printFormatted("Atom no. %2d, "%s"", ia + 1, atom_name);
  251.                         // Read in this receptor atom's coordinates,partial charges, and solvation parameters in PDBQS format...
  252.                         sscanf(&line[30], "%lf", &receptorAtom[ia].x);
  253.                         sscanf(&line[38], "%lf", &receptorAtom[ia].y);
  254.                         sscanf(&line[46], "%lf", &receptorAtom[ia].z);
  256.                         // Output the coordinates of this atom...
  257.                         logFile->printFormatted(" at (%.3lf, %.3lf, %.3lf), ", receptorAtom[ia].x, receptorAtom[ia].y, receptorAtom[ia].z);
  259.                         // 1:CHANGE HERE: need to set up vol and solpar
  260.                         sscanf(&line[70], "%lf", &charge[ia]);
  261.                         sscanf(&line[77], "%s", thisparm.autogridType);
  262.                         foundParam = parameterLibrary.findAtomParameter(thisparm.autogridType);
  263.                         if (foundParam != 0)
  264.                         {
  265.                             logFile->printFormatted("DEBUG: foundParam->recIndex = %d", foundParam->recIndex);
  266.                             if (foundParam->recIndex < 0)
  267.                             {
  268.                                 strncpy(receptorTypes[numReceptorTypes], foundParam->autogridType, 3);
  269.                                 foundParam->recIndex = numReceptorTypes++;
  270.                                 logFile->printFormatted("DEBUG: foundParam->recIndex => %d", foundParam->recIndex);
  271.                             }
  272.                             atomType[ia] = foundParam->recIndex;
  273.                             solpar[ia] = foundParam->solpar;
  274.                             vol[ia] = foundParam->vol;
  275.                             hbond[ia] = foundParam->hbond;  // NON=0, DS,D1, AS, A1, A2
  276.                             ++receptor_atom_type_count[foundParam->recIndex];
  277.                         }
  278.                         else
  279.                         {
  280.                             logFile->printFormatted("nnreceptor file contains unknown type: '%s'nadd parameters for it to the parameter library firstn", thisparm.autogridType);
  281.                             throw ExitProgram(-1);
  282.                         }
  284.                         // if from pdbqs: convert cal/molA**3 to kcal/molA**3
  285.                         // solpar[ia] *= 0.001;
  286.                         q_max = Mathd::Max(q_max, charge[ia]);
  287.                         q_min = Mathd::Min(q_min, charge[ia]);
  289.                         if (atom_name[0] == ' ')
  290.                         {
  291.                             // truncate the first character...
  292.                             atom_name[0] = atom_name[1];
  293.                             atom_name[1] = atom_name[2];
  294.                             atom_name[2] = atom_name[3];
  295.                             atom_name[3] = '';
  296.                         }
  297.                         else if ((atom_name[0] == '0' || atom_name[0] == '1' || atom_name[0] == '2' ||
  298.                                  atom_name[0] == '3' || atom_name[0] == '4' || atom_name[0] == '5' ||
  299.                                  atom_name[0] == '6' || atom_name[0] == '7' || atom_name[0] == '8' ||
  300.                                  atom_name[0] == '9') && atom_name[1] == 'H')
  301.                         {
  302.                             /* Assume this is the 'mangled' name of a hydrogen atom, after the atom name has been changed from 'HD21' to '1HD2' for example. [0-9]H(.)(.) 0 1 2 3 : : :
  303.                                : V V V V tmp 0 1 2 tmp : V 0 1 2 3 : : : : V V V V H(.)(.)[0-9] */
  304.                             temp_char = atom_name[0];
  305.                             atom_name[0] = atom_name[1];
  306.                             atom_name[1] = atom_name[2];
  307.                             atom_name[2] = atom_name[3];
  308.                             atom_name[3] = temp_char;
  309.                         }
  311.                         // Tell the user what you thought this atom was...
  312.                         logFile->printFormatted(" was assigned atom type "%s" (recIndex= %d, atomType= %d).n", foundParam->autogridType, foundParam->recIndex, atomType[ia]);
  314.                         // Count the number of each atom type
  315.                         // ++receptor_atom_type_count[ atomType[ia] ];
  317.                         // Keep track of the extents of the receptor
  318.                         for (int i = 0; i < 3; i++)
  319.                         {
  320.                             cmax[i] = Mathd::Max(cmax[i], receptorAtom[ia][i]);
  321.                             cmin[i] = Mathd::Min(cmin[i], receptorAtom[ia][i]);
  322.                         }
  323.                         csum += Vec3d(receptorAtom[ia]);
  324.                         // Total up the partial charges as we go...
  325.                         q_tot += charge[ia];
  327.                         // Increment the atom counter
  328.                         ia++;
  330.                         // Check that there aren't too many atoms...
  331.                         if (ia > AG_MAX_ATOMS)
  332.                         {
  333.                             logFile->printErrorFormatted(ERROR, "Too many atoms in receptor PDBQT file %s;", receptorFilename);
  334.                             logFile->printErrorFormatted(ERROR, "-- the maximum number of atoms, AG_MAX_ATOMS, allowed is %ul.", AG_MAX_ATOMS);
  335.                             logFile->printErrorFormatted(SUGGESTION, "Increase the value in the "#define AG_MAX_ATOMS %ul" line", AG_MAX_ATOMS);
  336.                             logFile->printError(SUGGESTION, "in the source file "autogrid.h", and re-compile " APPNAME ".");
  337.                             // FATAL_ERROR will cause this program to exit...
  338.                             logFile->printError(FATAL_ERROR, "Sorry, " APPNAME " cannot continue.");
  339.                         }           // endif
  340.                     }               // endif
  341.                 }                   // endwhile
  342.                 // Finished reading in the lines of the receptor file
  343.                 fclose(receptorFile);
  344.                 if (has_receptor_types_in_gpf == 1)
  345.                     // Check that the number of atom types found in the receptor PDBQT
  346.                     // file match the number parsed in by the "receptorTypes" command
  347.                     // in the GPF; if they do not match, exit!
  348.                     if (numReceptorTypes != receptor_types_gpf_ct)
  349.                     {
  350.                         logFile->printErrorFormatted(ERROR,
  351.                             "The number of atom types found in the receptor PDBQT (%d) does not match the number specified by the "receptor_types" command (%d) in the GPF!nn",
  352.                             numReceptorTypes, receptor_types_gpf_ct);
  353.                         // FATAL_ERROR will cause this program to exit...
  354.                         logFile->printError(FATAL_ERROR, "Sorry, " APPNAME " cannot continue.");
  355.                     }
  357.                 // Update the total number of atoms in the receptor
  358.                 numReceptorAtoms = ia;
  359.                 logFile->printFormatted("nMaximum partial atomic charge found = %+.3lf en", q_max);
  360.                 logFile->printFormatted("Minimum partial atomic charge found = %+.3lf enn", q_min);
  362.                 // Check there are partial charges...
  363.                 if (q_max == 0 && q_min == 0)
  364.                 {
  365.                     logFile->printErrorFormatted(ERROR, "No partial atomic charges were found in the receptor PDBQT file %s!nn", receptorFilename);
  366.                     // FATAL_ERROR will cause this program to exit...
  367.                     logFile->printError(FATAL_ERROR, "Sorry, " APPNAME " cannot continue.");
  368.                 }                   // if there are no charges EXIT
  370.                 logFile->print("AtomtAtomtNumber of this Typen"
  371.                               "Typet ID t in Receptorn"
  372.                               "____t____t___________________n");
  374.                 // 2. CHANGE HERE: need to count number of each receptor_type
  375.                 for (int ia = 0; ia < numReceptorTypes; ia++)
  376.                     if (receptor_atom_type_count[ia] != 0)
  377.                         logFile->printFormatted(" %dt %stt%6dn", (ia), receptorTypes[ia], receptor_atom_type_count[ia]);
  379.                 logFile->printFormatted("nTotal number of atoms :tt%d atoms n"
  380.                                        "Total charge :ttt%.2lf en"
  381.                                        "nnReceptor coordinates fit within the following volume:nn"
  382.                                        "                   _______(%.1lf, %.1lf, %.1lf)n"
  383.                                        "                  /|     /|n"
  384.                                        "                 / |    / |n"
  385.                                        "                /______/  |n"
  386.                                        "                |  |___|__| Midpoint = (%.1lf, %.1lf, %.1lf)n"
  387.                                        "                |  /   |  /n"
  388.                                        "                | /    | /n"
  389.                                        "                |/_____|/n"
  390.                                        "(%.1lf, %.1lf, %.1lf)      n"
  391.                                        "nMaximum coordinates :tt(%.3lf, %.3lf, %.3lf)n"
  392.                                        "Minimum coordinates :tt(%.3lf, %.3lf, %.3lf)nn"
  393.                                        "n",
  394.                                        numReceptorAtoms, q_tot, cmax.x, cmax.y, cmax.z,
  395.                                        (cmax.x + cmin.x) / 2, (cmax.y + cmin.y) / 2, (cmax.z + cmin.z) / 2,
  396.                                        cmin.x, cmin.y, cmin.z, cmax.x, cmax.y, cmax.z, cmin.x, cmin.y, cmin.z);
  398.                 cmean[0] = csum[0] / numReceptorAtoms;
  399.                 cmean[1] = csum[1] / numReceptorAtoms;
  400.                 cmean[2] = csum[2] / numReceptorAtoms;
  401.             }
  402.             break;
  404.         case GPF_GRIDFLD:
  405.             sscanf(GPFLine, "%*s %s", fldFilenameAVS);
  406.             infld = strIndex(fldFilenameAVS, ".fld");
  407.             if (infld == -1)
  408.                 logFile->printError(FATAL_ERROR, "Grid data file needs the extension ".fld" for AVS inputnn");
  409.             else
  410.             {
  411.                 infld = strIndex(fldFilenameAVS, "fld");
  412.                 strncpy(xyzFilename, fldFilenameAVS, MAX_CHARS);
  413.                 xyzFilename[infld] = 'x';
  414.                 xyzFilename[infld + 1] = 'y';
  415.                 xyzFilename[infld + 2] = 'z';
  416.             }
  417.             if ((xyz_fileptr = boincOpenFile(xyzFilename, "w")) == 0)
  418.             {
  419.                 logFile->printErrorFormatted(ERROR, "can't create grid extrema data file %sn", xyzFilename);
  420.                 logFile->printError(ERROR, "SORRY!    unable to create the ".xyz" file.nn");
  421.                 logFile->printError(FATAL_ERROR, "Unsuccessful completion.nn");
  422.             }
  423.             else
  424.                 logFile->printFormatted("nCreating (AVS-readable) grid-coordinates extrema file : %snn", xyzFilename);
  425.             break;
  427.         case GPF_NPTS:
  428.             {
  429.                 Vec3i numGridPointsMinusOne;
  430.                 sscanf(GPFLine, "%*s %d %d %d", &numGridPointsMinusOne.x, &numGridPointsMinusOne.y, &numGridPointsMinusOne.z);
  432.                 // numGridPointsMinusOne mustn't be negative, shouldn't be zero or larger than MAX_GRID_PTS and should be even
  433.                 for (int i = 0; i < 3; i++)
  434.                     numGridPointsMinusOne[i] = checkSize(numGridPointsMinusOne[i], xyz[i]);
  436.                 numGridPointsDiv2 = numGridPointsMinusOne / 2;
  437.                 numGridPoints = numGridPointsMinusOne + 1;
  438.                 logFile->printFormatted("nNumber of grid points in x-direction:t%dn"
  439.                                        "Number of grid points in y-direction:t%dn"
  440.                                        "Number of grid points in z-direction:t%dnn",
  441.                                        numGridPoints.x, numGridPoints.y, numGridPoints.z);
  443.                 numGridPointsPerMap = numGridPoints.Cube();
  444.             }
  445.             break;
  447.         case GPF_SPACING:
  448.             sscanf(GPFLine, "%*s %lf", &gridSpacing);
  449.             logFile->printFormatted("Grid Spacing :ttt%.3lf Angstromnn", gridSpacing);
  450.             break;
  452.         case GPF_GRIDCENTER:
  453.             sscanf(GPFLine, "%*s %s", token);
  454.             if (strncmp(token, "auto", 4) == 0)
  455.             {
  456.                 gridCenter = cmean;
  457.                 logFile->printFormatted("Grid maps will be centered on the center of mass.n"
  458.                                        "Coordinates of center of mass : (%.3lf, %.3lf, %.3lf)n", gridCenter.x, gridCenter.y, gridCenter.z);
  459.             }
  460.             else
  461.             {
  462.                 sscanf(GPFLine, "%*s %lf %lf %lf", &gridCenter.x, &gridCenter.y, &gridCenter.z);
  463.                 logFile->printFormatted("nGrid maps will be centered on user-defined coordinates:nntt(%.3lf, %.3lf, %.3lf)n", gridCenter.x, gridCenter.y, gridCenter.z);
  464.             }
  465.             // centering stuff...
  466.             for (int ia = 0; ia < numReceptorAtoms; ia++)
  467.                 *((Vec3d*)&receptorAtom[ia]) -= gridCenter;  // transform to center of gridmaps
  468.             gridExtent = Vec3d(numGridPointsDiv2) * gridSpacing;
  469.             gridCornerMax = gridCenter + gridExtent;
  470.             gridCornerMin = gridCenter - gridExtent;
  472.             logFile->printFormatted("nGrid maps will cover the following volume:nn"
  473.                                    "                   _______(%.1lf, %.1lf, %.1lf)n"
  474.                                    "                  /|     /|n"
  475.                                    "                 / |    / |n"
  476.                                    "                /______/  |n"
  477.                                    "                |  |___|__| Midpoint = (%.1lf, %.1lf, %.1lf)n"
  478.                                    "                |  /   |  /n"
  479.                                    "                | /    | /n"
  480.                                    "                |/_____|/n"
  481.                                    "(%.1lf, %.1lf, %.1lf)      nn", gridCornerMax.x, gridCornerMax.y, gridCornerMax.z,
  482.                                    gridCenter.x, gridCenter.y, gridCenter.z, gridCornerMin.x, gridCornerMin.y, gridCornerMin.z);
  484.             for (int i = 0; i < 3; i++)
  485.                 logFile->printFormatted("Grid map %c-dimension :tt%.1lf Angstromsn", xyz[i], 2 * gridExtent[i]);
  487.             logFile->printFormatted("nMaximum coordinates :tt(%.3lf, %.3lf, %.3lf)n"
  488.                                    "Minimum coordinates :tt(%.3lf, %.3lf, %.3lf)nn", gridCornerMax.x, gridCornerMax.y, gridCornerMax.z, gridCornerMin.x, gridCornerMin.y, gridCornerMin.z);
  489.             for (int i = 0; i < 3; i++)
  490.                 fprintf(xyz_fileptr, "%.3lf %.3lfn", gridCornerMin[i], gridCornerMax[i]);
  491.             fclose(xyz_fileptr);
  492.             break;
  494.         case GPF_LIGAND_TYPES:
  495.             {
  496.                 // Read in the list of atom types in the ligand.
  497.                 // GPFLine e.g.: "ligand_types N O A C HH NH"
  499.                 // array of ptrs used to parse input line
  500.                 char *ligandAtomTypes[MAX_MAPS];
  502.                 int numAtomMaps = parseTypes(GPFLine, ligandAtomTypes, MAX_ATOM_TYPES);
  503.                 for (int i = 0; i < numAtomMaps; i++)
  504.                     strncpy(ligandTypes[i], ligandAtomTypes[i], 3);
  505.                 for (int i = 0; i < numAtomMaps; i++)
  506.                 {
  507.                     foundParam = parameterLibrary.findAtomParameter(ligandTypes[i]);
  508.                     if (foundParam)
  509.                         foundParam->mapIndex = i;
  510.                     else
  511.                     {
  512.                         // return error here
  513.                         logFile->printFormatted("unknown ligand atom type %snadd parameters for it to the parameter library first!n", ligandAtomTypes[i]);
  514.                         throw ExitProgram(-1);
  515.                     }
  516.                 }
  518.                 // Check to see if there is enough memory to store these map objects
  519.                 gridmaps.setNumMaps(numAtomMaps+2);
  521.                 if (numGridPointsPerMap == INIT_NUM_GRID_PTS)
  522.                 {
  523.                     logFile->printError(ERROR, "You need to set the number of grid points using "npts" before setting the ligand atom types, using "ligand_types".n");
  524.                     logFile->printError(FATAL_ERROR, "Unsuccessful completion.nn");
  525.                 }
  527.                 for (int i = 0; i < numAtomMaps; i++)
  528.                 {
  529.                     strncpy(gridmaps[i].type, ligandTypes[i], 3);   // eg HD or OA or NA or N
  530.                     foundParam = parameterLibrary.findAtomParameter(ligandTypes[i]);
  531.                     if (strcmp(ligandTypes[i], "Z") == 0)
  532.                     {
  533.                         logFile->print("Found covalent map atomtypen");
  534.                         gridmaps[i].isCovalent = true;
  535.                     }
  536.                     gridmaps[i].atomType = foundParam->mapIndex;
  537.                     gridmaps[i].solparProbe = foundParam->solpar;
  538.                     gridmaps[i].volProbe = foundParam->vol;
  539.                     gridmaps[i].Rij = foundParam->Rij;
  540.                     gridmaps[i].epsij = foundParam->epsij;
  541.                     gridmaps[i].hbond = foundParam->hbond;
  542.                     gridmaps[i].RijHB = foundParam->RijHB;
  543.                     gridmaps[i].epsijHB = foundParam->epsijHB;
  544.                     if (gridmaps[i].hbond > 0)
  545.                         gridmaps[i].isHBonder = true;
  547.                     for (int j = 0; j < numReceptorTypes; j++)
  548.                     {
  549.                         foundParam = parameterLibrary.findAtomParameter(receptorTypes[j]);
  550.                         gridmaps[i].nbpR[j] = (gridmaps[i].Rij + foundParam->Rij) / 2;
  551.                         gridmaps[i].nbpEps[j] = sqrt(gridmaps[i].epsij * foundParam->epsij);
  552.                         // apply the vdW forcefield parameter/weight here
  553.                         // This was removed because "setup_p_l" does this for us... gridmaps[i].nbpEps[j] *= FE_coeff_vdW;
  554.                         gridmaps[i].xA[j] = 12;
  555.                         // setup hbond dependent stuff
  556.                         gridmaps[i].xB[j] = 6;
  557.                         gridmaps[i].hbonder[j] = 0;
  558.                         if (gridmaps[i].hbond > D1 && (foundParam->hbond == DS || foundParam->hbond == D1))
  559.                         {
  560.                             // AS,A1,A2 map vs DS,D1 probe
  561.                             gridmaps[i].xB[j] = 10;
  562.                             gridmaps[i].hbonder[j] = 1;
  563.                             gridmaps[i].isHBonder = true;
  564.                             // Rij and epsij for this hb interaction in parm_data.dat file as Rii and epsii for heavy atom hb factors
  565.                             gridmaps[i].nbpR[j] = gridmaps[i].RijHB;
  566.                             gridmaps[i].nbpEps[j] = gridmaps[i].epsijHB;
  568.                             // apply the hbond forcefield parameter/weight here
  569.                             // This was removed because "setup_p_l" does this for us... gridmaps[i].nbpEps[j] *= FE_coeff_hbond;
  570.                         }
  571.                         else if ((gridmaps[i].hbond == DS || gridmaps[i].hbond == D1) && foundParam->hbond > D1)
  572.                         {
  573.                             // DS,D1 map vs AS,A1,A2 probe
  574.                             gridmaps[i].xB[j] = 10;
  575.                             gridmaps[i].hbonder[j] = 1;
  576.                             gridmaps[i].isHBonder = true;
  577.                             // Rij and epsij for this hb interaction in parm_data.dat file as Rii and epsii for heavy atom hb factors
  578.                             gridmaps[i].nbpR[j] = foundParam->RijHB;
  579.                             gridmaps[i].nbpEps[j] = foundParam->epsijHB;
  581.                             // apply the hbond forcefield parameter/weight here
  582.                             // This was removed because "setup_p_l" does this for us... gridmaps[i].nbpEps[j] *= FE_coeff_hbond;
  583.                         }
  584.                     }              // initialize energy parms for each possible receptor type
  585.                 }                   // for each map
  587.                 logFile->printFormatted("nAtom type names for ligand atom types 1-%d used for ligand-atom affinity grid maps:nn", numAtomMaps);
  588.                 for (int i = 0; i < numAtomMaps; i++)
  589.                 {
  590.                     logFile->printFormatted("tttAtom type number %d corresponds to atom type name "%s".n", i, gridmaps[i].type);
  592.                     if (gridmaps[i].isCovalent)
  593.                         logFile->printFormatted("nAtom type number %d will be used to calculate a covalent affinity grid mapnn", i + 1);
  594.                 }
  595.                 logFile->print("nn");
  596.             }
  597.             break;
  599.         case GPF_RECEPTOR_TYPES:
  600.             {
  601.                 // Read in the list of atom types in the receptor.
  602.                 // GPFLine e.g.: "receptorTypes N O A C HH NH"
  603.                 //
  604.                 // NOTE: This line is not guaranteed to match the actual
  605.                 // atom types present in the receptor PDBQT file
  606.                 // specified by the "receptor" command.
  608.                 // array of ptrs used to parse input line
  609.                 char *receptor_atom_types[NUM_RECEPTOR_TYPES];
  611.                 numReceptorTypes = parseTypes(GPFLine, receptor_atom_types, MAX_ATOM_TYPES);
  612.                 receptor_types_gpf_ct = numReceptorTypes;
  613.                 has_receptor_types_in_gpf = 1;
  615.                 for (int i = 0; i < numReceptorTypes; i++)
  616.                     strncpy(receptorTypes[i], receptor_atom_types[i], 3);
  617.                 for (int i = 0; i < numReceptorTypes; i++)
  618.                 {
  619.                     foundParam = parameterLibrary.findAtomParameter(receptor_atom_types[i]);
  620.                     if (foundParam != 0)
  621.                         foundParam->recIndex = i;
  622.                     else
  623.                     {
  624.                         logFile->printErrorFormatted(ERROR, "Unknown receptor type: "%s"n -- Add parameters for it to the parameter library first!n", receptor_atom_types[i]);
  625.                         logFile->printError(FATAL_ERROR, "Unsuccessful completion.nn");
  626.                     }
  627.                 }
  628.             }
  629.             break;
  631.         case GPF_SOL_PAR:      // THIS IS OBSOLETE!!!
  632.             // Read volume and solvation parameter for probe:
  633.             sscanf(GPFLine, "%*s %s %lf %lf", thisparm.autogridType, &temp_vol, &temp_solpar);
  634.             foundParam = parameterLibrary.findAtomParameter(thisparm.autogridType);
  635.             if (foundParam != 0)
  636.             {
  637.                 foundParam->vol = temp_vol;
  638.                 foundParam->solpar = temp_solpar;
  639.                 int i = foundParam->mapIndex;
  640.                 if (i >= 0)
  641.                 {
  642.                     // DON'T!!!
  643.                     // convert cal/molA^3 to kcal/molA^3
  644.                     // gridmaps[i].solparProbe = temp_solpar * 0.001;
  645.                     gridmaps[i].solparProbe = temp_solpar;
  646.                     logFile->printFormatted("nProbe %s solvation parameters: nntatomic fragmental volume: %.2f A^3ntatomic solvation parameter: %.4f cal/mol A^3nn",
  647.                                   foundParam->autogridType, foundParam->vol, foundParam->solpar);
  648.                 }
  649.             }
  650.             else
  651.                 logFile->printFormatted("%s key not foundn", thisparm.autogridType);
  652.             break;              // end solvation parameter
  654.         case GPF_MAP:
  655.             /* The variable "mapIndex" is the 0-based index of the ligand atom type we are calculating a map for. If the "types" line was CNOSH, there would be 5 ligand atom maps to
  656.                calculate, and since "mapIndex" is initialized to -1, mapIndex will increment each time there is a "map" keyword in the GPF.  The value of mapIndex should therefore go
  657.                from 0 to 4 for each "map" keyword. In this example, numAtomMaps would be 5, and numAtomMaps-1 would be 4, so if mapIndex is > 4, there is something wrong in the
  658.                number of "map" keywords. */
  659.             ++mapIndex;
  660.             if (mapIndex >= gridmaps.getNumAtomMaps())
  661.             {
  662.                 logFile->printErrorFormatted(ERROR,
  663.                     "Too many "map" keywords (%d);  the "ligand_types" command declares only %d atom types.nRemove a "map" keyword from the GPF.n",
  664.                     mapIndex + 1, gridmaps.getNumAtomMaps());
  665.                 logFile->printError(FATAL_ERROR, "Unsuccessful completion.nn");
  666.             }
  667.             // Read in the gridParameterFilename for this grid map
  668.             sscanf(GPFLine, "%*s %s", gridmaps[mapIndex].filename);
  669.             logFile->printFormatted("nOutput Grid Map %d:   %snn", (mapIndex + 1), gridmaps[mapIndex].filename);
  670.             break;
  672.         case GPF_ELECMAP:
  673.             sscanf(GPFLine, "%*s %s", gridmaps.getElectrostaticMap().filename);
  674.             logFile->printFormatted("nOutput Electrostatic Potential Energy Grid Map: %snn", gridmaps.getElectrostaticMap().filename);
  675.             break;
  677.         case GPF_DSOLVMAP:
  678.             sscanf(GPFLine, "%*s %s", gridmaps.getDesolvationMap().filename);
  679.             logFile->printFormatted("nOutput Desolvation Free Energy Grid Map: %snn", gridmaps.getDesolvationMap().filename);
  680.             break;
  682.         case GPF_COVALENTMAP:
  683.             {
  684.                 double covHalfWidth;
  685.                 sscanf(GPFLine, "%*s %lf %lf %lf %lf %lf", &covHalfWidth, &covBarrier, &(covalentPoint.x), &(covalentPoint.y), &(covalentPoint.z));
  686.                 covHalfWidthSquaredInv = 1 / Mathd::Sqr(covHalfWidth);
  688.                 logFile->printFormatted("ncovalentmap <half-width in Angstroms> <barrier> <x> <y> <z>n"
  689.                                        "nCovalent well's half-width in Angstroms:         %8.3fn"
  690.                                        "nCovalent barrier energy in kcal/mol:             %8.3fn"
  691.                                        "nCovalent attachment point will be positioned at: (%8.3f, %8.3f, %8.3f)nn",
  692.                                        covHalfWidth, covBarrier, covalentPoint.x, covalentPoint.y, covalentPoint.z);
  694.                 // center covalentPoint in the grid maps frame of reference,
  695.                 covalentPoint -= gridCenter;
  696.             }
  697.             break;
  699.         case GPF_DISORDER:
  700.             disorderH = true;
  701.             logFile->print("nHydroxyls will be disordered nn");
  702.             break;
  704.         case GPF_SMOOTH:
  705.             sscanf(GPFLine, "%*s %lf", &rSmooth);
  706.             logFile->printFormatted("nPotentials will be smoothed by: %.3lf Angstromnn", rSmooth);
  707.             break;
  709.         case GPF_QASP:
  710.             sscanf(GPFLine, "%*s %lf", &solparQ);
  711.             logFile->printFormatted("nCharge component of the atomic solvation parameter: %.3lfnn", solparQ);
  712.             // Typical value of solparQ is 0.001118
  713.             break;
  715.         case GPF_DIEL:
  716.             sscanf(GPFLine, "%*s %lf", &diel);
  717.             if (diel < 0)
  718.             {
  719.                 // negative...
  720.                 distDepDiel = true;
  721.                 // calculate inverted ddd of Mehler & Solmajer
  722.                 for (int indexR = 0; indexR < MAX_DIST; indexR++)
  723.                     epsilon[indexR] = calculateDistDepDielInv(indexToAngstrom<double>(indexR));
  724.                 logFile->print("nUsing *distance-dependent* dielectric function of Mehler and Solmajer, Prot.Eng.4, 903-910.nn"
  725.                               "  d   Dielectricn ___  __________n");
  726.                 for (int i = 0; i <= 500; i += 10)
  727.                 {
  728.                     ri = indexToAngstrom<double>(i);
  729.                     logFile->printFormatted("%4.1lf%9.2lfn", ri, DDD_FACTOR / epsilon[i]);
  730.                 }
  731.                 logFile->print("n");
  732.             }
  733.             else
  734.             {
  735.                 // positive or zero...
  736.                 distDepDiel = false;
  737.                 if (diel <= APPROX_ZERO)
  738.                     diel = 40;
  739.                 logFile->printFormatted("Using a *constant* dielectric of:  %.2fn", diel);
  740.                 invDielCal = DDD_FACTOR / diel;
  741.             }
  742.             break;
  744.         case GPF_FMAP:
  745.             sscanf(GPFLine, "%*s %s", floatingGridFilename);
  746.             logFile->printFormatted("nFloating Grid file name = %sn", floatingGridFilename);
  747.             break;
  749.         case GPF_PARAM_FILE:
  750.             // open and read the AD4 parameters .dat file
  751.             sscanf(GPFLine, "%*s %s ", parameterLibraryFilename);
  752.             break;
  753.         }                       // second switch
  754.     }                           // while: finished reading gpf
  756.     // Map files checkpoint  (number of maps, desolv and elec maps) SF
  758.     // Number of maps defined for atom types
  759.     if (mapIndex < gridmaps.getNumAtomMaps() - 1)
  760.     {
  761.         logFile->printFormatted("Too few "map" keywords (%d);  the "ligand_types" command declares %d atom types.nAdd a "map" keyword from the GPF.n", mapIndex + 1, gridmaps.getNumAtomMaps());
  762.         logFile->printError(ERROR, "Not enough map keywords found.n");
  763.         logFile->printError(FATAL_ERROR, "Unsuccessful completion.nn");
  764. }
  766.     // Desolvation map
  767.     if (!gridmaps.getDesolvationMap().filename[0])
  768.     {
  769.         logFile->print("The desolvation map file is not defined in the GPF.n");
  770.         logFile->printError(ERROR, "No desolvation map file defined.n");
  771.         logFile->printError(FATAL_ERROR, "Unsuccessful completion.nn");
  772. }
  774.     // Electrostatic map
  775.     if (!gridmaps.getElectrostaticMap().filename[0])
  776.     {
  777.         logFile->print("The electrostatic map file is not defined in the GPF.n");
  778.         logFile->printError(ERROR, "No electrostatic map file defined.n");
  779.         logFile->printError(FATAL_ERROR, "Unsuccessful completion.nn");
  780. }
  781.     // End of map files checkpoint SF
  783.     logFile->print("n>>> Closing the grid parameter file (GPF)... <<<nn" UnderLine);
  784.     fclose(GPF);
  786.     if (!floatingGridFilename[0])
  787.         logFile->print("nnNo Floating Grid was requested.n");
  789.     // apply the estat forcefield coefficient/weight here
  790.     if (distDepDiel)
  791.         for (int ia = 0; ia < numReceptorAtoms; ia++)
  792.             receptorAtom[ia].w = charge[ia] * parameterLibrary.coeff_estat;
  793.     else
  794.         // apply the constant dielectric
  795.         for (int ia = 0; ia < numReceptorAtoms; ia++)
  796.             receptorAtom[ia].w = charge[ia] * parameterLibrary.coeff_estat * invDielCal;
  797. }
  799. int InputDataLoader::parseGPFLine(const char *line)
  800. {
  801.     size_t l;
  802. int i, token = -1; // return -1 if nothing is recognized
  803.     char c[LINE_LEN];
  805.     l = strIndex(line, " ");
  806.     if (l == -1)
  807.     {
  808.         l = strIndex(line, "t");
  809.         if (l == -1)
  810.             l = strlen(line);
  811.     }
  812.     for(i=0; i<l; i++)
  813.         c[i] = char(tolower(line[i]));
  815.     if ((c[0]=='n')||(c[0]==''))
  816.         token = GPF_NULL;
  817.     else if (c[0]=='#')
  818.         token = GPF_COMMENT;
  819.     else if (strncmp(c, "receptor_types", 14) == 0)
  820.         token = GPF_RECEPTOR_TYPES;
  821.     else if (strncmp(c, "receptor", 8) == 0)
  822.         token = GPF_RECEPTOR;
  823.     else if (strncmp(c, "gridfld", 7) == 0)
  824.         token = GPF_GRIDFLD;
  825.     else if (strncmp(c, "npts", 4) == 0)
  826.         token = GPF_NPTS;
  827.     else if (strncmp(c, "spacing", 7) == 0)
  828.         token = GPF_SPACING;
  829.     else if (strncmp(c, "gridcenter", 10) == 0)
  830.         token = GPF_GRIDCENTER;
  831.     else if (strncmp(c, "types", 5) == 0)
  832.         token = GPF_LIGAND_TYPES;
  833.     else if (strncmp(c, "ligand_types", 12) == 0)
  834.         token = GPF_LIGAND_TYPES;
  835.     else if (strncmp(c, "map", 3) == 0)
  836.         token = GPF_MAP;
  837.     else if (strncmp(c, "elecmap", 7) == 0)
  838.         token = GPF_ELECMAP;
  839.     else if (strncmp(c, "dsolvmap", 8) == 0)
  840.         token = GPF_DSOLVMAP;
  841.     else if (strncmp(c, "covalentmap", 11) == 0)
  842.         token = GPF_COVALENTMAP;
  843.     else if (strncmp(c, "nbp_coeffs", 10) == 0)
  844.         token = GPF_NBP_COEFFS;
  845.     else if (strncmp(c, "nbp_r_eps", 9) == 0)
  846.         token = GPF_NBP_R_EPS;
  847.     else if (strncmp(c, "dielectric", 10) == 0)
  848.         token = GPF_DIEL;
  849.     else if (strncmp(c, "qasp", 4) == 0)
  850.         token = GPF_QASP;
  851.     else if (strncmp(c, "fmap", 4) == 0)
  852.         token = GPF_FMAP;
  853.     else if (strncmp(c, "disorder_h", 10) == 0)
  854.         token = GPF_DISORDER;
  855.     else if (strncmp(c, "smooth", 6) == 0)
  856.         token = GPF_SMOOTH;
  857.     else if (strncmp(c, "sol_par", 7) == 0)
  858.         token = GPF_SOL_PAR;
  859.     else if (strncmp(c, "constant", 8) == 0)
  860.         token = GPF_CONSTANT;
  861.     else if (strncmp(c, "parameter_file", 14) == 0)
  862.         token = GPF_PARAM_FILE;
  864.     return token;
  865. }
  867. // checks that number of grid elements is valid
  868. int InputDataLoader::checkSize(int numGridPointsMinusOne, char axischar)
  869. {
  870.     // numGridPointsMinusOne mustn't be negative, shouldn't be zero or larger than MAX_GRID_PTS and should be even
  871.     if (numGridPointsMinusOne < 0)
  872.         logFile->printErrorFormatted(FATAL_ERROR, "Negative number of %c-grid elements!  Aborting.nn", axischar);
  873.     else if (numGridPointsMinusOne == 0)
  874.         logFile->printErrorFormatted(WARNING, "0 %c-grid elements!nn", axischar);
  875.     else if (numGridPointsMinusOne > MAX_GRID_PTS)
  876.     {
  877.         logFile->printErrorFormatted(WARNING, "Maximum number of %c-grid elements allowed is %d. Using this value.n", axischar, MAX_GRID_PTS);
  878.         numGridPointsMinusOne = MAX_GRID_PTS;
  879.     }
  880.     else if (numGridPointsMinusOne % 2 == 1)
  881.     {
  882.         logFile->printTitledFormatted("Number of grid elements must be even; %c-elements changed to: %dn", axischar, numGridPointsMinusOne);
  883.         numGridPointsMinusOne -= 1;
  884.     }
  886.     return numGridPointsMinusOne;
  887. }
  889. // utility func for parsing types
  890. int InputDataLoader::parseTypes(char *line, char **words, int maxwords)
  891. {
  892.     char *char_ptr = line;
  893.     int num_types = 0;
  894.     // flag for first word which is always a keyword
  895.     int found_keyword = 0;
  896.     int index = 0;
  898.     for (;;)
  899.     {
  900.         // skip spaces
  901.         while (isspace(*char_ptr))
  902.         {
  903.             char_ptr++;
  904.             index++;
  905.         }
  907.         // done parsing when get eol 'null' character
  908.         // could get null after a space
  909.         if (*char_ptr == '')
  910.             return num_types; // return number of 'types' found
  912.         // the first word is the keyword not a type
  913.         if (!found_keyword)
  914.             found_keyword++;
  915.         else
  916.             words[num_types++] = char_ptr; // words is a list of indicies of beginning of 1 or 2 char types
  918.         // once in a type, skip possible 2nd characters up to a space or null character
  919.         while (!isspace(*char_ptr) && *char_ptr != '')
  920.         {
  921.             char_ptr++;
  922.             index++;
  923.         }
  925.         // done parsing when get eol 'null' character
  926.         // could get null after a character
  927.         if (*char_ptr == '')
  928.             return num_types;
  930.         // make each 'type' a null terminated string
  931.         *char_ptr++ = '';
  932.         index++;
  934.         //if there are too many types, return
  935.         if (num_types >= maxwords)
  936.             return num_types;
  937.     }
  938. }
  940. // returns index of t in s, -1 if none.
  941. int InputDataLoader::strIndex(const char *s, const char *t)
  942. {
  943.     const char *r = strstr(s, t);
  944.     return r? int(r-s) : -1;
  945. }