libhoomd/analyzers/PDBDumpWriter.cc

   1 /*
   2 Highly Optimized Object-oriented Many-particle Dynamics -- Blue Edition
   3 (HOOMD-blue) Open Source Software License Copyright 2009-2014 The Regents of
   4 the University of Michigan All rights reserved.
   5
   6 HOOMD-blue may contain modifications ("Contributions") provided, and to which
   7 copyright is held, by various Contributors who have granted The Regents of the
   8 University of Michigan the right to modify and/or distribute such Contributions.
   9
  10 You may redistribute, use, and create derivate works of HOOMD-blue, in source
  11 and binary forms, provided you abide by the following conditions:
  12
  13 * Redistributions of source code must retain the above copyright notice, this
  14 list of conditions, and the following disclaimer both in the code and
  15 prominently in any materials provided with the distribution.
  16
  17 * Redistributions in binary form must reproduce the above copyright notice, this
  18 list of conditions, and the following disclaimer in the documentation and/or
  19 other materials provided with the distribution.
  20
  21 * All publications and presentations based on HOOMD-blue, including any reports
  22 or published results obtained, in whole or in part, with HOOMD-blue, will
  23 acknowledge its use according to the terms posted at the time of submission on:
  24 http://codeblue.umich.edu/hoomd-blue/citations.html
  25
  26 * Any electronic documents citing HOOMD-Blue will link to the HOOMD-Blue website:
  27 http://codeblue.umich.edu/hoomd-blue/
  28
  29 * Apart from the above required attributions, neither the name of the copyright
  30 holder nor the names of HOOMD-blue's contributors may be used to endorse or
  31 promote products derived from this software without specific prior written
  32 permission.
  33
  34 Disclaimer
  35
  36 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS ``AS IS'' AND
  37 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  38 WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, AND/OR ANY
  39 WARRANTIES THAT THIS SOFTWARE IS FREE OF INFRINGEMENT ARE DISCLAIMED.
  40
  41 IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
  42 INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  43 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  44 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  45 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
  46 OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
  47 ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  48 */
  49
  50 // Maintainer: joaander
  51
  52 /*! \file PDBDumpWriter.cc
  53     \brief Defines the PDBDumpWriter class
  54 */
  55
  56 // this file was written from scratch, but some error checks and PDB snprintf strings were copied from VMD's molfile plugin
  57 // it is used under the following license
  58
  59 // University of Illinois Open Source License
  60 // Copyright 2003 Theoretical and Computational Biophysics Group,
  61 // All rights reserved.
  62
  63 // Developed by:       Theoretical and Computational Biophysics Group
  64 //             University of Illinois at Urbana-Champaign
  65 //            http://www.ks.uiuc.edu/
  66
  67 #ifdef WIN32
  68 #pragma warning( push )
  69 #pragma warning( disable : 4244 )
  70 #endif
  71
  72 #include <boost/python.hpp>
  73 using namespace boost::python;
  74
  75 #include <string>
  76 #include <fstream>
  77 #include <iomanip>
  78 #include <stdio.h>
  79
  80 #include "PDBDumpWriter.h"
  81 #include "BondedGroupData.h"
  82
  83 using namespace std;
  84 using namespace boost;
  85
  86 /*! \param sysdef System definition containing particle data to write
  87     \param base_fname Base filename to expand with **timestep**.pdb when writing
  88 */
  89 PDBDumpWriter::PDBDumpWriter(boost::shared_ptr<SystemDefinition> sysdef, std::string base_fname)
  90         : Analyzer(sysdef), m_base_fname(base_fname), m_output_bond(false)
  91     {
  92     m_exec_conf->msg->notice(5) << "Constructing PDBDumpWriter: " << base_fname << endl;
  93     }
  94
  95 PDBDumpWriter::~PDBDumpWriter()
  96     {
  97     m_exec_conf->msg->notice(5) << "Destroying PDBDumpWriter" << endl;
  98     }
  99
 100 /*! \param timestep Current time step of the simulation
 101
 102     analzye() constructs af file name m_base_fname.**timestep**.pdb and writes out the
 103     current state of the system to that file.
 104 */
 105 void PDBDumpWriter::analyze(unsigned int timestep)
 106     {
 107     if (m_prof)
 108         m_prof->push("Dump PDB");
 109
 110     ostringstream full_fname;
 111     string filetype = ".pdb";
 112
 113     // Generate a filename with the timestep padded to ten zeros
 114     full_fname << m_base_fname << "." << setfill('0') << setw(10) << timestep << filetype;
 115     // then write the file
 116     writeFile(full_fname.str());
 117
 118     if (m_prof)
 119         m_prof->pop();
 120     }
 121
 122 /*! \param fname File name to write
 123
 124     Writes the current state of the system to the pdb file \a fname
 125 */
 126 void PDBDumpWriter::writeFile(std::string fname)
 127     {
 128     // open the file for writing
 129     ofstream f(fname.c_str());
 130     f.exceptions ( ifstream::eofbit | ifstream::failbit | ifstream::badbit );
 131
 132     if (!f.good())
 133         {
 134         m_exec_conf->msg->error() << "dump.pdb: Unable to open file for writing: " << fname << endl;
 135         throw runtime_error("Error writting pdb dump file");
 136         }
 137
 138     // acquire the particle data
 139     ArrayHandle<Scalar4> h_pos(m_pdata->getPositions(), access_location::host, access_mode::read);
 140     ArrayHandle<unsigned int> h_rtag(m_pdata->getRTags(), access_location::host, access_mode::read);
 141
 142
 143     // get the box dimensions
 144     BoxDim box = m_pdata->getBox();
 145
 146     // start writing the heinous PDB format
 147     const int linesize = 82;
 148     char buf[linesize];
 149
 150     // output the box dimensions
 151     Scalar a,b,c,alpha,beta,gamma;
 152     Scalar3 va = box.getLatticeVector(0);
 153     Scalar3 vb = box.getLatticeVector(1);
 154     Scalar3 vc = box.getLatticeVector(2);
 155     a = sqrt(dot(va,va));
 156     b = sqrt(dot(vb,vb));
 157     c = sqrt(dot(vc,vc));
 158     alpha = 90.0 - asin(dot(vb,vc)/(b*c)) * 90.0 / M_PI_2;
 159     beta = 90.0 - asin(dot(va,vc)/(a*c)) * 90.0 / M_PI_2;
 160     gamma = 90.0 - asin(dot(va,vb)/(a*b)) * 90.0 / M_PI_2;
 161
 162     snprintf(buf, linesize, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n", a,b,c, alpha, beta, gamma);
 163     f << buf;
 164
 165     // write out all the atoms
 166     for (unsigned int j = 0; j < m_pdata->getN(); j++)
 167         {
 168         int i;
 169         i= h_rtag.data[j];
 170
 171         // first check that everything will fit into the PDB output
 172         if (h_pos.data[i].x < -999.9994f || h_pos.data[i].x > 9999.9994f || h_pos.data[i].y < -999.9994f || h_pos.data[i].y > 9999.9994f || h_pos.data[i].z < -999.9994f || h_pos.data[i].z > 9999.9994f)
 173             {
 174             m_exec_conf->msg->error() << "dump.pdb: Coordinate " << h_pos.data[i].x << " " << h_pos.data[i].y << " " << h_pos.data[i].z << " is out of range for PDB writing" << endl;
 175             throw runtime_error("Error writing PDB file");
 176             }
 177         // check the length of the type name
 178         const string &type_name = m_pdata->getNameByType(__scalar_as_int(h_pos.data[i].w));
 179         if (type_name.size() > 4)
 180             {
 181             m_exec_conf->msg->error() << "dump.pdb: Type " << type_name << " is too long for PDB writing" << endl;
 182             throw runtime_error("Error writing PDB file");
 183             }
 184
 185         // start preparing the stuff to write (copied from VMD's molfile plugin)
 186         char indexbuf[32];
 187         char residbuf[32];
 188         char segnamebuf[5];
 189         char altlocchar;
 190
 191         /* XXX                                                          */
 192         /* if the atom or residue indices exceed the legal PDB spec, we */
 193         /* start emitting asterisks or hexadecimal strings rather than  */
 194         /* aborting.  This is not really legal, but is an accepted hack */
 195         /* among various other programs that deal with large PDB files  */
 196         /* If we run out of hexadecimal indices, then we just print     */
 197         /* asterisks.                                                   */
 198         if (i < 100000)
 199             {
 200             sprintf(indexbuf, "%5d", i);
 201             }
 202         else if (i < 1048576)
 203             {
 204             sprintf(indexbuf, "%05x", i);
 205             }
 206         else
 207             {
 208             sprintf(indexbuf, "*****");
 209             }
 210
 211         /*if (resid < 10000) {
 212         sprintf(residbuf, "%4d", resid);
 213         } else if (resid < 65536) {
 214         sprintf(residbuf, "%04x", resid);
 215         } else {
 216         sprintf(residbuf, "****");
 217         }*/
 218         sprintf(residbuf, "%4d", 1);
 219
 220         //altlocchar = altloc[0];
 221         //if (altlocchar == '\0') {
 222         altlocchar = ' ';
 223         //}
 224
 225         /* make sure the segname does not overflow the format */
 226         sprintf(segnamebuf, "SEG ");
 227
 228         snprintf(buf, linesize, "%-6s%5s %4s%c%-4s%c%4s%c   %8.3f%8.3f%8.3f%6.2f%6.2f      %-4s%2s\n", "ATOM  ", indexbuf, type_name.c_str(), altlocchar, "RES", ' ',
 229                  residbuf, ' ', h_pos.data[i].x, h_pos.data[i].y, h_pos.data[i].z, 0.0f, 0.0f, segnamebuf, "  ");
 230         f << buf;
 231         }
 232
 233     if (m_output_bond)
 234         {
 235         // error check: pdb files cannot contain bonds with 100,000 or more atom records
 236         if (m_pdata->getN() >= 100000)
 237             {
 238             m_exec_conf->msg->error() << "dump.pdb: PDB files with bonds cannot hold more than 99,999 atoms!" << endl;
 239             throw runtime_error("Error dumping PDB file");
 240             }
 241
 242         // grab the bond data
 243         boost::shared_ptr<BondData> bond_data = m_sysdef->getBondData();
 244         for (unsigned int i = 0; i < bond_data->getN(); i++)
 245             {
 246             BondData::members_t bond = bond_data->getMembersByIndex(i);
 247             snprintf(buf, linesize, "CONECT%1$5d%2$5d\n", bond.tag[0], bond.tag[1]);
 248             f << buf;
 249             }
 250         }
 251
 252     }
 253
 254 void export_PDBDumpWriter()
 255     {
 256     class_<PDBDumpWriter, boost::shared_ptr<PDBDumpWriter>, bases<Analyzer>, boost::noncopyable>
 257     ("PDBDumpWriter", init< boost::shared_ptr<SystemDefinition>, std::string >())
 258     .def("setOutputBond", &PDBDumpWriter::setOutputBond)
 259     .def("writeFile", &PDBDumpWriter::writeFile)
 260     ;
 261     }
 262
 263 #ifdef WIN32
 264 #pragma warning( pop )
 265 #endif