dcd.hpp

/*
  This file is part of LOOS.
 
  LOOS (Lightweight Object-Oriented Structure library)
  Copyright (c) 2008, Tod D. Romo, Alan Grossfield
  Department of Biochemistry and Biophysics
  School of Medicine & Dentistry, University of Rochester
 
  This package (LOOS) is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation under version 3 of the License.
 
  This package 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 General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/
 
 
 
#if !defined(LOOS_DCD_HPP)
#define LOOS_DCD_HPP
 
 
#include <iostream>
#include <string>
#include <stdexcept>
#include <exception>
 
#include <loos_defs.hpp>
 
#include <Trajectory.hpp>
 
 
namespace loos {
 
 
 
    class DCD : public Trajectory {
        static bool suppress_warnings;
 
 
        // Use a union to convert data to appropriate type...
        typedef union { unsigned int ui; int i; char c[4]; float f; } DataOverlay;
 
    public:
        class EndOfFile : public LOOSError {
        public:
            EndOfFile() : LOOSError("unexpected end of file while reading DCD") { }
        };
 
        explicit DCD(const std::string s) :  Trajectory(s), _natoms(0), _nframes(0),
                                             qcrys(std::vector<double>(6)),
                                             frame_size(0), first_frame_pos(0),
                                             swabbing(false) { initTrajectory(); }
 
        explicit DCD(const char* s) :  Trajectory(s), _natoms(0), _nframes(0),
                                       qcrys(std::vector<double>(6)), frame_size(0),
                                       first_frame_pos(0), swabbing(false) { initTrajectory(); }
 
        explicit DCD(std::istream& fs) : Trajectory(fs), _natoms(0), _nframes(0),
                                         qcrys(std::vector<double>(6)), frame_size(0), first_frame_pos(0),
                                         swabbing(false) { initTrajectory(); };
 
        std::string description() const { return("CHARMM/NAMD DCD"); }
 
        static pTraj create(const std::string& fname, const AtomicGroup& model) {
            return(pTraj(new DCD(fname)));
        }
 
 
 
        // Accessor methods...
 
        virtual uint natoms(void) const;
        virtual bool hasPeriodicBox(void) const;
        virtual GCoord periodicBox(void) const;
 
        virtual bool hasVelocities() const { return(false); }
        virtual double velocityConversionFactor() const { return(20.45482706); }
 
        std::vector<std::string> titles(void) const;
 
        int icntrl(const int) const;
        void icntrl(const int, const int);
 
        // * legacy *
        std::vector<double> crystalParams(void) const;
        bool hasCrystalParams(void) const;
 
        virtual float timestep(void) const;
        virtual uint nframes(void) const;
 
        std::vector<dcd_real> xcoords(void) const;
        std::vector<dcd_real> ycoords(void) const;
        std::vector<dcd_real> zcoords(void) const;
 
        // The following track CHARMm names (more or less...)
        unsigned int nsteps(void) const;
        float delta(void) const;
        int nsavc(void) const;
        int nfile(void) const;
        int nfixed(void) const;
 
        bool nativeFormat(void) const;
 
 
        virtual std::vector<GCoord> coords(void) const;
 
        // This is slated to go away...
        std::vector<GCoord> mappedCoords(const std::vector<int>& map);
 
 
 
        static void setSuppression(const bool b) { suppress_warnings = b; }
 
        virtual bool parseFrame(void);
 
    private:
 
        void readHeader(void);
 
        void initTrajectory();
 
        uint calculateNumberOfFrames();
 
        // Trajectory member functions we must provide...
        virtual void seekNextFrameImpl(void) { }    // DCD frames are always contiguous, so do nothing...
        virtual void seekFrameImpl(const uint);
 
        virtual void rewindImpl(void);
 
        virtual void updateGroupCoordsImpl(AtomicGroup& g);
 
 
 
        void allocateSpace(const int n);
        bool readCrystalParams(void);
        bool readCoordLine(std::vector<float>& v);
 
        void endianMatch(pStream& fsw);
 
        // For reading F77 I/O
        unsigned int readRecordLen(void);
        DataOverlay* readF77Line(unsigned int *len);
 
 
 
    private:
        int _icntrl[20];          // DCD header data
        uint _natoms;              // # of atoms
        uint _nframes;
        std::vector<std::string> _titles;   // Vector of title lines from DCD
        std::vector<double> qcrys;     // Crystal params
        float _delta;             // Timestep (extracted from _icntrl)
 
        std::streamoff frame_size;                    // *Internal* size (in bytes) of each frame
        std::streampos first_frame_pos;     // *Internal* location in file of start of data frames
 
        bool swabbing;            // DCD being read is not in native format...
 
        std::vector<dcd_real> xcrds, ycrds, zcrds;
 
    };
 
}
 
 
 
#endif