Program Listing for File PDBReader.h¶
↰ Return to documentation for file (include/gamer/PDBReader.h)
// This file is part of the GAMer software.
// Copyright (C) 2016-2021
// by Christopher T. Lee and contributors
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, see <http://www.gnu.org/licenses/>
// or write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330,
// Boston, MA 02111-1307 USA
/*
* Parts of this are adapted from the PDBParser developed in Chandrajit
* Bajaj's group at The University of Texas at Austin for GAMer by Zeyun Yu.
*/
#pragma once
#include <algorithm>
#include <array>
#include <fstream>
#include <iostream>
#include <regex>
#include <string>
#include "gamer/Vertex.h"
#include "gamer/gamer.h"
#include "gamer/SurfaceMesh.h"
namespace gamer {
namespace pdbreader_detail {
const double EPSILON = 1e-3;
static const std::regex PDB(".*.pdb", std::regex::icase | std::regex::optimize);
static const std::regex PQR(".*.pqr", std::regex::icase | std::regex::optimize);
static const std::regex XYZR(".*.xyzr",
std::regex::icase | std::regex::optimize);
static const std::regex atom("ATOM.*\n*", std::regex::optimize);
struct PDBelementInformation {
const char *atomName;
const char *residueName;
float radius;
float red;
float green;
float blue;
int hydrophobicity;
unsigned char residueIndex;
};
static const std::size_t MAX_BIOCHEM_ELEMENTS = 167;
static PDBelementInformation PDBelementTable[MAX_BIOCHEM_ELEMENTS] = {
{" N ", "GLY", 1.625f, 0.0f, 0.0f, 1.0f, 1, 10},
{" CA ", "GLY", 1.750f, 0.3f, 0.3f, 0.3f, -1, 10},
{" C ", "GLY", 1.875f, 0.3f, 0.3f, 0.3f, 1, 10},
{" O ", "GLY", 1.480f, 1.0f, 0.0f, 0.0f, 1, 10},
{" N ", "ALA", 1.625f, 0.0f, 0.0f, 1.0f, 1, 1},
{" CA ", "ALA", 1.750f, 0.3f, 0.3f, 0.3f, -1, 1},
{" C ", "ALA", 1.875f, 0.3f, 0.3f, 0.3f, 1, 1},
{" O ", "ALA", 1.480f, 1.0f, 0.0f, 0.0f, 1, 1},
{" CB ", "ALA", 1.750f, 0.3f, 0.3f, 0.3f, -1, 1},
{" N ", "VAL", 1.625f, 0.0f, 0.0f, 1.0f, 1, 23},
{" CA ", "VAL", 1.750f, 0.3f, 0.3f, 0.3f, -1, 23},
{" C ", "VAL", 1.875f, 0.3f, 0.3f, 0.3f, 1, 23},
{" O ", "VAL", 1.480f, 1.0f, 0.0f, 0.0f, 1, 23},
{" CB ", "VAL", 1.750f, 0.3f, 0.3f, 0.3f, -1, 23},
{" CG1", "VAL", 1.750f, 0.3f, 0.3f, 0.3f, -1, 23},
{" CG2", "VAL", 1.750f, 0.3f, 0.3f, 0.3f, -1, 23},
{" N ", "LEU", 1.625f, 0.0f, 0.0f, 1.0f, 1, 13},
{" CA ", "LEU", 1.750f, 0.3f, 0.3f, 0.3f, -1, 13},
{" C ", "LEU", 1.875f, 0.3f, 0.3f, 0.3f, 1, 13},
{" O ", "LEU", 1.480f, 1.0f, 0.0f, 0.0f, 1, 13},
{" CB ", "LEU", 1.750f, 0.3f, 0.3f, 0.3f, -1, 13},
{" CG ", "LEU", 1.750f, 0.3f, 0.3f, 0.3f, -1, 13},
{" CD1", "LEU", 1.750f, 0.3f, 0.3f, 0.3f, -1, 13},
{" CD2", "LEU", 1.750f, 0.3f, 0.3f, 0.3f, -1, 13},
{" N ", "ILE", 1.625f, 0.0f, 0.0f, 1.0f, 1, 12},
{" CA ", "ILE", 1.750f, 0.3f, 0.3f, 0.3f, -1, 12},
{" C ", "ILE", 1.875f, 0.3f, 0.3f, 0.3f, 1, 12},
{" O ", "ILE", 1.480f, 1.0f, 0.0f, 0.0f, 1, 12},
{" CB ", "ILE", 1.750f, 0.3f, 0.3f, 0.3f, -1, 12},
{" CG1", "ILE", 1.750f, 0.3f, 0.3f, 0.3f, -1, 12},
{" CG2", "ILE", 1.750f, 0.3f, 0.3f, 0.3f, -1, 12},
{" CD1", "ILE", 1.750f, 0.3f, 0.3f, 0.3f, -1, 12},
{" N ", "MET", 1.625f, 0.0f, 0.0f, 1.0f, 1, 15},
{" CA ", "MET", 1.750f, 0.3f, 0.3f, 0.3f, -1, 15},
{" C ", "MET", 1.875f, 0.3f, 0.3f, 0.3f, 1, 15},
{" O ", "MET", 1.480f, 1.0f, 0.0f, 0.0f, 1, 15},
{" CB ", "MET", 1.750f, 0.3f, 0.3f, 0.3f, -1, 15},
{" CG ", "MET", 1.750f, 0.3f, 0.3f, 0.3f, -1, 15},
{" SD ", "MET", 1.775f, 1.0f, 1.0f, 0.0f, 1, 15},
{" CE ", "MET", 1.750f, 0.3f, 0.3f, 0.3f, -1, 15},
{" N ", "PRO", 1.625f, 0.0f, 0.0f, 1.0f, -1, 17},
{" CA ", "PRO", 1.750f, 0.3f, 0.3f, 0.3f, -1, 17},
{" C ", "PRO", 1.875f, 0.3f, 0.3f, 0.3f, 1, 17},
{" O ", "PRO", 1.480f, 1.0f, 0.0f, 0.0f, 1, 17},
{" CB ", "PRO", 1.750f, 0.3f, 0.3f, 0.3f, -1, 17},
{" CG ", "PRO", 1.750f, 0.3f, 0.3f, 0.3f, -1, 17},
{" CD ", "PRO", 1.750f, 0.3f, 0.3f, 0.3f, -1, 17},
{" N ", "PHE", 1.625f, 0.0f, 0.0f, 1.0f, 1, 16},
{" CA ", "PHE", 1.750f, 0.3f, 0.3f, 0.3f, -1, 16},
{" C ", "PHE", 1.875f, 0.3f, 0.3f, 0.3f, 1, 16},
{" O ", "PHE", 1.480f, 1.0f, 0.0f, 0.0f, 1, 16},
{" CB ", "PHE", 1.750f, 0.3f, 0.3f, 0.3f, -1, 16},
{" CG ", "PHE", 1.775f, 0.3f, 0.3f, 0.3f, -1, 16},
{" CD1", "PHE", 1.775f, 0.3f, 0.3f, 0.3f, -1, 16},
{" CD2", "PHE", 1.775f, 0.3f, 0.3f, 0.3f, -1, 16},
{" CE1", "PHE", 1.775f, 0.3f, 0.3f, 0.3f, -1, 16},
{" CE2", "PHE", 1.775f, 0.3f, 0.3f, 0.3f, -1, 16},
{" CZ ", "PHE", 1.775f, 0.3f, 0.3f, 0.3f, -1, 16},
{" N ", "TRP", 1.625f, 0.0f, 0.0f, 1.0f, 1, 20},
{" CA ", "TRP", 1.750f, 0.3f, 0.3f, 0.3f, -1, 20},
{" C ", "TRP", 1.875f, 0.3f, 0.3f, 0.3f, 1, 20},
{" O ", "TRP", 1.480f, 1.0f, 0.0f, 0.0f, 1, 20},
{" CB ", "TRP", 1.750f, 0.3f, 0.3f, 0.3f, -1, 20},
{" CG ", "TRP", 1.775f, 0.3f, 0.3f, 0.3f, -1, 20},
{" CD1", "TRP", 1.775f, 0.3f, 0.3f, 0.3f, -1, 20},
{" CD2", "TRP", 1.775f, 0.3f, 0.3f, 0.3f, -1, 20},
{" NE1", "TRP", 1.625f, 0.2f, 0.2f, 1.0f, 1, 20},
{" CE2", "TRP", 1.775f, 0.3f, 0.3f, 0.3f, -1, 20},
{" CE3", "TRP", 1.775f, 0.3f, 0.3f, 0.3f, -1, 20},
{" CZ2", "TRP", 1.775f, 0.3f, 0.3f, 0.3f, -1, 20},
{" CZ3", "TRP", 1.775f, 0.3f, 0.3f, 0.3f, -1, 20},
{" CH2", "TRP", 1.775f, 0.3f, 0.3f, 0.3f, -1, 20},
{" N ", "SER", 1.625f, 0.0f, 0.0f, 1.0f, 1, 18},
{" CA ", "SER", 1.750f, 0.3f, 0.3f, 0.3f, -1, 18},
{" C ", "SER", 1.875f, 0.3f, 0.3f, 0.3f, 1, 18},
{" O ", "SER", 1.480f, 1.0f, 0.0f, 0.0f, 1, 18},
{" CB ", "SER", 1.750f, 0.3f, 0.3f, 0.3f, -1, 18},
{" OG ", "SER", 1.560f, 1.0f, 0.0f, 0.0f, 1, 18},
{" N ", "THR", 1.625f, 0.0f, 0.0f, 1.0f, 1, 19},
{" CA ", "THR", 1.750f, 0.3f, 0.3f, 0.3f, -1, 19},
{" C ", "THR", 1.875f, 0.3f, 0.3f, 0.3f, 1, 19},
{" O ", "THR", 1.480f, 1.0f, 0.0f, 0.0f, 1, 19},
{" CB ", "THR", 1.750f, 0.3f, 0.3f, 0.3f, -1, 19},
{" OG1", "THR", 1.560f, 1.0f, 0.0f, 0.0f, 1, 19},
{" CG2", "THR", 1.750f, 0.3f, 0.3f, 0.3f, -1, 19},
{" N ", "ASN", 1.625f, 0.0f, 0.0f, 1.0f, 1, 3},
{" CA ", "ASN", 1.750f, 0.3f, 0.3f, 0.3f, -1, 3},
{" C ", "ASN", 1.875f, 0.3f, 0.3f, 0.3f, 1, 3},
{" O ", "ASN", 1.480f, 1.0f, 0.0f, 0.0f, 1, 3},
{" CB ", "ASN", 1.750f, 0.3f, 0.3f, 0.3f, -1, 3},
{" CG ", "ASN", 1.875f, 0.3f, 0.3f, 0.3f, 1, 3},
{" OD1", "ASN", 1.480f, 1.0f, 0.0f, 0.0f, 1, 3},
{" ND2", "ASN", 1.625f, 0.2f, 0.2f, 1.0f, 1, 3},
{" N ", "GLN", 1.625f, 0.0f, 0.0f, 1.0f, 1, 7},
{" CA ", "GLN", 1.750f, 0.3f, 0.3f, 0.3f, -1, 7},
{" C ", "GLN", 1.875f, 0.3f, 0.3f, 0.3f, 1, 7},
{" O ", "GLN", 1.480f, 1.0f, 0.0f, 0.0f, 1, 7},
{" CB ", "GLN", 1.750f, 0.3f, 0.3f, 0.3f, -1, 7},
{" CG ", "GLN", 1.750f, 0.3f, 0.3f, 0.3f, -1, 7},
{" CD ", "GLN", 1.875f, 0.3f, 0.3f, 0.3f, 1, 7},
{" OE1", "GLN", 1.480f, 1.0f, 0.0f, 0.0f, 1, 7},
{" NE2", "GLN", 1.625f, 0.2f, 0.2f, 1.0f, 1, 7},
{" N ", "TYR", 1.625f, 0.0f, 0.0f, 1.0f, 1, 21},
{" CA ", "TYR", 1.750f, 0.3f, 0.3f, 0.3f, -1, 21},
{" C ", "TYR", 1.875f, 0.3f, 0.3f, 0.3f, 1, 21},
{" O ", "TYR", 1.480f, 1.0f, 0.0f, 0.0f, 1, 21},
{" CB ", "TYR", 1.750f, 0.3f, 0.3f, 0.3f, -1, 21},
{" CG ", "TYR", 1.775f, 0.3f, 0.3f, 0.3f, -1, 21},
{" CD1", "TYR", 1.775f, 0.3f, 0.3f, 0.3f, -1, 21},
{" CD2", "TYR", 1.775f, 0.3f, 0.3f, 0.3f, -1, 21},
{" CE1", "TYR", 1.775f, 0.3f, 0.3f, 0.3f, -1, 21},
{" CE2", "TYR", 1.775f, 0.3f, 0.3f, 0.3f, -1, 21},
{" CZ ", "TYR", 1.775f, 0.3f, 0.3f, 0.3f, -1, 21},
{" OH ", "TYR", 1.535f, 1.0f, 0.0f, 0.0f, 1, 21},
{" N ", "CYS", 1.625f, 0.0f, 0.0f, 1.0f, 1, 6},
{" CA ", "CYS", 1.750f, 0.3f, 0.3f, 0.3f, -1, 6},
{" C ", "CYS", 1.875f, 0.3f, 0.3f, 0.3f, 1, 6},
{" O ", "CYS", 1.480f, 1.0f, 0.0f, 0.0f, 1, 6},
{" CB ", "CYS", 1.750f, 0.3f, 0.3f, 0.3f, -1, 6},
{" SG ", "CYS", 1.775f, 1.0f, 1.0f, 0.0f, 1, 6},
{" N ", "LYS", 1.625f, 0.0f, 0.0f, 1.0f, 1, 14},
{" CA ", "LYS", 1.750f, 0.3f, 0.3f, 0.3f, -1, 14},
{" C ", "LYS", 1.875f, 0.3f, 0.3f, 0.3f, 1, 14},
{" O ", "LYS", 1.480f, 1.0f, 0.0f, 0.0f, 1, 14},
{" CB ", "LYS", 1.750f, 0.3f, 0.3f, 0.3f, -1, 14},
{" CG ", "LYS", 1.750f, 0.3f, 0.3f, 0.3f, -1, 14},
{" CD ", "LYS", 1.750f, 0.3f, 0.3f, 0.3f, -1, 14},
{" CE ", "LYS", 1.750f, 0.3f, 0.3f, 0.3f, -1, 14},
{" NZ ", "LYS", 1.625f, 0.2f, 0.2f, 1.0f, 1, 14},
{" N ", "ARG", 1.625f, 0.0f, 0.0f, 1.0f, 1, 2},
{" CA ", "ARG", 1.750f, 0.3f, 0.3f, 0.3f, -1, 2},
{" C ", "ARG", 1.875f, 0.3f, 0.3f, 0.3f, 1, 2},
{" O ", "ARG", 1.480f, 1.0f, 0.0f, 0.0f, 1, 2},
{" CB ", "ARG", 1.750f, 0.3f, 0.3f, 0.3f, -1, 2},
{" CG ", "ARG", 1.750f, 0.3f, 0.3f, 0.3f, -1, 2},
{" CD ", "ARG", 1.750f, 0.3f, 0.3f, 0.3f, -1, 2},
{" NE ", "ARG", 1.625f, 0.2f, 0.2f, 1.0f, 1, 2},
{" CZ ", "ARG", 1.125f, 0.3f, 0.3f, 0.3f, 1, 2},
{" NH1", "ARG", 1.625f, 0.2f, 0.2f, 1.0f, 1, 2},
{" NH2", "ARG", 1.625f, 0.2f, 0.2f, 1.0f, 1, 2},
{" N ", "HIS", 1.625f, 0.0f, 0.0f, 1.0f, 1, 11},
{" CA ", "HIS", 1.750f, 0.3f, 0.3f, 0.3f, -1, 11},
{" C ", "HIS", 1.875f, 0.3f, 0.3f, 0.3f, 1, 11},
{" O ", "HIS", 1.480f, 1.0f, 0.0f, 0.0f, 1, 11},
{" CB ", "HIS", 1.750f, 0.3f, 0.3f, 0.3f, -1, 11},
{" CG ", "HIS", 1.775f, 0.3f, 0.3f, 0.3f, -1, 11},
{" ND1", "HIS", 1.625f, 0.2f, 0.2f, 1.0f, 1, 11},
{" CD2", "HIS", 1.775f, 0.3f, 0.3f, 0.3f, -1, 11},
{" CE1", "HIS", 1.775f, 0.3f, 0.3f, 0.3f, 1, 11},
{" NE2", "HIS", 1.625f, 0.2f, 0.2f, 1.0f, 1, 11},
{" N ", "ASP", 1.625f, 0.0f, 0.0f, 1.0f, 1, 4},
{" CA ", "ASP", 1.750f, 0.3f, 0.3f, 0.3f, -1, 4},
{" C ", "ASP", 1.875f, 0.3f, 0.3f, 0.3f, 1, 4},
{" O ", "ASP", 1.480f, 1.0f, 0.0f, 0.0f, 1, 4},
{" CB ", "ASP", 1.750f, 0.3f, 0.3f, 0.3f, -1, 4},
{" CG ", "ASP", 1.875f, 0.3f, 0.3f, 0.3f, 1, 4},
{" OD1", "ASP", 1.480f, 1.0f, 1.0f, 1.0f, 1, 4},
{" OD2", "ASP", 1.480f, 1.0f, 0.0f, 0.0f, 1, 4},
{" N ", "GLU", 1.625f, 0.0f, 0.0f, 1.0f, 1, 8},
{" CA ", "GLU", 1.750f, 0.3f, 0.3f, 0.3f, -1, 8},
{" C ", "GLU", 1.875f, 0.3f, 0.3f, 0.3f, 1, 8},
{" O ", "GLU", 1.480f, 1.0f, 0.0f, 0.0f, 1, 8},
{" CB ", "GLU", 1.750f, 0.3f, 0.3f, 0.3f, -1, 8},
{" CG ", "GLU", 1.750f, 0.3f, 0.3f, 0.3f, -1, 8},
{" CD ", "GLU", 1.875f, 0.3f, 0.3f, 0.3f, 1, 8},
{" OE1", "GLU", 1.480f, 1.0f, 0.0f, 0.0f, 1, 8},
{" OE2", "GLU", 1.480f, 1.0f, 0.0f, 0.0f, 1, 8}
// {"SI ", "UNL", 1.875f, 1.0f, 1.0f, 1.0f, 1, 27 },
// {" O ", "UNL", 1.480f, 1.0f, 0.0f, 0.0f, 1, 27 }
};
static std::map<std::string, std::map<std::string, PDBelementInformation>>
PDBelementMap;
} // End namespace pdbreader_detail
struct Atom {
Vector3f pos;
double radius;
};
static void initElementMap() {
// If the element map is empty build it...
if (pdbreader_detail::PDBelementMap.size() == 0) {
for (int i = 0; i < pdbreader_detail::MAX_BIOCHEM_ELEMENTS; ++i) {
const std::string atomName =
pdbreader_detail::PDBelementTable[i].atomName;
const std::string residueName =
pdbreader_detail::PDBelementTable[i].residueName;
pdbreader_detail::PDBelementMap[residueName][atomName] =
pdbreader_detail::PDBelementTable[i];
}
}
}
template <typename Inserter>
bool readPDB(const std::string &filename, Inserter inserter) {
std::ifstream infile(filename);
std::string line;
initElementMap(); // Init the element map if it
if (infile.is_open()) {
while (std::getline(infile, line)) {
std::smatch match;
if (std::regex_match(line, match, pdbreader_detail::atom)) {
Atom atom;
// See PDB file formatting guidelines
float x = std::atof(line.substr(30, 8).c_str());
float y = std::atof(line.substr(38, 8).c_str());
float z = std::atof(line.substr(46, 8).c_str());
atom.pos = Vector({x, y, z});
atom.radius = 1.0f; // default radius
std::string atomName = line.substr(12, 4);
std::string residueName = line.substr(17, 3);
auto innerMapIT = pdbreader_detail::PDBelementMap.find(residueName);
if (innerMapIT != pdbreader_detail::PDBelementMap.end()) {
auto innerMap = innerMapIT->second;
auto typeIT = innerMap.find(atomName);
if (typeIT != innerMap.end()) {
atom.radius = typeIT->second.radius;
} else {
std::cout << "Could not find atomtype of '" << atomName
<< "' in residue '" << residueName << "'. "
<< "Using default radius." << std::endl;
}
} else {
std::cout << "Could not find ResidueName '" << residueName
<< "' in table. "
<< "Using default radius." << std::endl;
}
*inserter++ = atom;
}
}
return true;
} else {
std::cerr << "Unable to open \"" << filename << "\"" << std::endl;
return false;
}
}
template <typename Inserter>
bool readPQR(const std::string &filename, Inserter inserter) {
std::ifstream infile(filename);
std::string line;
initElementMap(); // Init the element map if it
if (infile.is_open()) {
while (std::getline(infile, line)) {
std::smatch match;
if (std::regex_match(line, match, pdbreader_detail::atom)) {
Atom atom;
// See PDB file formatting guidelines
float x = std::atof(line.substr(30, 8).c_str());
float y = std::atof(line.substr(38, 8).c_str());
float z = std::atof(line.substr(46, 8).c_str());
atom.pos = Vector({x, y, z});
atom.radius = std::atof(line.substr(62, 7).c_str());
*inserter++ = atom;
}
}
return true;
} else {
std::cerr << "Unable to open \"" << filename << "\"" << std::endl;
return false;
}
}
template <typename Iterator, typename BlurFunc>
void getMinMax(Iterator begin, Iterator end, Vector3f &min, Vector3f &max,
BlurFunc &&f) {
float maxRad = 0.0;
float tmpRad;
min[0] = min[1] = min[2] = std::numeric_limits<float>::infinity();
max[0] = max[1] = max[2] = -std::numeric_limits<float>::infinity();
for (auto curr = begin; curr != end; ++curr) {
float x = curr->pos[0];
float y = curr->pos[1];
float z = curr->pos[2];
if (min[0] > x)
min[0] = x;
if (max[0] < x)
max[0] = x;
if (min[1] > y)
min[1] = y;
if (max[1] < y)
max[1] = y;
if (min[2] > z)
min[2] = z;
if (max[2] < z)
max[2] = z;
tmpRad = f(curr->radius); // * sqrt(1.0 + log(pdbreader_detail::EPSILON)
// / blobbyness);
if (maxRad < tmpRad)
maxRad = tmpRad;
}
min -= Vector3f({maxRad, maxRad, maxRad});
max += Vector3f({maxRad, maxRad, maxRad});
}
template <typename Iterator>
void blurAtoms(Iterator begin, Iterator end, float *dataset,
const Vector3f &min, const Vector3f &maxMin, const Vector3i &dim,
float blobbyness) {
// Functor to calculate gaussian blur
auto evalDensity = [blobbyness](const Atom &atom, Vector3f &pnt,
float maxRadius) -> float {
double expval;
Vector3f tmp = atom.pos - pnt;
double r = tmp | tmp;
double r0 = atom.radius * atom.radius;
// expval = BLOBBYNESS*(r/r0 - 1.0);
expval = blobbyness * (r - r0);
// Truncate gaussian
if (sqrt(r) > maxRadius) {
return 0.0;
}
return (float)exp(expval);
};
Vector3f span;
span = (maxMin).ElementwiseDivision(
static_cast<Vector3f>((dim - Vector3i({1, 1, 1}))));
float radFactor =
sqrt(1.0 + log(pdbreader_detail::EPSILON) / (2.0 * blobbyness));
for (auto curr = begin; curr != end; ++curr) {
float maxRad = curr->radius * radFactor;
// compute the dataset coordinates of the atom's center
Vector3f tmpVec = (curr->pos - min).ElementwiseDivision(span);
Vector3i c;
std::transform(tmpVec.begin(), tmpVec.end(), c.begin(),
[](float v) -> int { return round(v); });
// std::cout << "Max Radius: " << maxRad << std::endl;
// compute the bounding box of the atom (maxRad^3)
Vector3i amin;
Vector3i amax;
for (int j = 0; j < 3; ++j) {
int tmp;
float tmpRad = maxRad / span[j];
tmp = (int)(c[j] - tmpRad - 1);
amin[j] = (tmp < 0) ? 0 : tmp; // check if tmp is < 0
tmp = (int)(c[j] + tmpRad + 1);
amax[j] = (tmp > (dim[j] - 1)) ? (dim[j] - 1) : tmp;
}
// std::cout << amin << " " << amax << std::endl;
// Blur kernel in bounding box
for (int k = amin[2]; k <= amax[2]; k++) {
for (int j = amin[1]; j <= amax[1]; j++) {
for (int i = amin[0]; i <= amax[0]; i++) {
Vector3f pnt =
min + Vector3f({static_cast<float>(i), static_cast<float>(j),
static_cast<float>(k)})
.ElementwiseProduct(span);
dataset[Vect2Index(i, j, k, dim)] += evalDensity(*curr, pnt, maxRad);
}
}
}
}
}
template <typename Iterator>
void gridSAS(const Iterator begin, const Iterator end, const Vector3i &dim,
float *dataset) {
// For atom in atoms :
for (auto curr = begin; curr != end; ++curr) {
float radius = curr->radius;
Vector3f pos = curr->pos;
// compute the dataset coordinates of the atom's center
Vector3i c;
std::transform(pos.begin(), pos.end(), c.begin(),
[](float v) -> int { return round(v); });
// compute bounding box for atom
Vector3i amin;
Vector3i amax;
for (int j = 0; j < 3; ++j) {
int tmp;
tmp = (int)(c[j] - radius - 1);
amin[j] = (tmp < 0) ? 0 : tmp; // check if tmp is < 0
tmp = (int)(c[j] + radius + 1);
amax[j] = (tmp > (dim[j] - 1)) ? (dim[j] - 1) : tmp;
}
// Blur kernel in bounding box
for (int k = amin[2]; k <= amax[2]; k++) {
for (int j = amin[1]; j <= amax[1]; j++) {
for (int i = amin[0]; i <= amax[0]; i++) {
Vector3f coord =
Vector3f({static_cast<float>(i), static_cast<float>(j),
static_cast<float>(k)});
coord -= curr->pos;
float dist = -(std::sqrt(coord | coord) - radius); // inside
// is
// positive
int idx = Vect2Index(i, j, k, dim);
if (dist > dataset[idx]) {
dataset[idx] = dist;
}
}
}
}
}
}
template <typename Iterator>
void gridSES(const Iterator begin, const Iterator end, const Vector3i &dim,
float *dataset, const float radius) {
for (auto curr = begin; curr != end; ++curr) {
Vector3f pos = (*curr).position;
// compute bounding box for atom
Vector3i amin;
Vector3i amax;
for (int i = 0; i < 3; ++i) {
int tmp;
tmp = (int)(pos[i] - radius - 1);
amin[i] = (tmp < 0) ? 0 : tmp; // check if tmp is < 0
tmp = (int)(pos[i] + radius + 1);
amax[i] = (tmp > (dim[i] - 1)) ? (dim[i] - 1) : tmp;
}
// Blur kernel in bounding box
for (int k = amin[2]; k <= amax[2]; k++) {
for (int j = amin[1]; j <= amax[1]; j++) {
for (int i = amin[0]; i <= amax[0]; i++) {
Vector3f coord =
Vector3f({static_cast<float>(i), static_cast<float>(j),
static_cast<float>(k)});
coord -= pos;
float dist = -(std::sqrt(coord | coord) - radius);
int idx = Vect2Index(i, j, k, dim);
if (dist > dataset[idx]) {
dataset[idx] = dist;
}
}
}
}
}
}
// TODO: (1) these functions should all take arrays of x,y,z,r instead of
// filenames
std::unique_ptr<SurfaceMesh> readPDB_molsurf(const std::string &filename);
std::unique_ptr<SurfaceMesh> readPDB_gauss(const std::string &filename,
float blobbyness, float isovalue);
std::unique_ptr<SurfaceMesh> readPDB_distgrid(const std::string &filename,
const float radius);
std::unique_ptr<SurfaceMesh> readPQR_molsurf(const std::string &filename);
std::unique_ptr<SurfaceMesh> readPQR_gauss(const std::string &filename,
float blobbyness, float isovalue);
} // end namespace gamer