diff --git a/doc/CpptrajManual.pdf b/doc/CpptrajManual.pdf index 77c319d9b3..396e825409 100644 Binary files a/doc/CpptrajManual.pdf and b/doc/CpptrajManual.pdf differ diff --git a/doc/DocumentChecksums.txt b/doc/DocumentChecksums.txt index 3a73bd341d..bb2f2a27e3 100644 --- a/doc/DocumentChecksums.txt +++ b/doc/DocumentChecksums.txt @@ -1,3 +1,3 @@ f6f8cb1a79951d80a9d2656fd9c30f55 CpptrajDevelopmentGuide.lyx -30a19d148a0011fdb02ecc0d2f26e81f cpptraj.lyx +8ba86aca8b37cdcb5b01ef05141e5dce cpptraj.lyx 5d9b5b5ed47a3ded57b6464df99b3585 CpptrajManual.lyx diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index b3733716b1..a750220c87 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -20551,6 +20551,10 @@ ensemble {[] [ | last] [offset]} | lastframe bycrdidx | \end_layout +\begin_layout LyX-Code + bytemp | +\end_layout + \begin_layout LyX-Code remlog [nstlim ntwx ]}] \begin_inset Separator latexpar @@ -20628,6 +20632,10 @@ last] Frame to stop reading ensemble at; This is preferred over sorting via 'remlog'. \end_layout +\begin_layout Description +[bytemp] If specified the ensemble will be sorted by temperature values instead of replica indices. +\end_layout + \begin_layout Description [remlog \begin_inset space ~ diff --git a/src/CpptrajState.cpp b/src/CpptrajState.cpp index 1144e785d8..aba2c0565d 100644 --- a/src/CpptrajState.cpp +++ b/src/CpptrajState.cpp @@ -682,6 +682,8 @@ int CpptrajState::RunEnsemble() { EnsembleParm[0].Nframes() ); lastPindex = currentTop->Pindex(); } + // Check output trajectories against current coord info + ensembleOut_.CheckEnsembleOutCoordInfo(currentCoordInfo, FrameEnsemble); # ifdef TIMER setup_time.Stop(); # endif @@ -895,6 +897,8 @@ int CpptrajState::RunParaEnsemble() { return 1; // Allocate FrameEnsemble here in case preload is needed. FrameEnsemble.SetupFrames( NAV.FirstParm()->Atoms(), NAV.EnsCoordInfo() ); + // Check output trajectories against current coord info + ensembleOut_.CheckEnsembleOutCoordInfo(NAV.EnsCoordInfo(), FrameEnsemble); // Figure out if any frames need to be preloaded on ranks int preload_err = 0; if (!TrajComm.Master()) { @@ -931,6 +935,13 @@ int CpptrajState::RunParaEnsemble() { break; } if (!NAV.CurrentEns()->BadEnsemble()) { + // Has the ensemble changed? + if (set != my_start && NAV.EnsembleHasChanged()) { + if (debug_ > 0) + rprintf("DEBUG: Ensemble has changed at set %i\n", set); + // Check output trajectories against current input ensemble + ensembleOut_.CheckEnsembleOutCoordInfo(NAV.CurrentEns()->EnsembleCoordInfo(), FrameEnsemble); + } // Since Frame can be modified by actions, save original and use currentFrame ActionFrame currentFrame( SortedFrames[0], set ); if ( currentFrame.Frm().CheckCoordsInvalid() ) @@ -1337,6 +1348,8 @@ int CpptrajState::RunNormal() { currentSetup.Nframes() ); lastPindex = currentSetup.Top().Pindex(); } + // Check output trajectories against current coord info + trajoutList_.CheckTrajoutCoordInfo(currentSetup.CoordInfo(), TrajFrame); # ifdef TIMER setup_time.Stop(); # endif diff --git a/src/DataSet_Coords.cpp b/src/DataSet_Coords.cpp index ed097df7ac..64146e7f19 100644 --- a/src/DataSet_Coords.cpp +++ b/src/DataSet_Coords.cpp @@ -10,6 +10,7 @@ Frame DataSet_Coords::AllocateFrame() const { void DataSet_Coords::CommonInfo() const { if (cInfo_.HasBox()) mprintf(" Box Coords,"); if (cInfo_.HasVel()) mprintf(" Velocities,"); + if (cInfo_.HasForce()) mprintf(" Forces,"); mprintf(" %i atoms", top_.Natom()); } diff --git a/src/DataSet_Coords_TRJ.cpp b/src/DataSet_Coords_TRJ.cpp index 6b43eeb7e1..7b87947b45 100644 --- a/src/DataSet_Coords_TRJ.cpp +++ b/src/DataSet_Coords_TRJ.cpp @@ -111,8 +111,18 @@ void DataSet_Coords_TRJ::GetFrame(int idx, Frame& fIn) { Traj_ = trajinList_[ IDX_.CurrentTrajNum() ]; // NOTE: Currently enforcing all traj have same # atoms, no need to check topology. // See if frame needs (re-)allocation. - if (previousTraj == 0 || previousTraj->TrajCoordInfo() != Traj_->TrajCoordInfo()) + if (previousTraj == 0 || previousTraj->TrajCoordInfo() != Traj_->TrajCoordInfo()) { readFrame_.SetupFrameV( Traj_->Traj().Parm()->Atoms(), Traj_->TrajCoordInfo() ); + // If previously had vel/frc info but now do not, zero it out. + if (CoordsInfo().HasVel() && !Traj_->TrajCoordInfo().HasVel()) + readFrame_.AddVelocities(Frame::Darray(readFrame_.size(), 0.0)); + if (fIn.HasVelocity() && !Traj_->TrajCoordInfo().HasVel()) + fIn.ZeroVelocities(); + if (CoordsInfo().HasForce() && !Traj_->TrajCoordInfo().HasForce()) + readFrame_.AddForces(Frame::Darray(readFrame_.size(), 0.0)); + if (fIn.HasForce() && !Traj_->TrajCoordInfo().HasForce()) + fIn.ZeroForces(); + } // Open traj. if (Traj_->BeginTraj()) { mprinterr("Error: Could not open trajectory %i '%s'\n", IDX_.CurrentTrajNum(), diff --git a/src/EnsembleIn.cpp b/src/EnsembleIn.cpp index 5d43debfff..d27ac2eab7 100644 --- a/src/EnsembleIn.cpp +++ b/src/EnsembleIn.cpp @@ -1,5 +1,6 @@ #include "EnsembleIn.h" #include "CpptrajStdio.h" +#include "Constants.h" // FNE EnsembleIn::EnsembleIn() : targetType_(ReplicaInfo::NONE), badEnsemble_(0), debug_(0) @@ -35,6 +36,27 @@ int EnsembleIn::GatherIndices(int* iAddress, std::vector& allIndice return 0; } +int EnsembleIn::GatherValues(Frame const& frameIn, ReplicaDimArray const& replicaDims, + std::vector& allValues, + Parallel::Comm const& commIn) +{ + int Ndims = replicaDims.Ndims(); + std::vector all_values(allValues.size() * Ndims, 0); + std::vector dtemp(Ndims); + for (int idim = 0; idim != replicaDims.Ndims(); idim++) + dtemp[idim] = frameIn.RemdValue(idim, replicaDims); + if (commIn.AllGather( &dtemp[0], Ndims, MPI_DOUBLE, &all_values[0])) { + rprinterr("Error: Gathering replica values.\n"); + return 1; + } + std::vector::const_iterator idx_it = all_values.begin(); + for (std::vector::iterator it = allValues.begin(); + it != allValues.end(); + ++it, idx_it += Ndims) + it->assign(idx_it, idx_it + Ndims); + return 0; +} + #ifdef TIMER double EnsembleIn::total_mpi_allgather_ = 0.0; double EnsembleIn::total_mpi_sendrecv_ = 0.0; @@ -93,3 +115,92 @@ int EnsembleIn::SetIndicesMap(std::vector const& allIndices) { } return 0; } + +int EnsembleIn::SetIdxValMap(std::vector const& allIndices, + std::vector const& allValues) +{ + if (allIndices.size() != allValues.size()) { + rprinterr("Internal Error: EnsembleIn::SetIdxValMap(): allIndices size %zu != allValues size %zu\n", + allIndices.size(), allValues.size()); + return 1; + } + + if (debug_ > 0) { + mprintf("DEBUG: index/value map:\n"); + for (unsigned int member = 0; member != allIndices.size(); member++) { + mprintf("\t"); + for (Frame::RemdIdxType::const_iterator idx = allIndices[member].begin(); + idx != allIndices[member].end(); ++idx) + mprintf(" %6i", *idx); + mprintf(" : "); + for (Darray::const_iterator dval = allValues[member].begin(); + dval != allValues[member].end(); ++dval) + mprintf(" %8.3f", *dval); + mprintf("\n"); + } + } + + typedef std::pair ivPair; + for (unsigned int idx = 0; idx != allIndices.size(); idx++) + RemdIdxValMap_.insert( ivPair(allIndices[idx], allValues[idx]) ); + return 0; +} + +static inline void indexprint(Frame::RemdIdxType const& IDXin) { + for (Frame::RemdIdxType::const_iterator idx = IDXin.begin(); + idx != IDXin.end(); ++idx) + mprintf(" %6i", *idx); +} + +static inline void valprint(std::vector const& VALin) { + for (std::vector::const_iterator dval = VALin.begin(); + dval != VALin.end(); ++dval) + mprintf(" %8.3f", *dval); +} + +/// Floating point not equals. +static inline bool FNE(double v1, double v2) { + double delta = v1 - v2; + if (delta < 0.0) delta = -delta; + return (delta > Constants::SMALL); +} + +// Compare Darrays +static inline bool DNE(std::vector const& lhs, + std::vector const& rhs) +{ + if (lhs.size() != rhs.size()) return true; + for (unsigned int idx = 0; idx < lhs.size(); idx++) + if (FNE(lhs[idx], rhs[idx])) return true; + return false; +} + +int EnsembleIn::CheckIdxValMap(RemdIdxValMapType const& mapIn) const { + if (mapIn.size() != RemdIdxValMap_.size()) { + mprintf("Warning: REMD index/value map size %zu does not match previous map size %zu\n", + RemdIdxValMap_.size(), mapIn.size()); + return 1; + } + int match = 0; + for (RemdIdxValMapType::const_iterator it = RemdIdxValMap_.begin(); + it != RemdIdxValMap_.end(); ++it) + { + RemdIdxValMapType::const_iterator jt = mapIn.find( it->first ); + if (jt == mapIn.end()) { + mprintf("Warning: Index"); + indexprint( it->first ); + mprintf(" for this ensemble not found in previous ensemble.\n"); + match = 1; + } else if (DNE(jt->second, it->second)) { + mprintf("Warning: Index"); + indexprint( it->first ); + mprintf(" values "); + valprint( it->second ); + mprintf(" do not match previous map values "); + valprint( jt->second ); + mprintf("\n"); + match = 1; + } + } + return match; +} diff --git a/src/EnsembleIn.h b/src/EnsembleIn.h index 7cca6da76a..11d8b99d4d 100644 --- a/src/EnsembleIn.h +++ b/src/EnsembleIn.h @@ -16,7 +16,12 @@ class FrameArray; class CoordinateInfo; /// Read in an array of frames at a time. class EnsembleIn { + protected: + typedef Frame::RemdIdxType RemdIdxType; + typedef std::vector Darray; public: + typedef std::map RemdIdxValMapType; + EnsembleIn(); virtual ~EnsembleIn() {} virtual int SetupEnsembleRead(FileName const&, ArgList&, Topology*) = 0; @@ -36,19 +41,28 @@ class EnsembleIn { InputTrajCommon const& Traj() const { return traj_; } /// \return true if there was a problem reading the ensemble. int BadEnsemble() const { return badEnsemble_; } + /// \return 1 if this index/value map does not match the given map + int CheckIdxValMap(RemdIdxValMapType const&) const; + /// \return index/value map + RemdIdxValMapType const& RemdIdxValMap() const { return RemdIdxValMap_; } + /// \return The sorting mode + ReplicaInfo::TargetType TargetMode() const { return targetType_; } + /// Set the debug level void SetDebug(int d) { debug_ = d; } /// Write current replica mapping to STDOUT. void PrintReplicaInfo() const; protected: - typedef Frame::RemdIdxType RemdIdxType; int SetTemperatureMap(std::vector const&); int SetIndicesMap(std::vector const&); + int SetIdxValMap(std::vector const&, std::vector const&); InputTrajCommon& SetTraj() { return traj_; } /// For converting temperature to replica index ReplicaInfo::Map TemperatureMap_; /// For converting indices to replica index ReplicaInfo::Map IndicesMap_; + + RemdIdxValMapType RemdIdxValMap_; ReplicaInfo::TargetType targetType_; ///< Hold type of REMD frame being searched for. int badEnsemble_; ///< Set to 1 if problem reading ensemble, 0 otherwise. int debug_; @@ -60,6 +74,8 @@ class EnsembleIn { Parallel::Comm const& EnsembleComm() const { return Parallel::EnsembleComm(); } static int GatherTemperatures(double*, std::vector&, Parallel::Comm const&); static int GatherIndices(int*, std::vector&, int, Parallel::Comm const&); + static int GatherValues(Frame const&, ReplicaDimArray const&, + std::vector&, Parallel::Comm const&); # ifdef TIMER Timer mpi_allgather_timer_; Timer mpi_sendrecv_timer_; @@ -68,8 +84,8 @@ class EnsembleIn { # endif # endif private: - static void PrintReplicaTmap(ReplicaInfo::Map const&); - static void PrintReplicaImap(ReplicaInfo::Map const&); + //static void PrintReplicaTmap(ReplicaInfo::Map const&); + //static void PrintReplicaImap(ReplicaInfo::Map const&); InputTrajCommon traj_; # ifdef MPI diff --git a/src/EnsembleIn_Multi.cpp b/src/EnsembleIn_Multi.cpp index feff51061b..b2af096e7f 100644 --- a/src/EnsembleIn_Multi.cpp +++ b/src/EnsembleIn_Multi.cpp @@ -23,6 +23,7 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn, double remlog_ntwx = argIn.getKeyDouble("ntwx", 1.0); bool no_sort = argIn.hasKey("nosort"); bool by_crdidx = argIn.hasKey("bycrdidx"); + bool by_temp = argIn.hasKey("bytemp"); if (no_sort && by_crdidx) { mprinterr("Error: Should only specify either 'nosort' or 'bycrdidx', not both.\n"); return 1; @@ -35,6 +36,12 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn, mprinterr("Error: Should only specify either 'bycrdidx' or 'remlog', not both.\n"); return 1; } + if (by_temp) { + if (by_crdidx || no_sort) { + mprinterr("Error: 'bycrdidx'/'nosort' should not be specified with 'bytemp'\n"); + return 1; + } + } // CRDIDXARG: Parse out 'crdidx ' now so it is not processed // by SetupTrajIO. ArgList crdidxarg; @@ -127,8 +134,14 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn, } } else { // If dimensions are present index by replica indices, otherwise index - // by temperature. - if (cInfo_.ReplicaDimensions().Ndims() > 0) + // by temperature. + if (by_temp) { + targetType_ = ReplicaInfo::TEMP; + if (cInfo_.ReplicaDimensions().Ndims() > 1) { + mprinterr("Error: 'bytemp' not compatible with multiple dimensions.\n"); + return 1; + } + } else if (cInfo_.ReplicaDimensions().Ndims() > 0) targetType_ = ReplicaInfo::INDICES; else targetType_ = ReplicaInfo::TEMP; @@ -147,10 +160,14 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn, frameIn.SetupFrameV( Traj().Parm()->Atoms(), cInfo_ ); std::vector allTemps; std::vector allIndices; + std::vector allValues; if (targetType_ == ReplicaInfo::TEMP) allTemps.assign(REMDtraj_.size(), -1.0); - else if (targetType_ == ReplicaInfo::INDICES) + else if (targetType_ == ReplicaInfo::INDICES) { allIndices.resize( REMDtraj_.size() ); + if (cInfo_.UseRemdValues()) + allValues.resize( REMDtraj_.size(), Darray(cInfo_.ReplicaDimensions().Ndims(), -1.0) ); + } # ifdef MPI int err = 0; if (Parallel::TrajComm().Master()) { @@ -174,6 +191,10 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn, if (GatherIndices(frameIn.iAddress(), allIndices, cInfo_.ReplicaDimensions().Ndims(), EnsembleComm())) err = 4; + if (cInfo_.UseRemdValues()) { + if (GatherValues(frameIn, cInfo_.ReplicaDimensions(), allValues, EnsembleComm())) + err = 5; + } } } } @@ -196,6 +217,18 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn, Parallel::TrajComm().MasterBcast( &((*it)[0]), it->size(), MPI_INT ); } } + if (cInfo_.UseRemdValues()) { + for (std::vector::iterator it = allValues.begin(); + it != allValues.end(); ++it) + { + if (Parallel::TrajComm().Master()) + Parallel::TrajComm().MasterBcast( &((*it)[0]), it->size(), MPI_DOUBLE ); + else { + it->resize( cInfo_.ReplicaDimensions().Ndims() ); + Parallel::TrajComm().MasterBcast( &((*it)[0]), it->size(), MPI_DOUBLE ); + } + } + } } } # else /* MPI */ @@ -205,14 +238,24 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn, REMDtraj_[member]->closeTraj(); if (targetType_ == ReplicaInfo::TEMP) allTemps[member] = frameIn.Temperature(); - else if (targetType_ == ReplicaInfo::INDICES) + else if (targetType_ == ReplicaInfo::INDICES) { allIndices[member] = frameIn.RemdIndices(); + if (cInfo_.UseRemdValues()) { + for (int idim = 0; idim != cInfo_.ReplicaDimensions().Ndims(); idim++) { + double dval = frameIn.RemdValue(idim, cInfo_.ReplicaDimensions()); + allValues[member][idim] = dval; + } + } // END if UseRemdValues + } } # endif /* MPI */ if (targetType_ == ReplicaInfo::TEMP) { if (SetTemperatureMap(allTemps)) return 1; } else if (targetType_ == ReplicaInfo::INDICES) { if (SetIndicesMap(allIndices)) return 1; + if (cInfo_.UseRemdValues()) { + if (SetIdxValMap(allIndices, allValues)) return 1; + } } } // Otherwise NONE, no sorting diff --git a/src/EnsembleIn_Multi.h b/src/EnsembleIn_Multi.h index dc19fc90a7..ef1ac4663f 100644 --- a/src/EnsembleIn_Multi.h +++ b/src/EnsembleIn_Multi.h @@ -20,7 +20,6 @@ class EnsembleIn_Multi : public EnsembleIn { std::string FinalCrdIndices() const; // TODO Remove /// \return true if sorting crdidx by remlog bool SortingByRemlog() const; - ReplicaInfo::TargetType TargetMode() const { return targetType_; } // TODO Remove private: TrajIOarray REMDtraj_; CoordinateInfo cInfo_; diff --git a/src/EnsembleNavigator.h b/src/EnsembleNavigator.h index 4377265076..e90167be22 100644 --- a/src/EnsembleNavigator.h +++ b/src/EnsembleNavigator.h @@ -9,13 +9,14 @@ */ class EnsembleNavigator { public: - EnsembleNavigator() : currentEns_(0), FirstParm_(0) {} + EnsembleNavigator() : currentEns_(0), FirstParm_(0), ensembleHasChanged_(false) {} int AddEnsembles(TrajinList::ensemble_it const&, TrajinList::ensemble_it const&); inline int GetEnsemble(int, FrameArray&, FramePtrArray&); EnsembleIn* CurrentEns() { return currentEns_; } Topology* FirstParm() { return FirstParm_; } CoordinateInfo const& EnsCoordInfo() const { return ensCoordInfo_; } TrajFrameIndex const& IDX() const { return IDX_; } + bool EnsembleHasChanged() const { return ensembleHasChanged_; } private: typedef std::vector Earray; Earray Ensembles_; ///< Array of input ensembles @@ -23,14 +24,17 @@ class EnsembleNavigator { TrajFrameIndex IDX_; ///< Used to convert global index to individual index EnsembleIn* currentEns_; ///< Currently open ensemble Topology* FirstParm_; ///< Topology associated with first ensemble + bool ensembleHasChanged_; ///< Set to true if the last call to GetEnsemble() was the start of another ensemble. }; /// Get ensemble set specified by global set index int EnsembleNavigator::GetEnsemble(int set, FrameArray& FrameEnsemble, FramePtrArray& SortedFrames) { + ensembleHasChanged_ = false; int internalIdx = IDX_.FindIndex( set ); // If desired ensemble is different than current, open desired ensemble. if (IDX_.TrajHasChanged()) { + ensembleHasChanged_ = true; if (currentEns_ != 0) currentEns_->EndEnsemble(); currentEns_ = Ensembles_[ IDX_.CurrentTrajNum() ]; // NOTE: Need to check for reallocation? TODO better error check diff --git a/src/EnsembleOutList.cpp b/src/EnsembleOutList.cpp index 371c707ccd..2306e8afab 100644 --- a/src/EnsembleOutList.cpp +++ b/src/EnsembleOutList.cpp @@ -4,6 +4,7 @@ #include "EnsembleOut_Single.h" #include "EnsembleOut_Multi.h" #include "ArgList.h" +#include "FrameArray.h" /// CONSTRUCTOR EnsembleOutList::EnsembleOutList() : debug_(0) {} @@ -90,6 +91,51 @@ int EnsembleOutList::SetupEnsembleOut(Topology* CurrentParm, CoordinateInfo cons return 0; } +/** Check active output ensembles against the current coordinate info. */ +int EnsembleOutList::CheckEnsembleOutCoordInfo(CoordinateInfo const& currentCoordInfo, + FrameArray& FrameEnsemble) +const +{ + for (EnsArray::const_iterator it = active_.begin(); it != active_.end(); ++it) + { + EnsembleOut* activeTraj = *it; + CoordinateInfo const& ensembleoutCoordInfo = activeTraj->Traj().CoordInfo(); + // Velocities + if (currentCoordInfo.HasVel() != ensembleoutCoordInfo.HasVel()) { + if (currentCoordInfo.HasVel()) { + mprintf("Warning: Input ensemble has velocity information but output ensemble was set up without velocities.\n"); + mprintf("Warning: Velocity information will not be written to output ensemble.\n"); + } else { + mprintf("Warning: Input ensemble has no velocity information but output ensemble was set up with velocities.\n"); + mprintf("Warning: All zeroes will be written for velocities.\n"); + for (unsigned int idx = 0; idx < FrameEnsemble.Size(); idx++) { + if (FrameEnsemble[idx].HasVelocity()) + FrameEnsemble[idx].ZeroVelocities(); + else + FrameEnsemble[idx].AddVelocities(Frame::Darray(FrameEnsemble[idx].size(), 0.0)); + } + } + } + // Forces + if (currentCoordInfo.HasForce() != ensembleoutCoordInfo.HasForce()) { + if (currentCoordInfo.HasForce()) { + mprintf("Warning: Input ensemble has force information but output ensemble was set up without forces.\n"); + mprintf("Warning: Force information will not be written to output ensemble.\n"); + } else { + mprintf("Warning: Input ensemble has no force information but output ensemble was set up with forces.\n"); + mprintf("Warning: All zeroes will be written for forces.\n"); + for (unsigned int idx = 0; idx < FrameEnsemble.Size(); idx++) { + if (FrameEnsemble[idx].HasForce()) + FrameEnsemble[idx].ZeroForces(); + else + FrameEnsemble[idx].AddForces(Frame::Darray(FrameEnsemble[idx].size(), 0.0)); + } + } + } + } + return 0; +} + /** List only active output trajectories. */ void EnsembleOutList::ListActive() const { if (!ensout_.empty()) { diff --git a/src/EnsembleOutList.h b/src/EnsembleOutList.h index dcafea62eb..0bf5f7e2ac 100644 --- a/src/EnsembleOutList.h +++ b/src/EnsembleOutList.h @@ -13,6 +13,7 @@ class ArgList; class DataSetList; class Topology; class CoordinateInfo; +class FrameArray; /// Hold output ensembles. class EnsembleOutList { public: @@ -24,6 +25,8 @@ class EnsembleOutList { int AddEnsembleOut(std::string const&, ArgList const&, DataSetList const&, Topology*, int); /// Set up ensembles for given topology. int SetupEnsembleOut(Topology*, CoordinateInfo const&, int); //TODO const Topology? + /// Check active ensembles against given CoordinateInfo + int CheckEnsembleOutCoordInfo(CoordinateInfo const&, FrameArray&) const; /// Write array of Frames to ensemble int WriteEnsembleOut(int, FramePtrArray const&); /// Close all ensembles diff --git a/src/Exec_Help.cpp b/src/Exec_Help.cpp index 635b557b26..23fdcb5f2a 100644 --- a/src/Exec_Help.cpp +++ b/src/Exec_Help.cpp @@ -149,6 +149,9 @@ int Exec_Help::Math(ArgList& argIn) const { "\tsin() : Sine\n" "\tcos() : Cosine\n" "\ttan() : Tangent\n" + "\tsinh() : Hyperbolic Sine\n" + "\tcosh() : Hyperbolic Cosine\n" + "\ttanh() : Hyperbolic Tangent\n" " Data set functions\n" "\tsum() : Summation\n" "\tavg() : Average\n" diff --git a/src/Exec_Traj.cpp b/src/Exec_Traj.cpp index 55cc32df02..1b5beeec07 100644 --- a/src/Exec_Traj.cpp +++ b/src/Exec_Traj.cpp @@ -37,12 +37,17 @@ void Exec_Ensemble::Help() const { mprintf("\t [trajnames ,,...,]\n" "\t [{nosort |\n" "\t bycrdidx |\n" + "\t bytemp |\n" "\t remlog [nstlim ntwx ]}]\n" " Load an ensemble of trajectories starting with that will be\n" " processed together as an ensemble.\n" " The default behavior is to sort the ensemble by replica. If 'nosort'\n" " is specified the incoming ensemble will not be sorted; this is useful\n" " if the ensemble is already sorted or the trajectories are independent.\n" + " If 'bycrdidx' is specified the ensemble will be sorted by coordinate index\n" + " if it is present in the incoming trajectories.\n" + " If 'bytemp' is specified the ensemble will be sorted by temperature values\n" + " instead of replica indices.\n" " If 'remlog' is specified the ensemble will be sorted by coordinate index;\n" " 'nstlim' specifies the number of steps between exchanges and 'ntwx'\n" " specifies the number of steps between trajectory writes.\n" diff --git a/src/Frame.cpp b/src/Frame.cpp index fe1207a06f..265f69dc7c 100644 --- a/src/Frame.cpp +++ b/src/Frame.cpp @@ -262,6 +262,21 @@ Frame &Frame::operator=(Frame rhs) { } // ---------- ACCESS INTERNAL DATA --------------------------------------------- +/** \return Replica value for specified dimension */ +double Frame::RemdValue(int idim, ReplicaDimArray const& repDim) const { + if (remd_indices_.empty()) return -1.0; + double dval = -1.0; + switch (repDim.DimType(idim)) { + case ReplicaDimArray::RXSGLD : // FIXME double check this + case ReplicaDimArray::TEMPERATURE : dval = Temperature(); break; + case ReplicaDimArray::PH : dval = pH(); break; + case ReplicaDimArray::REDOX : dval = RedOx(); break; + // Default to replica index + default : dval = (double)RemdIndices()[idim]; break; + } + return dval; +} + /** \return CoordinateInfo describing the Frame. */ CoordinateInfo Frame::CoordsInfo() const { // TODO no good way to tell about time yet. diff --git a/src/Frame.h b/src/Frame.h index d45d0221ad..7fb807db01 100644 --- a/src/Frame.h +++ b/src/Frame.h @@ -78,6 +78,8 @@ class Frame { double pH() const { return pH_; } double RedOx() const { return redox_; } double Time() const { return time_; } + /// \return Replica value for specified dimension. + double RemdValue(int, ReplicaDimArray const&) const; /// \return CoordinateInfo that describes the Frame CoordinateInfo CoordsInfo() const; /// \return pointer to start of XYZ coords for given atom. diff --git a/src/RPNcalc.cpp b/src/RPNcalc.cpp index 6af5bf28f7..9d68305949 100644 --- a/src/RPNcalc.cpp +++ b/src/RPNcalc.cpp @@ -62,6 +62,9 @@ const RPNcalc::FnIdType RPNcalc::FnIdArray_[] = { { FN_SIN, 4, "sin(" }, { FN_COS, 4, "cos(" }, { FN_TAN, 4, "tan(" }, + { FN_SINH, 5, "sinh(" }, + { FN_COSH, 5, "cosh(" }, + { FN_TANH, 5, "tanh(" }, { FN_SUM, 4, "sum(" }, { FN_AVG, 4, "avg(" }, { FN_STDEV, 6, "stdev(" }, @@ -372,6 +375,9 @@ double RPNcalc::DoOperation(double d1, double d2, TokenType op_type) { case FN_SIN: return sin(d1); case FN_COS: return cos(d1); case FN_TAN: return tan(d1); + case FN_SINH: return sinh(d1); + case FN_COSH: return cosh(d1); + case FN_TANH: return tanh(d1); case FN_STDEV: return 0.0; case FN_SUM: case FN_AVG: @@ -1142,6 +1148,9 @@ const RPNcalc::OpType RPNcalc::Token::OpArray_[] = { { 0, 1, NO_A, FN, 0, "Sine" }, // FN_SIN { 0, 1, NO_A, FN, 0, "Cosine" }, // FN_COS { 0, 1, NO_A, FN, 0, "Tangent" }, // FN_TAN + { 0, 1, NO_A, FN, 0, "Hyperbolic Sine" }, // FN_SINH + { 0, 1, NO_A, FN, 0, "Hyperbolic Cosine" }, // FN_COSH + { 0, 1, NO_A, FN, 0, "Hyperbolic Tangent" }, // FN_TANH { 0, 1, NO_A, FN, 1, "Summation" }, // FN_SUM { 0, 1, NO_A, FN, 1, "Average" }, // FN_AVG { 0, 1, NO_A, FN, 1, "Standard Dev"}, // FN_STDEV diff --git a/src/RPNcalc.h b/src/RPNcalc.h index b57bd147ee..aebc15f00e 100644 --- a/src/RPNcalc.h +++ b/src/RPNcalc.h @@ -38,6 +38,8 @@ class RPNcalc { FN_SQRT, FN_EXP, FN_LN, FN_ABS, // Trig functions FN_SIN, FN_COS, FN_TAN, + // Hyperbolic Trig Functions + FN_SINH, FN_COSH, FN_TANH, // Functions that take a data set FN_SUM, FN_AVG, FN_STDEV, FN_MIN, FN_MAX, // Parentheses (for infix conversion only) diff --git a/src/TrajinList.cpp b/src/TrajinList.cpp index cabcc56f1b..790688d2b4 100644 --- a/src/TrajinList.cpp +++ b/src/TrajinList.cpp @@ -87,6 +87,16 @@ int TrajinList::AddEnsembleIn(std::string const& fname, Topology* topIn, ArgList err++; continue; } + // Check that replica indices match to the same values as previous ensemble + if (ensemble->EnsembleCoordInfo().UseRemdValues() && !ensemble_.empty()) { + if (ensemble->CheckIdxValMap( ensemble_.back()->RemdIdxValMap() )) { + if (ensemble->TargetMode() == ReplicaInfo::INDICES) { + mprinterr("Error: Ensemble replica indices do not map to the same values as previous ensemble.\n" + "Error: Use the 'bytemp' keyword for this ensemble.\n"); + return 1; + } + } + } if (args.CheckForMoreArgs()) { delete ensemble; return 1; diff --git a/src/TrajoutList.cpp b/src/TrajoutList.cpp index 671f1a88b3..f5dff1b551 100644 --- a/src/TrajoutList.cpp +++ b/src/TrajoutList.cpp @@ -96,6 +96,47 @@ int TrajoutList::SetupTrajout(Topology* CurrentParm, CoordinateInfo const& cInfo return 0; } +/** Check active output trajectories against the current input coordinate info. */ +int TrajoutList::CheckTrajoutCoordInfo(CoordinateInfo const& currentCoordInfo, + Frame& currentFrame) +const +{ + for (ListType::const_iterator it = active_.begin(); it != active_.end(); ++it) + { + Trajout_Single* activeTraj = *it; + CoordinateInfo const& trajoutCoordInfo = activeTraj->Traj().CoordInfo(); + // Velocities + if (currentCoordInfo.HasVel() != trajoutCoordInfo.HasVel()) { + if (currentCoordInfo.HasVel()) { + mprintf("Warning: Input trajectory has velocity information but output trajectory was set up without velocities.\n"); + mprintf("Warning: Velocity information will not be written to output trajectory.\n"); + } else { + mprintf("Warning: Input trajectory has no velocity information but output trajectory was set up with velocities.\n"); + mprintf("Warning: All zeroes will be written for velocities.\n"); + if (currentFrame.HasVelocity()) + currentFrame.ZeroVelocities(); + else + currentFrame.AddVelocities(Frame::Darray(currentFrame.size(), 0.0)); + } + } + // Forces + if (currentCoordInfo.HasForce() != trajoutCoordInfo.HasForce()) { + if (currentCoordInfo.HasForce()) { + mprintf("Warning: Input trajectory has force information but output trajectory was set up without forces.\n"); + mprintf("Warning: Force information will not be written to output trajectory.\n"); + } else { + mprintf("Warning: Input trajectory has no force information but output trajectory was set up with forces.\n"); + mprintf("Warning: All zeroes will be written for forces.\n"); + if (currentFrame.HasForce()) + currentFrame.ZeroForces(); + else + currentFrame.AddForces(Frame::Darray(currentFrame.size(), 0.0)); + } + } + } + return 0; +} + /** List only active output trajectories. */ void TrajoutList::ListActive() const { if (!trajout_.empty()) { diff --git a/src/TrajoutList.h b/src/TrajoutList.h index 292e3ba4d2..6a75bc2c96 100644 --- a/src/TrajoutList.h +++ b/src/TrajoutList.h @@ -25,6 +25,8 @@ class TrajoutList { int AddTrajout(std::string const&, ArgList const&, DataSetList const&, Topology*); /// Set up trajectories for given topology. int SetupTrajout(Topology*, CoordinateInfo const&, int); + /// Check active trajectories against given CoordinateInfo + int CheckTrajoutCoordInfo(CoordinateInfo const&, Frame&) const; /// Write frame to normal output trajectories. int WriteTrajout(int, Frame const&); /// Call end for all trajectories diff --git a/src/Trajout_Single.h b/src/Trajout_Single.h index e1988f2d76..10292acf2d 100644 --- a/src/Trajout_Single.h +++ b/src/Trajout_Single.h @@ -33,7 +33,7 @@ class Trajout_Single { /// Print information on trajectory to be written. void PrintInfo(int) const; // ------------------------------------------- - OutputTrajCommon Traj() const { return traj_; } + OutputTrajCommon const& Traj() const { return traj_; } bool IsInitialized() const { return trajio_ != 0; } /// Init and setup/open traj. int PrepareTrajWrite(FileName const&, ArgList const&, DataSetList const&, Topology*, diff --git a/src/Version.h b/src/Version.h index b58fa75374..de05d7979a 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 "V7.9.1" +#define CPPTRAJ_INTERNAL_VERSION "V7.10.0" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/src/cpptrajdepend b/src/cpptrajdepend index b1256246a6..0ec3c0e47b 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -307,12 +307,12 @@ Energy/PME_RecipParams.o : Energy/PME_RecipParams.cpp Box.h CpptrajStdio.h Energ Energy/VDW_LongRange_Correction.o : Energy/VDW_LongRange_Correction.cpp Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h Energy/VDW_LongRange_Correction.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h Unit.h Vec3.h EnergyArray.o : EnergyArray.cpp CpptrajFile.h CpptrajStdio.h EnergyArray.h FileIO.h FileName.h Parallel.h Energy_Sander.o : Energy_Sander.cpp ArgList.h Atom.h AtomMask.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h Energy_Sander.h FileName.h FileTypes.h File_TempName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h ParmFile.h Range.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h Topology.h Unit.h Vec3.h -EnsembleIn.o : EnsembleIn.cpp Atom.h AtomMask.h Box.h CoordinateInfo.h CpptrajStdio.h EnsembleIn.h FileName.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h Timer.h TrajFrameCounter.h Unit.h Vec3.h +EnsembleIn.o : EnsembleIn.cpp Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h EnsembleIn.h FileName.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h Timer.h TrajFrameCounter.h Unit.h Vec3.h EnsembleIn_Multi.o : EnsembleIn_Multi.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_RemLog.h Dimension.h EnsembleIn.h EnsembleIn_Multi.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajIOarray.h TrajectoryIO.h Unit.h Vec3.h EnsembleIn_Single.o : EnsembleIn_Single.cpp ArgList.h Atom.h AtomMask.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h EnsembleIn.h EnsembleIn_Single.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Unit.h Vec3.h EnsembleNavigator.o : EnsembleNavigator.cpp Atom.h AtomMask.h Box.h CoordinateInfo.h CpptrajStdio.h EnsembleIn.h EnsembleNavigator.h FileName.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h Timer.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h TrajinList.h Unit.h Vec3.h EnsembleOut.o : EnsembleOut.cpp ActionFrameCounter.h BaseIOtype.h Box.h CoordinateInfo.h CpptrajStdio.h EnsembleOut.h FileName.h FileTypes.h FramePtrArray.h Matrix_3x3.h OutputTrajCommon.h Parallel.h Range.h ReplicaDimArray.h TrajectoryFile.h Vec3.h -EnsembleOutList.o : EnsembleOutList.cpp ActionFrameCounter.h ArgList.h Atom.h AtomMask.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h EnsembleOut.h EnsembleOutList.h EnsembleOut_Multi.h EnsembleOut_Single.h FileName.h FileTypes.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TrajectoryFile.h Unit.h Vec3.h +EnsembleOutList.o : EnsembleOutList.cpp ActionFrameCounter.h ArgList.h Atom.h AtomMask.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h EnsembleOut.h EnsembleOutList.h EnsembleOut_Multi.h EnsembleOut_Single.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TrajectoryFile.h Unit.h Vec3.h EnsembleOut_Multi.o : EnsembleOut_Multi.cpp ActionFrameCounter.h ArgList.h Atom.h AtomMask.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h EnsembleOut.h EnsembleOut_Multi.h FileName.h FileTypes.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h Topology.h TrajectoryFile.h TrajectoryIO.h Unit.h Vec3.h EnsembleOut_Single.o : EnsembleOut_Single.cpp ActionFrameCounter.h ArgList.h Atom.h AtomMask.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h EnsembleOut.h EnsembleOut_Single.h FileName.h FileTypes.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TrajectoryFile.h TrajectoryIO.h Unit.h Vec3.h EwaldOptions.o : EwaldOptions.cpp ArgList.h CpptrajStdio.h EwaldOptions.h diff --git a/test/Test_VelFrc/RST/001.rst7 b/test/Test_VelFrc/RST/001.rst7 new file mode 100644 index 0000000000..b28e2a8ca9 Binary files /dev/null and b/test/Test_VelFrc/RST/001.rst7 differ diff --git a/test/Test_VelFrc/RST/002.rst7 b/test/Test_VelFrc/RST/002.rst7 new file mode 100644 index 0000000000..149042c2a2 Binary files /dev/null and b/test/Test_VelFrc/RST/002.rst7 differ diff --git a/test/Test_VelFrc/RunTest.sh b/test/Test_VelFrc/RunTest.sh index 5184cf5ae4..890c070209 100755 --- a/test/Test_VelFrc/RunTest.sh +++ b/test/Test_VelFrc/RunTest.sh @@ -3,70 +3,89 @@ . ../MasterTest.sh CleanFiles cpptraj.in CrdFrcVel.nc Vel.crd Frc.crd Vel1.crd Frc1.crd \ - CrdFrcVel.ncrst.? fromncrst.nc trpzip2.*.crd + CrdFrcVel.ncrst.? fromncrst.nc trpzip2.*.crd combined.nc \ + combined2.nc ensemble.nc.? TESTNAME='Read separate velocity/force trajectory data tests' -Requires notparallel netcdf +Requires netcdf INPUT="-i cpptraj.in" -cat > cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in <