diff --git a/.travis.yml b/.travis.yml index ebbe1bc143..71169d0766 100644 --- a/.travis.yml +++ b/.travis.yml @@ -36,7 +36,8 @@ matrix: before_install: - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then - brew install netcdf fftw; + brew update; + brew install gcc netcdf fftw; curl -L https://anaconda.org/AmberMD/amber_phenix/0.9.6/download/osx-64/amber_phenix-0.9.6-0.tar.bz2 > $HOME/osx_amber.tar.bz2; tar jxf $HOME/osx_amber.tar.bz2 lib/libsander.dylib AmberTools/src/sander/sander.h; diff --git a/doc/CpptrajManual.pdf b/doc/CpptrajManual.pdf index 1a2e5036c0..4848d14e50 100644 Binary files a/doc/CpptrajManual.pdf and b/doc/CpptrajManual.pdf differ diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 15b0c9cd9a..39bf9e9a82 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -13976,6 +13976,10 @@ Options for pdb format: [bfacscale] [occscale] [bfacmax ] [occmax ] \end_layout +\begin_layout LyX-Code + [adpdata ] +\end_layout + \begin_deeper \begin_layout Description dumpq PQR format; write charges (in units of e-) and GB radii to occupancy @@ -13992,6 +13996,8 @@ vdw PQR format; write charges and vdW radii to occupancy/B-factor columns. \begin_layout Description pdbres: Use PDB V3 residue names. + Will write a default chain ID ('Z') for each residue if the corresponding + topology does not have chain ID information. \end_layout \begin_layout Description @@ -14000,6 +14006,7 @@ pdbatom: Use PDB V3 atom names. \begin_layout Description pdbv3: Use PDB V3 residue/atom names. + Same as specifying 'pdbres' and 'pdbatom'. \end_layout \begin_layout Description @@ -14112,13 +14119,32 @@ occmax Max value for occscale. \end_layout -\end_deeper -\begin_layout Standard -Note that if the input topology does not contain chain IDs, a default chain - ID of 'Z' will be used if no other chain ID option is specified in order - to remain consistent with the PDB standard. +\begin_layout Description +adpdata +\begin_inset space ~ +\end_inset + + Use data in (e.g. + from the +\series bold +\emph on +atomicfluct +\series default +\emph default + command, +\begin_inset CommandInset ref +LatexCommand vpageref +reference "subsec:cpptraj_atomicfluct" +plural "false" +caps "false" +noprefix "false" + +\end_inset + +) for anisotropic B-factors. \end_layout +\end_deeper \begin_layout Subsubsection \shape italic @@ -17643,11 +17669,11 @@ atomicfluct \end_inset - [out ] [] [byres | byatom | bymask] [bfactor] + [] [out ] [] [byres | byatom | bymask] \end_layout \begin_layout LyX-Code - [calcadp [adpout ]] + [bfactor] [calcadp [adpout ]] \end_layout \begin_layout LyX-Code @@ -17659,6 +17685,10 @@ atomicfluct \end_layout \begin_deeper +\begin_layout Description + Output data set name. +\end_layout + \begin_layout Description out \begin_inset space ~ @@ -17722,6 +17752,18 @@ graphic B-factors. [] Frames to skip between calculations (default 1). \end_layout +\begin_layout Standard +DataSets created +\end_layout + +\begin_layout Description + Hold atomic fluctuations. +\end_layout + +\begin_layout Description +[ADP] Hold anisotropic displacement parameters if 'calcadp' specified. +\end_layout + \end_deeper \begin_layout Standard Compute the atomic positional fluctuations (also referred to as root-mean-square @@ -17796,6 +17838,7 @@ calcadp adpout \series default (or STDOUT if not specified) using PDB ANISOU record format. + The displacement factors will be saved to a data set. Note that \series bold calcadp @@ -17862,6 +17905,39 @@ To write the RMSF or atomic positional fluctuations of the same atoms, use atomicfluct out backbone-atoms.agr @C,CA,N \end_layout +\begin_layout Standard +To write a PDB of averaged coordinates (after fitting to the first frame) + with both B-factors and anisotropic temperature factors: +\end_layout + +\begin_layout LyX-Code +parm myparm.parm7 +\end_layout + +\begin_layout LyX-Code +trajin mytraj.nc +\end_layout + +\begin_layout LyX-Code +rms first +\end_layout + +\begin_layout LyX-Code +average crdset MyAvg +\end_layout + +\begin_layout LyX-Code +atomicfluct MyFluct calcadp +\end_layout + +\begin_layout LyX-Code +run +\end_layout + +\begin_layout LyX-Code +crdout MyAvg mypdb.pdb adpdata MyFluct[ADP] bfacdata MyFluct +\end_layout + \begin_layout Subsection atommap \end_layout diff --git a/src/Action_AtomicFluct.cpp b/src/Action_AtomicFluct.cpp index 591d245880..8d270636f8 100644 --- a/src/Action_AtomicFluct.cpp +++ b/src/Action_AtomicFluct.cpp @@ -14,12 +14,13 @@ Action_AtomicFluct::Action_AtomicFluct() : adpoutfile_(0), fluctParm_(0), outtype_(BYATOM), - dataout_(0) + dataout_(0), + adpset_(0) {} void Action_AtomicFluct::Help() const { - mprintf("\t[out ] [] [byres [pdbres] | byatom | bymask] [bfactor]\n" - "\t[calcadp [adpout ]]\n" + mprintf("\t[] [out ] [] [byres [pdbres] | byatom | bymask]\n" + "\t[bfactor] [calcadp [adpout ]]\n" "\t%s\n" " Calculate atomic fluctuations of atoms in \n", ActionFrameCounter::HelpText); } @@ -50,8 +51,10 @@ Action::RetType Action_AtomicFluct::Init(ArgList& actionArgs, ActionInit& init, std::string setname = actionArgs.GetStringNext(); // Add output dataset MetaData md( setname, "", MetaData::NOT_TS ); + bool setLegend = false; if (setname.empty()) { // Only overwrite legend if no name specified. + setLegend = true; if (bfactor_) md.SetLegend("B-factors"); else @@ -62,8 +65,20 @@ Action::RetType Action_AtomicFluct::Init(ArgList& actionArgs, ActionInit& init, mprinterr("Error: AtomicFluct: Could not allocate dataset for output.\n"); return Action::ERR; } + if (calc_adp_) { + MetaData md2(dataout_->Meta().Name(), "ADP"); + if (setLegend) + md2.SetLegend("Aniso. B-factors"); + adpset_ = init.DSL().AddSet( DataSet::TENSOR, md2 ); + if (adpset_ == 0) { + mprinterr("Error: Could not allocate ADP dataset.\n"); + return Action::ERR; + } + adpset_->ModifyDim(Dimension::X).SetLabel("Atom"); + } # ifdef MPI dataout_->SetNeedsSync( false ); // Not a time series + if (adpset_ != 0) adpset_->SetNeedsSync( false ); trajComm_ = init.TrajComm(); # endif if (outfile != 0) @@ -201,24 +216,30 @@ void Action_AtomicFluct::Print() { int idx = (i/3); int atom = Mask_[idx]; int resnum = (*fluctParm_)[atom].ResNum(); - int u11 = (int)(SumCoords2_[i ] * 10000); - int u22 = (int)(SumCoords2_[i+1] * 10000); - int u33 = (int)(SumCoords2_[i+2] * 10000); + // Calculate the anisotropic factors; + double anisou[6]; + anisou[0] = SumCoords2_[i ]; // u11 + anisou[1] = SumCoords2_[i+1]; // u22 + anisou[2] = SumCoords2_[i+2]; // u33 // Calculate covariance: - etc. - int u12 = (int)((Cross_[i ] - SumCoords_[i ] * SumCoords_[i+1]) * 10000); - int u13 = (int)((Cross_[i+1] - SumCoords_[i ] * SumCoords_[i+2]) * 10000); - int u23 = (int)((Cross_[i+2] - SumCoords_[i+1] * SumCoords_[i+2]) * 10000); - // To make PDB as compliant as possible, ensure there is a chain ID. + anisou[3] = (Cross_[i ] - SumCoords_[i ] * SumCoords_[i+1]); // u12 + anisou[4] = (Cross_[i+1] - SumCoords_[i ] * SumCoords_[i+2]); // u13 + anisou[5] = (Cross_[i+2] - SumCoords_[i+1] * SumCoords_[i+2]); // u23 + // Store as data; index by actual atom number + adpset_->Add(atom+1, (const void*)&anisou); + // Default to a blank chain ID char chainid; if (fluctParm_->Res(resnum).HasChainID()) chainid = fluctParm_->Res(resnum).ChainID(); else - chainid = Residue::DefaultChainID(); - PDBfile& adpout = static_cast( *adpoutfile_ ); - adpout.WriteANISOU( - atom+1, (*fluctParm_)[atom].c_str(), fluctParm_->Res(resnum).c_str(), - chainid, fluctParm_->Res(resnum).OriginalResNum(), - u11, u22, u33, u12, u13, u23, (*fluctParm_)[atom].ElementName(), 0 ); + chainid = ' '; + if (adpoutfile_ != 0) { + PDBfile& adpout = static_cast( *adpoutfile_ ); + adpout.WriteANISOU( + atom+1, (*fluctParm_)[atom].c_str(), fluctParm_->Res(resnum).c_str(), + chainid, fluctParm_->Res(resnum).OriginalResNum(), + anisou, (*fluctParm_)[atom].ElementName(), 0 ); + } } } } else { diff --git a/src/Action_AtomicFluct.h b/src/Action_AtomicFluct.h index 7f956fb84f..8541739988 100644 --- a/src/Action_AtomicFluct.h +++ b/src/Action_AtomicFluct.h @@ -32,5 +32,6 @@ class Action_AtomicFluct : public Action, ActionFrameCounter { Topology *fluctParm_; outputType outtype_; DataSet* dataout_; + DataSet* adpset_; ///< DataSet for holding anisotropic temp. factors }; #endif diff --git a/src/DataSet.cpp b/src/DataSet.cpp index 3fa4af9bfa..22cc7a4c1f 100644 --- a/src/DataSet.cpp +++ b/src/DataSet.cpp @@ -28,7 +28,8 @@ const char* DataSet::Descriptions_[] = { "pH", // PH "pH REMD (explicit)", // PH_EXPL "pH REMD (implicit)", // PH_IMPL - "parameters" // PARAMETERS + "parameters", // PARAMETERS + "tensor" // TENSOR }; // CONSTRUCTOR diff --git a/src/DataSet.h b/src/DataSet.h index 54a36b5b4d..f4c3714792 100644 --- a/src/DataSet.h +++ b/src/DataSet.h @@ -22,12 +22,13 @@ class DataSet { public: typedef DataSet* (*AllocatorType)(); typedef std::vector SizeArray; - /// DataSet base type. + /// DataSet base type. + /** When adding new entries make sure that Descriptions_ is updated. */ enum DataType { UNKNOWN_DATA=0, DOUBLE, FLOAT, INTEGER, STRING, MATRIX_DBL, MATRIX_FLT, COORDS, VECTOR, MODES, GRID_FLT, GRID_DBL, REMLOG, XYMESH, TRAJ, REF_FRAME, MAT3X3, TOPOLOGY, CMATRIX, CMATRIX_NOMEM, CMATRIX_DISK, PH, PH_EXPL, PH_IMPL, - PARAMETERS + PARAMETERS, TENSOR }; /// Group DataSet belongs to. enum DataGroup { diff --git a/src/DataSetList.cpp b/src/DataSetList.cpp index 16ba773e0d..b1a225ceef 100644 --- a/src/DataSetList.cpp +++ b/src/DataSetList.cpp @@ -30,6 +30,7 @@ #include "DataSet_PHREMD_Explicit.h" #include "DataSet_PHREMD_Implicit.h" #include "DataSet_Parameters.h" +#include "DataSet_Tensor.h" bool DataSetList::useDiskCache_ = false; @@ -76,6 +77,7 @@ DataSet* DataSetList::NewSet(DataSet::DataType typeIn) { case DataSet::PH_EXPL : ds = DataSet_PHREMD_Explicit::Alloc(); break; case DataSet::PH_IMPL : ds = DataSet_PHREMD_Implicit::Alloc(); break; case DataSet::PARAMETERS : ds = DataSet_Parameters::Alloc(); break; + case DataSet::TENSOR : ds = DataSet_Tensor::Alloc(); break; // Sanity check default: mprinterr("Internal Error: No allocator for DataSet type '%s'\n", diff --git a/src/DataSet_Tensor.cpp b/src/DataSet_Tensor.cpp new file mode 100644 index 0000000000..b572e63d5b --- /dev/null +++ b/src/DataSet_Tensor.cpp @@ -0,0 +1,68 @@ +#include // std::copy +#include "DataSet_Tensor.h" + +/** CONSTRUCTOR. */ +DataSet_Tensor::DataSet_Tensor() : + DataSet(TENSOR, GENERIC, TextFormat(TextFormat::DOUBLE, 8, 4, 6), 1) +{} + +#ifdef MPI +int DataSet_Tensor::Sync(size_t, std::vector const&, Parallel::Comm const&) +{ + return 1; +} +#endif + +/** Reserve memory for given input size. */ +int DataSet_Tensor::Allocate( SizeArray const& sizeIn ) { + if (!sizeIn.empty()) { + Xvals_.reserve( sizeIn[0] ); + Data_.reserve( sizeIn[0] ); + } + return 0; +} + +/** Add index and tensor. */ +void DataSet_Tensor::Add(size_t frame, const void* tIn) { + Xvals_.push_back( (double)frame ); + Data_.push_back( Ttype( (const double*)tIn ) ); +} + +/** Write tensor. */ +void DataSet_Tensor::WriteBuffer(CpptrajFile &cbuffer, SizeArray const& pIn) const { + if (pIn[0] >= Data_.size()) { + cbuffer.Printf(format_.fmt(), 0.0, 0.0, 0.0, 0.0, 0.0, 0.0); + } else { + const double* ptr = Data_[ pIn[0] ].Ptr(); + cbuffer.Printf(format_.fmt(), ptr[0], ptr[1], ptr[2], + ptr[3], ptr[4], ptr[5]); + } +} + +/** \return specified index value. */ +double DataSet_Tensor::Coord(unsigned int d, size_t p) const { + if (d > 0) return 0.0; + return Xvals_[p]; +} + +/** Append given Tensor set to this one. */ +int DataSet_Tensor::Append(DataSet* dsIn) { + if (dsIn == 0 || dsIn->Empty()) return 0; + if (dsIn->Type() != TENSOR) return 1; + Tarray const& tIn = ((DataSet_Tensor*)dsIn)->Data_; + Farray const& fIn = ((DataSet_Tensor*)dsIn)->Xvals_; + size_t oldsize = Data_.size(); + Data_.resize( oldsize + tIn.size() ); + std::copy( tIn.begin(), tIn.end(), Data_.begin() + oldsize ); + Xvals_.resize( oldsize + fIn.size() ); // TODO check that tIn and fIn size match? + std::copy( fIn.begin(), fIn.end(), Xvals_.begin() + oldsize ); + + return 0; +} + +/** \return Memory usage in bytes. */ +size_t DataSet_Tensor::MemUsageInBytes() const { + return ( sizeof(Tarray) + + sizeof(Farray) + + (Data_.size() + Xvals_.size()) * sizeof(double) ); +} diff --git a/src/DataSet_Tensor.h b/src/DataSet_Tensor.h new file mode 100644 index 0000000000..e22956392f --- /dev/null +++ b/src/DataSet_Tensor.h @@ -0,0 +1,43 @@ +#ifndef INC_DATASET_TENSOR_H +#define INC_DATASET_TENSOR_H +#include +#include "SymmetricTensor.h" +#include "DataSet.h" +/// Hold an array of 3x3 symmetric tensor values. The array can be sparse. +/** NOTE: We use a floating point type array for indices in case we want to + * hold time values at some point. + */ +class DataSet_Tensor : public DataSet { + public: + typedef SymmetricTensor Ttype; + + /// CONSTRUCTOR + DataSet_Tensor(); + /// Allocator + static DataSet* Alloc() { return (DataSet*)new DataSet_Tensor(); } + + /// \return Index at specified position + double Xvals(unsigned int idx) const { return Xvals_[idx]; } + /// \return Tensor at specified position + Ttype Tensor(unsigned int idx) const { return Data_[idx]; } + + // ----- DataSet functions ------------------- + size_t Size() const { return Data_.size(); } +# ifdef MPI + int Sync(size_t, std::vector const&, Parallel::Comm const&); +# endif + void Info() const { return; } + int Allocate(SizeArray const&); + void Add( size_t, const void* ); + void WriteBuffer(CpptrajFile&, SizeArray const&) const; + double Coord(unsigned int, size_t) const; + int Append(DataSet*); + size_t MemUsageInBytes() const; + private: + typedef std::vector Tarray; + typedef std::vector Farray; + + Farray Xvals_; /// < Hold the index values + Tarray Data_; ///< Hold the tensor values +}; +#endif diff --git a/src/PDBfile.cpp b/src/PDBfile.cpp index e2ae7ab65c..888b97352a 100644 --- a/src/PDBfile.cpp +++ b/src/PDBfile.cpp @@ -441,6 +441,28 @@ static inline void x_to_buf(char* coord_buf, double X, bool& coordOverflow) sprintf(coord_buf, "%8.3f", X); } +static inline void SetChargeString(char* charge_buf, int charge) { + charge_buf[0] = ' '; + charge_buf[1] = ' '; + charge_buf[2] = '\0'; + if (charge > 0) { + if (charge > 9) { + mprintf("Warning: Charge %i is too large. Not printing.\n", charge); + } else { + charge_buf[0] = (char)(charge + '0'); + charge_buf[1] = '+'; + } + } else if (charge < 0) { + if (charge < -9) { + mprintf("Warning: Charge %i is too large. Not printing.\n", charge); + } else { + charge = -charge; + charge_buf[0] = (char)(charge + '0'); + charge_buf[1] = '-'; + } + } +} + // PDBfile::WriteCoord() void PDBfile::WriteCoord(PDB_RECTYPE Record, int anum, NameType const& name, char altLoc, @@ -456,21 +478,32 @@ void PDBfile::WriteCoord(PDB_RECTYPE Record, int anum, NameType const& name, x_to_buf(coord_buf+8 , Y, coordOverflow_); x_to_buf(coord_buf+16, Z, coordOverflow_); WriteRecordHeader(Record, anum, name, altLoc, resnameIn, chain, resnum, icode, Elt); + static char charge_buf[3]; + SetChargeString(charge_buf, charge); if (highPrecision) - Printf(" %24s%8.4f%8.4f %2s%2s\n", coord_buf, Occ, B, Elt, ""); + Printf(" %24s%8.4f%8.4f %2s%2s\n", coord_buf, Occ, B, Elt, charge_buf); else - Printf(" %24s%6.2f%6.2f %2s%2s\n", coord_buf, Occ, B, Elt, ""); + Printf(" %24s%6.2f%6.2f %2s%2s\n", coord_buf, Occ, B, Elt, charge_buf); } // PDBfile::WriteANISOU() void PDBfile::WriteANISOU(int anum, NameType const& name, NameType const& resnameIn, char chain, int resnum, - int u11, int u22, int u33, int u12, int u13, int u23, + const double* anisou, const char* Elt, int charge) { // TODO icode, altLoc + // Convert to integers + int u11 = (int)(anisou[0] * 10000); + int u22 = (int)(anisou[1] * 10000); + int u33 = (int)(anisou[2] * 10000); + int u12 = (int)(anisou[3] * 10000); + int u13 = (int)(anisou[4] * 10000); + int u23 = (int)(anisou[5] * 10000); + static char charge_buf[3]; + SetChargeString(charge_buf, charge); WriteRecordHeader(ANISOU, anum, name, ' ', resnameIn, chain, resnum, ' ', Elt); - Printf(" %7i%7i%7i%7i%7i%7i %2s%2i\n", u11, u22, u33, - u12, u13, u23, Elt, charge); + Printf(" %7i%7i%7i%7i%7i%7i %2s%2s\n", u11, u22, u33, + u12, u13, u23, Elt, charge_buf); } // PDBfile::WriteTITLE() diff --git a/src/PDBfile.h b/src/PDBfile.h index b691ff59e2..2dc4199b93 100644 --- a/src/PDBfile.h +++ b/src/PDBfile.h @@ -64,7 +64,7 @@ class PDBfile : public CpptrajFile { bool CoordOverflow() { bool stat = coordOverflow_; coordOverflow_ = false; return stat; } /// Write ANISOU record. void WriteANISOU(int, NameType const&, NameType const&, char, int, - int, int, int, int, int, int, const char *, int); + const double*, const char *, int); /// Write TITLE void WriteTITLE(std::string const&); /// Write CRYST1 diff --git a/src/SymmetricTensor.h b/src/SymmetricTensor.h new file mode 100644 index 0000000000..b6c10b7640 --- /dev/null +++ b/src/SymmetricTensor.h @@ -0,0 +1,35 @@ +#ifndef INC_SYMMETRICTENSOR_H +#define INC_SYMMETRICTENSOR_H +/// Class for holding elements of a 3x3 symmetric tensor (6 elements total) +/** U11 U12 U13 + * U12 U22 U23 + * U13 U23 U33 + * stored as + * U11 U22 U33 U12 U13 U23 + */ +template class SymmetricTensor { + public: + /// CONSTRUCTOR + SymmetricTensor() {} + /// CONSTRUCTOR - take pointer to array of size 6 + SymmetricTensor(T const* p) { + U_[0] = p[0]; + U_[1] = p[1]; + U_[2] = p[2]; + U_[3] = p[3]; + U_[4] = p[4]; + U_[5] = p[5]; + } + + T const& U11() const { return U_[0]; } + T const& U22() const { return U_[1]; } + T const& U33() const { return U_[2]; } + T const& U12() const { return U_[3]; } + T const& U13() const { return U_[4]; } + T const& U23() const { return U_[5]; } + + T const* Ptr() const { return U_; } + private: + T U_[6]; +}; +#endif diff --git a/src/Traj_PDBfile.cpp b/src/Traj_PDBfile.cpp index 2bd0cb662d..82674a43c8 100644 --- a/src/Traj_PDBfile.cpp +++ b/src/Traj_PDBfile.cpp @@ -7,6 +7,7 @@ #include "DistRoutines.h" #include "Constants.h" #include "DataSet_1D.h" // for bfacdata, occdata +#include "DataSet_Tensor.h" // for adpdata // CONSTRUCTOR Traj_PDBfile::Traj_PDBfile() : @@ -32,6 +33,7 @@ Traj_PDBfile::Traj_PDBfile() : chainchar_(Residue::BlankChainID()), bfacdata_(0), occdata_(0), + adpdata_(0), bfacmax_(99.99), occmax_(99.99) {} @@ -204,6 +206,7 @@ void Traj_PDBfile::WriteHelp() { "\toccscale : If specified scale values in occupancy column between 0 and .\n" "\tbfacmax : Max value for bfacscale.\n" "\toccmax : Max value for occscale.\n" + "\tadpdata : Use data in for anisotropic B-factors.\n" ); } @@ -285,6 +288,18 @@ int Traj_PDBfile::processWriteArgs(ArgList& argIn, DataSetList const& DSLin) { if (bfacscale_) bfacmax_ = argIn.getKeyDouble("bfacmax", 99.99); occscale_ = argIn.hasKey("occscale"); if (occscale_) occmax_ = argIn.getKeyDouble("occmax", 99.99); + temp = argIn.GetStringKey("adpdata"); + if (!temp.empty()) { + adpdata_ = DSLin.GetDataSet( temp ); + if (adpdata_ == 0) { + mprinterr("Error: No data set selected for 'adpdata %s'\n", temp.c_str()); + return 1; + } + if (adpdata_->Type() != DataSet::TENSOR) { + mprinterr("Error: Only TENSOR data can be used for 'adpdata'\n"); + return 1; + } + } return 0; } @@ -385,13 +400,20 @@ int Traj_PDBfile::setupTrajout(FileName const& fname, Topology* trajParm, // Set a chainID for each residue // TODO: Set different chain ID for solute mols and solvent chainID_.clear(); + // Default to a blank chain ID unless user requested PDB v3 compliance + char def_chainid; + if (pdbres_) + def_chainid = Residue::DefaultChainID(); + else + def_chainid = ' '; + // If no chain ID specified, determine chain ID. if (chainchar_ == Residue::BlankChainID()) { chainID_.reserve( trajParm->Nres() ); for (Topology::res_iterator res = trajParm->ResStart(); res != trajParm->ResEnd(); ++res) - if (dumpq_ || res->HasChainID()) + if (res->HasChainID()) chainID_.push_back( res->ChainID() ); else - chainID_.push_back( Residue::DefaultChainID() ); + chainID_.push_back( def_chainid); } else chainID_.resize(trajParm->Nres(), chainchar_); @@ -742,6 +764,7 @@ int Traj_PDBfile::writeFrame(int set, Frame const& frameOut) { if (pdbWriteMode_==MODEL) file_.WriteMODEL(set + 1); + unsigned int adpidx = 0; // Index into adpout_ float Occ = 1.0; float Bfac = 0.0; char altLoc = ' '; @@ -784,11 +807,25 @@ int Traj_PDBfile::writeFrame(int set, Frame const& frameOut) { else if (pdbTop_->Res(res).Name() == "ILE" && atomName == "CD") atomName = "CD1"; } + // TODO determine formal charges? file_.WriteCoord(rectype, anum, atomName, altLoc, resNames_[res], chainID_[res], pdbTop_->Res(res).OriginalResNum(), pdbTop_->Res(res).Icode(), Xptr[0], Xptr[1], Xptr[2], Occ, Bfac, atom.ElementName(), 0, dumpq_); + if (adpdata_ != 0 && adpidx < adpdata_->Size()) { + // Does this internal atom number match current X value? + DataSet_Tensor const& ADP = static_cast( *adpdata_ ); + unsigned int currentIdx = (unsigned int)ADP.Xvals(adpidx); + //mprintf("DEBUG: currentIdx %u aidx+1=%i\n", currentIdx, aidx+1); + if ( currentIdx == (unsigned int)(aidx + 1) ) { + DataSet_Tensor::Ttype const& UM = ADP.Tensor(adpidx); + file_.WriteANISOU( anum, atomName, resNames_[res], chainID_[res], + pdbTop_->Res(res).OriginalResNum(), + UM.Ptr(), atom.ElementName(), 0 ); + adpidx++; + } + } if (conectMode_ != NO_CONECT) atrec_[aidx] = anum; // Store ATOM record # } diff --git a/src/Traj_PDBfile.h b/src/Traj_PDBfile.h index 261b49d7d7..f340d3e1e6 100644 --- a/src/Traj_PDBfile.h +++ b/src/Traj_PDBfile.h @@ -88,6 +88,7 @@ class Traj_PDBfile: public TrajectoryIO { char chainchar_; DataSet* bfacdata_; DataSet* occdata_; + DataSet* adpdata_; ///< Hold anisotropic B-factor data for writing. double bfacmax_; double occmax_; }; diff --git a/src/Version.h b/src/Version.h index a9f19a4e52..70aabf9ce3 100644 --- a/src/Version.h +++ b/src/Version.h @@ -12,7 +12,7 @@ * Whenever a number that precedes is incremented, all subsequent * numbers should be reset to 0. */ -#define CPPTRAJ_INTERNAL_VERSION "V4.19.5" +#define CPPTRAJ_INTERNAL_VERSION "V4.20.0" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 83f099c656..a0428ce6e3 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -192,7 +192,7 @@ DataIO_VecTraj.o : DataIO_VecTraj.cpp ActionFrameCounter.h ArgList.h ArrayIterat DataIO_XVG.o : DataIO_XVG.cpp ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h BufferedLine.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_XVG.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_double.h Dimension.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h DataIO_Xplor.o : DataIO_Xplor.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h BufferedLine.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Xplor.h DataSet.h DataSetList.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h DataSet.o : DataSet.cpp ArgList.h AssociatedData.h CpptrajFile.h CpptrajStdio.h DataSet.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h TextFormat.h -DataSetList.o : DataSetList.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AtomType.h AtomTypeArray.h Box.h CharMask.h ClusterDist.h ClusterSieve.h ComplexArray.h Constants.h CoordinateInfo.h Cph.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Cmatrix.h DataSet_Cmatrix_DISK.h DataSet_Cmatrix_MEM.h DataSet_Cmatrix_NOMEM.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_Coords_TRJ.h DataSet_GridDbl.h DataSet_GridFlt.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_PHREMD.h DataSet_PHREMD_Explicit.h DataSet_PHREMD_Implicit.h DataSet_Parameters.h DataSet_RemLog.h DataSet_Topology.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_integer_disk.h DataSet_integer_mem.h DataSet_pH.h DataSet_string.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h Hungarian.h InputTrajCommon.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NC_Cmatrix.h NameType.h Parallel.h ParameterHolders.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h Vec3.h +DataSetList.o : DataSetList.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AtomType.h AtomTypeArray.h Box.h CharMask.h ClusterDist.h ClusterSieve.h ComplexArray.h Constants.h CoordinateInfo.h Cph.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Cmatrix.h DataSet_Cmatrix_DISK.h DataSet_Cmatrix_MEM.h DataSet_Cmatrix_NOMEM.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_Coords_TRJ.h DataSet_GridDbl.h DataSet_GridFlt.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_PHREMD.h DataSet_PHREMD_Explicit.h DataSet_PHREMD_Implicit.h DataSet_Parameters.h DataSet_RemLog.h DataSet_Tensor.h DataSet_Topology.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_integer_disk.h DataSet_integer_mem.h DataSet_pH.h DataSet_string.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h Hungarian.h InputTrajCommon.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NC_Cmatrix.h NameType.h Parallel.h ParameterHolders.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h SymbolExporting.h SymmetricRmsdCalc.h SymmetricTensor.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h Vec3.h DataSet_1D.o : DataSet_1D.cpp ArgList.h ArrayIterator.h AssociatedData.h ComplexArray.h Constants.h Corr.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_1D.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h PubFFT.h Range.h TextFormat.h DataSet_3D.o : DataSet_3D.cpp ArgList.h AssociatedData.h Box.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_3D.h Dimension.h FileIO.h FileName.h GridBin.h Matrix_3x3.h MetaData.h Parallel.h Range.h TextFormat.h Vec3.h DataSet_Cmatrix.o : DataSet_Cmatrix.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h Box.h CharMask.h ClusterDist.h ClusterSieve.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_1D.h DataSet_Cmatrix.h DataSet_Coords.h Dimension.h FileIO.h FileName.h Frame.h Hungarian.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Topology.h Vec3.h @@ -214,6 +214,7 @@ DataSet_PHREMD_Explicit.o : DataSet_PHREMD_Explicit.cpp ArgList.h AssociatedData DataSet_PHREMD_Implicit.o : DataSet_PHREMD_Implicit.cpp ArgList.h AssociatedData.h Cph.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_PHREMD.h DataSet_PHREMD_Implicit.h Dimension.h FileIO.h FileName.h MetaData.h NameType.h Parallel.h Range.h TextFormat.h DataSet_Parameters.o : DataSet_Parameters.cpp ArgList.h AssociatedData.h AtomType.h AtomTypeArray.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_Parameters.h Dimension.h FileIO.h FileName.h MetaData.h NameType.h Parallel.h ParameterHolders.h ParameterTypes.h Range.h TextFormat.h DataSet_RemLog.o : DataSet_RemLog.cpp ArgList.h AssociatedData.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_RemLog.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h ReplicaDimArray.h TextFormat.h +DataSet_Tensor.o : DataSet_Tensor.cpp ArgList.h AssociatedData.h CpptrajFile.h DataSet.h DataSet_Tensor.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h SymmetricTensor.h TextFormat.h DataSet_Topology.o : DataSet_Topology.cpp ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h BondSearch.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_Topology.h Dimension.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h ParmFile.h ParmIO.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Topology.h Vec3.h DataSet_Vector.o : DataSet_Vector.cpp ArgList.h ArrayIterator.h AssociatedData.h ComplexArray.h Constants.h Corr.h CpptrajFile.h DataSet.h DataSet_1D.h DataSet_Vector.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h PubFFT.h Range.h SymbolExporting.h TextFormat.h Vec3.h DataSet_double.o : DataSet_double.cpp ArgList.h AssociatedData.h CpptrajFile.h DataSet.h DataSet_1D.h DataSet_double.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h TextFormat.h @@ -359,7 +360,7 @@ Traj_GmxXtc.o : Traj_GmxXtc.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOty Traj_Gro.o : Traj_Gro.cpp Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h BufferedLine.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h Topology.h Traj_Gro.h TrajectoryIO.h Vec3.h Traj_Mol2File.o : Traj_Mol2File.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Mol2File.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Topology.h Traj_Mol2File.h TrajectoryIO.h Vec3.h Traj_NcEnsemble.o : Traj_NcEnsemble.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FrameArray.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NC_Routines.h NameType.h NetcdfFile.h Parallel.h ParallelNetcdf.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Topology.h Traj_NcEnsemble.h TrajectoryIO.h Vec3.h -Traj_PDBfile.o : Traj_PDBfile.cpp ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DistRoutines.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h PDBfile.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h Traj_PDBfile.h TrajectoryIO.h Vec3.h +Traj_PDBfile.o : Traj_PDBfile.cpp ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Tensor.h Dimension.h DistRoutines.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h PDBfile.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h SymmetricTensor.h TextFormat.h Timer.h Topology.h Traj_PDBfile.h TrajectoryIO.h Vec3.h Traj_SDF.o : Traj_SDF.cpp Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SDFfile.h SymbolExporting.h Topology.h Traj_SDF.h TrajectoryIO.h Vec3.h Traj_SQM.o : Traj_SQM.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h Topology.h Traj_SQM.h TrajectoryIO.h Vec3.h Traj_Tinker.o : Traj_Tinker.cpp Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h BufferedLine.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h TinkerFile.h Topology.h Traj_Tinker.h TrajectoryIO.h Vec3.h diff --git a/src/cpptrajfiles b/src/cpptrajfiles index a29bf7dfe6..d898082042 100644 --- a/src/cpptrajfiles +++ b/src/cpptrajfiles @@ -223,6 +223,7 @@ COMMON_SOURCES= \ DataSet_integer_disk.cpp \ DataSet_integer_mem.cpp \ DataSet_string.cpp \ + DataSet_Tensor.cpp \ Deprecated.cpp \ DihedralSearch.cpp \ DistRoutines.cpp \ diff --git a/test/Test_AtomicFluct/RunTest.sh b/test/Test_AtomicFluct/RunTest.sh index e230e80185..a1efd09830 100755 --- a/test/Test_AtomicFluct/RunTest.sh +++ b/test/Test_AtomicFluct/RunTest.sh @@ -3,7 +3,8 @@ . ../MasterTest.sh CleanFiles atomic.in fluct.*.dat dpdp.fluct.dat dpdp.adp.dat \ - fluct.2.pdb occ.2.pdb scale.2.pdb fluct.1.pdb + fluct.2.pdb occ.2.pdb scale.2.pdb fluct.1.pdb \ + dpdp.adp.pdb myfluct.adp.dat heavy.adp.pdb TESTNAME='Atomic fluctuations tests' Requires netcdf @@ -28,11 +29,20 @@ TOP=../DPDP.parm7 cat > $INPUT < $INPUT <