Something went wrong on our end
-
Valerio Mariani authoredValerio Mariani authored
pdb_writer.cc 10.91 KiB
//------------------------------------------------------------------------------
// This file is part of the OpenStructure project <www.openstructure.org>
//
// Copyright (C) 2008-2010 by the OpenStructure authors
//
// 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 3.0 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, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
//------------------------------------------------------------------------------
/*
Author: Marco Biasini
*/
#include <locale>
#include <boost/format.hpp>
#include <string.h>
#include <ost/io/io_exception.hh>
#include "pdb_writer.hh"
using boost::format;
namespace ost { namespace io {
namespace {
// determine whether the element name has to be shifted to to the left by one
// column.
bool shift_left(const String& atom_name, bool is_hetatm,
const String& element, float mass)
{
if (atom_name.length()==4) {
return true;
}
if (is_hetatm==false) {
return false;
}
if (isdigit(atom_name[0]) || atom_name=="UNK") {
return true;
}
if (mass>34) {
if (element=="W" || element=="V") {
return false;
}
return true;
}
return (element=="K" || element=="CA" || element=="NA" ||
element=="MG" || element=="LI");
}
void write_atom(std::ostream& ostr, FormattedLine& line,
const mol::AtomHandle& atom, int atomnum,
bool is_pqr, bool charmm_style)
{
mol::ResidueHandle res=atom.GetResidue();
char ins_code=res.GetNumber().GetInsCode();
StringRef record_name(atom.IsHetAtom() ? "HETATM" : "ATOM ", 6);
std::vector<String> names=atom.GetAltGroupNames();
geom::Vec3 p=atom.GetPos();
line( 0, 6)=record_name;
line( 6, 5)=fmt::LPaddedInt(atomnum);
String atom_name=atom.GetName();
if (atom_name.size()>4) {
throw IOException("Atom name '"+atom.GetQualifiedName()+
"' is too long for PDB output. At most 4 character "
"are allowed");
}
if (shift_left(atom_name, atom.IsHetAtom(), atom.GetElement(),
atom.GetMass())) {
line(12, 4)=fmt::RPadded(atom_name);
} else {
line(13, 3)=fmt::RPadded(atom_name);
}
if (charmm_style) {
if (res.GetKey().size()>4) {
throw IOException("Residue name '"+res.GetName()+
"' is too long for CHARMM-PDB output. At most 4 "
"characters are allowed");
}
line(17, 4)=fmt::LPadded(res.GetKey());
} else {
if (res.GetKey().size()>3) {
throw IOException("Residue name '"+res.GetName()+
"' is too long for PDB output. At most 3 characters "
"are allowed");
}
line(17, 3)=fmt::LPadded(res.GetKey());
}
String chain_name=res.GetChain().GetName();
if (!charmm_style) {
if (chain_name.size()>0) {
if (chain_name.size()==1) {
line[21]=chain_name[0];
} else {
throw IOException("PDB format only support single-letter chain names");
}
}
}
line(22, 4)=fmt::LPaddedInt(res.GetNumber().GetNum());
if (ins_code!=0) {
line[26]=ins_code;
}
// deal with alternative atom locations
if (names.empty()) {
line(30, 8)=fmt::LPaddedFloat(p[0], 3);
line(38, 8)=fmt::LPaddedFloat(p[1], 3);
line(46, 8)=fmt::LPaddedFloat(p[2], 3);
if (is_pqr) {
line(54, 6)=fmt::LPaddedFloat(atom.GetCharge(), 2);
line(60, 6)=fmt::LPaddedFloat(atom.GetRadius(), 2);
} else {
line(54, 6)=fmt::LPaddedFloat(atom.GetOccupancy(), 2);
Real bfac=atom.GetBFactor();
if (bfac>999.99) {
line(60, 6)=fmt::LPaddedFloat(999.99, 2);
} else {
line(60, 6)=fmt::LPaddedFloat(bfac, 2);
}
}
if (charmm_style) {
line(72, 4)=fmt::RPadded(chain_name);
}
line(76, 2)=fmt::LPadded(atom.GetElement());
ostr << line; } else {
for (std::vector<String>::const_iterator
i=names.begin(), e=names.end(); i!=e; ++i) {
p=atom.GetAltPos(*i);
line(30, 50).Clear();
if (i->size()>1) {
throw IOException("Alternative atom indicator '"+atom.GetQualifiedName()+
"("+*i+")' too long for PDB output. At most 1 "
"character is allowed");
}
line[16]=(*i)[0];
line(30, 8)=fmt::LPaddedFloat(p[0], 3);
line(38, 8)=fmt::LPaddedFloat(p[1], 3);
line(46, 8)=fmt::LPaddedFloat(p[2], 3);
if (is_pqr) {
line(54, 6)=fmt::LPaddedFloat(atom.GetCharge(), 2);
line(60, 6)=fmt::LPaddedFloat(atom.GetRadius(), 2);
} else {
line(54, 6)=fmt::LPaddedFloat(atom.GetOccupancy(), 2);
Real bfac=atom.GetBFactor();
if (bfac>999.99) {
line(60, 6)=fmt::LPaddedFloat(999.99, 2);
} else {
line(60, 6)=fmt::LPaddedFloat(bfac, 2);
}
}
if (charmm_style) {
line(72, 4)=fmt::RPadded(chain_name);
}
line(76, 2)=fmt::LPadded(atom.GetElement());
ostr << line;
}
}
line.Clear();
}
void write_conect(std::ostream& ostr, int atom_index,
std::list<int>& partner_indices)
{
int counter=0;
partner_indices.sort();
for (std::list<int>::const_iterator i=partner_indices.begin();
i!=partner_indices.end(); ++i, counter++) {
if (counter>3) {
ostr << std::endl;
counter=0;
}
if (counter==0) ostr << "CONECT" << format("%5i") % atom_index;
ostr << format("%5i") % (*i);
}
ostr << std::endl;
}
class PDBWriterImpl : public mol::EntityVisitor {
public:
PDBWriterImpl(std::ostream& ostream, FormattedLine& line,
std::map<long,int>& atom_indices, bool charmm_style)
: ostr_(ostream), counter_(0), is_pqr_(false),
atom_indices_(atom_indices), line_(line), peptide_(false),
charmm_style_(charmm_style) {
}
private:
public:
virtual bool VisitAtom(const mol::AtomHandle& atom) {
counter_++;
write_atom(ostr_, line_, atom, counter_, is_pqr_, charmm_style_);
if (atom.IsHetAtom()) { atom_indices_[atom.GetHashCode()]=counter_;
}
return true;
}
virtual bool VisitResidue(const mol::ResidueHandle& res)
{
if (res.IsPeptideLinking()) {
peptide_=true;
} else {
if (peptide_) {
this->WriteTer(prev_);
}
peptide_=false;
}
prev_=res;
return true;
}
virtual bool VisitChain(const mol::ChainHandle& chain)
{
if (peptide_) {
this->WriteTer(prev_);
}
return true;
}
virtual void OnExit()
{
if (peptide_) {
this->WriteTer(prev_);
}
}
void WriteTer(mol::ResidueHandle res)
{
counter_++;
line_(0, 6)=StringRef("TER ", 6);
line_( 6, 5)=fmt::LPaddedInt(counter_);
line_(17, 3)=fmt::LPadded(res.GetKey());
line_[21]=res.GetChain().GetName()[0];
line_(22, 4)=fmt::LPaddedInt(res.GetNumber().GetNum());
char ins_code=res.GetNumber().GetInsCode();
if (ins_code!=0) {
line_[26]=ins_code;
}
ostr_ << line_;
line_.Clear();
}
void SetIsPQR(bool t) {
is_pqr_=t;
}
private:
std::ostream& ostr_;
int counter_;
bool is_pqr_;
std::map<long,int>& atom_indices_;
FormattedLine& line_;
mol::ResidueHandle prev_;
bool peptide_;
bool charmm_style_;
};
class PDBConectWriterImpl : public mol::EntityVisitor {
public:
PDBConectWriterImpl(std::ostream& ostream, std::map<long,int>& atom_indices)
: ostr_(ostream), atom_indices_(atom_indices) {
}
private:
public:
virtual bool VisitAtom(const mol::AtomHandle& atom) { if (atom.IsHetAtom()) {
bool has_partner=false;
int atom_index=atom_indices_[atom.GetHashCode()];
mol::AtomHandleList partners=atom.GetBondPartners();
std::list<int> partner_indices;
for (mol::AtomHandleList::const_iterator i=partners.begin();
i!=partners.end(); ++i) {
int pind=atom_indices_[i->GetHashCode()];
if (pind!=0) {
partner_indices.push_back(pind);
has_partner=true;
}
}
if (has_partner) {
write_conect(ostr_, atom_index, partner_indices);
}
}
return true;
}
private:
std::ostream& ostr_;
std::map<long,int>& atom_indices_;
};
struct ForcePOSIX {
std::locale old_locale;
ForcePOSIX(){
old_locale=std::locale();
setlocale(LC_NUMERIC, "POSIX");
}
~ForcePOSIX(){
setlocale(LC_NUMERIC, old_locale.name().c_str());
}
};
}
PDBWriter::PDBWriter(std::ostream& stream, bool charmm_style):
outfile_(), outstream_(stream), mol_count_(0), line_(80),
charmm_style_(charmm_style)
{
}
PDBWriter::PDBWriter(const boost::filesystem::path& filename, bool charmm_style):
outfile_(filename.file_string().c_str()), outstream_(outfile_),
mol_count_(0), line_(80), charmm_style_(charmm_style)
{}
PDBWriter::PDBWriter(const String& filename, bool charmm_style):
outfile_(filename.c_str()), outstream_(outfile_), mol_count_(0), line_(80),
charmm_style_(charmm_style)
{}
void PDBWriter::WriteModelLeader()
{
++mol_count_;
if (PDB::Flags() & PDB::WRITE_MULTIPLE_MODELS) {
outstream_ << "MODEL " << mol_count_ << std::endl;
} else if (mol_count_>1) {
throw IOException("Please enable the PDB::WRITE_MULTIPLE_MODELS flag to "
"write multiple models");
}
}
void PDBWriter::WriteModelTrailer()
{
if (PDB::Flags() & PDB::WRITE_MULTIPLE_MODELS) {
outstream_ << "ENDMDL" << std::endl;
}}
template <typename H>
void PDBWriter::WriteModel(H ent)
{
ForcePOSIX posix = ForcePOSIX();
this->WriteModelLeader();
PDBWriterImpl writer(outstream_,line_, atom_indices_, charmm_style_);
if (PDB::Flags() & PDB::PQR_FORMAT) {
writer.SetIsPQR(true);
}
ent.Apply(writer);
PDBConectWriterImpl con_writer(outstream_,atom_indices_);
ent.Apply(con_writer);
this->WriteModelTrailer();
}
void PDBWriter::Write(const mol::EntityView& ent)
{
this->WriteModel(ent);
}
void PDBWriter::Write(const mol::EntityHandle& ent)
{
this->WriteModel(ent);
}
void PDBWriter::Write(const mol::AtomHandleList& atoms)
{
ForcePOSIX posix = ForcePOSIX();
this->WriteModelLeader();
int counter=1;
mol::ChainHandle last_chain;
for (mol::AtomHandleList::const_iterator i=atoms.begin(),
e=atoms.end(); i!=e; ++i, ++counter) {
write_atom(outstream_, line_, *i, counter,
(PDB::Flags() & PDB::PQR_FORMAT) != 0, charmm_style_);
}
this->WriteModelTrailer();
}
PDBWriter::~PDBWriter()
{
outstream_ << "END ";
}
}}