From ebf924d914ddca929ebf138518eba5994ff954d7 Mon Sep 17 00:00:00 2001
From: Gabriel Studer <gabriel.studer@unibas.ch>
Date: Thu, 15 Jun 2023 14:23:33 +0200
Subject: [PATCH] read charges in PDB reader if present
---
modules/io/src/mol/pdb_reader.cc | 34 ++++++++++++++++++++++++--------
1 file changed, 26 insertions(+), 8 deletions(-)
diff --git a/modules/io/src/mol/pdb_reader.cc b/modules/io/src/mol/pdb_reader.cc
index 3b69915de..930ffbcf4 100644
--- a/modules/io/src/mol/pdb_reader.cc
+++ b/modules/io/src/mol/pdb_reader.cc
@@ -697,8 +697,11 @@ void PDBReader::ParseAndAddAtom(const StringRef& line, int line_num,
atom_name, alt_loc, record_type, serial)) {
return;
}
- std::pair<bool, Real> charge, radius;
- std::pair<bool, Real> occ, temp;
+
+ std::pair<bool, Real> charge = std::make_pair(false, 0.0);
+ std::pair<bool, Real> radius = std::make_pair(false, 0.0);
+ std::pair<bool, Real> occ = std::make_pair(false, 0.0);
+ std::pair<bool, Real> temp = std::make_pair(false, 0.0);
geom::Vec3 apos;
for (int i=0;i<3; ++i) {
@@ -728,6 +731,20 @@ void PDBReader::ParseAndAddAtom(const StringRef& line, int line_num,
if (line.length()>=66) {
temp=line.substr(60, 6).ltrim().to_float();
}
+ if (line.length()>=80) {
+ // example charge formatting: 1+, 2- etc.
+ charge = line.substr(78,1).trim().to_float();
+ if(charge.first) {
+ if(line[79] != '-' && line[79] != '+') {
+ std::stringstream ss;
+ ss << "error on line " << line_num << ": "
+ << "expect charge in format 1+, 2-, etc. got: "
+ << line.substr(78, 2);
+ throw IOException(ss.str());
+ }
+ if(line[79] == '-') charge.second *= (-1);
+ }
+ }
}
LOG_TRACE( "line: [" << line << "]" );
String s_ele;
@@ -838,6 +855,7 @@ void PDBReader::ParseAndAddAtom(const StringRef& line, int line_num,
}
Real b=temp.first ? temp.second : 0.0;
Real o=occ.first ? occ.second : 1.0;
+
if (!profile_.quack_mode && alt_loc!=' ') {
// Check if there is already a atom with the same name.
mol::AtomHandle me=curr_residue_.FindAtom(aname);
@@ -869,16 +887,16 @@ void PDBReader::ParseAndAddAtom(const StringRef& line, int line_num,
ah=editor.InsertAtom(curr_residue_, aname, apos, s_ele);
++atom_count_;
}
- if(is_pqr_) {
- if (radius.first) {
- ah.SetRadius(radius.second);
- }
+
+ if (radius.first) {
+ ah.SetRadius(radius.second);
}
- ah.SetBFactor(b);
- ah.SetOccupancy(o);
if (charge.first) {
ah.SetCharge(charge.second);
}
+ ah.SetBFactor(b);
+ ah.SetOccupancy(o);
+
ah.SetHetAtom(record_type[0]=='H');
if(profile_.read_conect && serial != -1) {
this->amap_[serial] = ah;
--
GitLab