diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index fa81932b2a..9d59717d98 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -18088,6 +18088,18 @@ ucell Write older (v21) format trajectory that stores unit cell params instead veltraj Write velocity trajectory instead of coordinates. \end_layout +\begin_layout Description +dt Set trajectory time step in ps. +\end_layout + +\begin_layout Description +nstep # steps between frames. +\end_layout + +\begin_layout Description +step0 Initial step. +\end_layout + \end_deeper \begin_layout Standard Note that by default CPPTRAJ will try to write the symmetric shape matrix @@ -38950,7 +38962,7 @@ unwrap \end_layout \begin_layout LyX-Code -unwrap [center] [{bymol | byres | byatom}] +unwrap [center] [{byatom | byres | bymol}] \end_layout \begin_layout LyX-Code @@ -38963,11 +38975,12 @@ unwrap [center] [{bymol | byres | byatom}] \begin_deeper \begin_layout Description -[center] Unwrap by center of mass; otherwise unwrap by first atom position. +[center] Unwrap by center of mass; otherwise unwrap by first atom position + (byres and bymol). \end_layout \begin_layout Description -bymol Unwrap by molecule (default). +byatom Unwrap by atom (default). \end_layout \begin_layout Description @@ -38975,7 +38988,7 @@ byres Unwrap by residue. \end_layout \begin_layout Description -byatom Unwrap by atom. +bymol Unwrap by molecule. \end_layout \begin_layout Description @@ -39048,28 +39061,20 @@ image \series default , but this command can also be used to place molecules side by side, for example, two strands of a DNA duplex. - However, this command fails when the + However, this command may fail if the \shape italic mask \shape default -ed molecules travel more than half of the box size within a single frame. +ed entities travel more than half of the box size within a single frame. + Note that uwrapping by atom (the default behavior) is slower than unwrapping + by residue or molecule, but is usually the safer method, especially if + it is unknown if the original imaging was done by atom, by residue, or + by molecule. \end_layout \begin_layout Standard -If the optional argument -\begin_inset Quotes eld -\end_inset - - -\family typewriter -reference -\family default - -\begin_inset Quotes erd -\end_inset - - is specified, then the first frame is unwrapped according to the reference - structure. +If the optional reference arguments are specified, then the first frame + is unwrapped according to the reference structure. Otherwise, the first frame is not modified. \end_layout diff --git a/src/Action_Diffusion.cpp b/src/Action_Diffusion.cpp index 4363f04112..e21605b7bf 100644 --- a/src/Action_Diffusion.cpp +++ b/src/Action_Diffusion.cpp @@ -1,8 +1,7 @@ -#include // sqrt +#include // sqrt, round #include "Action_Diffusion.h" #include "CpptrajStdio.h" #include "StringRoutines.h" // validDouble -#include "Unwrap.h" #include "DataSet_1D.h" // LinearRegression #ifdef TIMER # include "Timer.h" @@ -11,8 +10,6 @@ # include #endif -using namespace Cpptraj; - // CONSTRUCTOR Action_Diffusion::Action_Diffusion() : avg_x_(0), avg_y_(0), avg_z_(0), avg_r_(0), avg_a_(0), @@ -251,8 +248,9 @@ Action::RetType Action_Diffusion::Setup(ActionSetup& setup) { } else mprintf("\tImaging disabled.\n"); - // Reserve space for the previous coordinates array - previous_.reserve( mask_.Nselected() * 3 ); + // If imaging, reserve space for the previous fractional coordinates array + if (imageOpt_.ImagingEnabled()) + previousFrac_.reserve( mask_.Nselected() ); // If initial frame already set and current # atoms > # atoms in initial // frame this will probably cause an error. @@ -310,7 +308,11 @@ Action::RetType Action_Diffusion::Setup(ActionSetup& setup) { // Action_Diffusion::DoAction() Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) { - // Load initial frame if necessary + if (imageOpt_.ImagingEnabled()) { + imageOpt_.SetImageType( frm.Frm().BoxCrd().Is_X_Aligned_Ortho() ); + } + + // ----- Load initial frame if necessary ------- if (initial_.empty()) { initial_ = frm.Frm(); # ifdef MPI @@ -322,21 +324,15 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) { initial_.RecvFrame( 0, trajComm_ ); } # endif - for (AtomMask::const_iterator atom = mask_.begin(); atom != mask_.end(); ++atom) - { - const double* XYZ = initial_.XYZ(*atom); - previous_.push_back( XYZ[0] ); - previous_.push_back( XYZ[1] ); - previous_.push_back( XYZ[2] ); + if (imageOpt_.ImagingEnabled()) { + // Imaging. Store initial fractional coordinates. + for (AtomMask::const_iterator atom = mask_.begin(); atom != mask_.end(); ++atom) + { + previousFrac_.push_back( frm.Frm().BoxCrd().FracCell() * Vec3( initial_.XYZ(*atom) ) ); + } } - } - // Diffusion calculation - Vec3 boxLengths(0.0), limxyz(0.0); - if (imageOpt_.ImagingEnabled()) { - imageOpt_.SetImageType( frm.Frm().BoxCrd().Is_X_Aligned_Ortho() ); - boxLengths = frm.Frm().BoxCrd().Lengths(); - limxyz = boxLengths / 2.0; - } + } // END initial frame load + // ----- Diffusion calculation ----------------- // For averaging over selected atoms double average2 = 0.0; double avgx = 0.0; @@ -351,40 +347,30 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) { for (imask = 0; imask < mask_.Nselected(); imask++) { int at = mask_[imask]; - int idx = imask * 3; // Index into previous_ - // Get current and initial coords for this atom. + // Get current coords for this atom. const double* XYZ = frm.Frm().XYZ(at); - double fixedXYZ[3]; - fixedXYZ[0] = XYZ[0]; - fixedXYZ[1] = XYZ[1]; - fixedXYZ[2] = XYZ[2]; + // Get initial coords for this atom const double* iXYZ = initial_.XYZ(at); // Calculate distance from initial position. double delx, dely, delz; - if ( imageOpt_.ImagingType() == ImageOption::ORTHO ) { - // Orthorhombic imaging - // Calculate vector needed to correct XYZ for imaging. - Vec3 transVec = Unwrap::UnwrapVec_Ortho(Vec3(XYZ), Vec3((&previous_[0])+idx), boxLengths, limxyz); - fixedXYZ[0] += transVec[0]; - fixedXYZ[1] += transVec[1]; - fixedXYZ[2] += transVec[2]; - // Calculate the distance between this "fixed" coordinate - // and the reference (initial) frame. - delx = fixedXYZ[0] - iXYZ[0]; - dely = fixedXYZ[1] - iXYZ[1]; - delz = fixedXYZ[2] - iXYZ[2]; - } else if ( imageOpt_.ImagingType() == ImageOption::NONORTHO ) { - // Non-orthorhombic imaging - // Calculate vector needed to correct XYZ for imaging. - Vec3 transVec = Unwrap::UnwrapVec_Nonortho(Vec3(XYZ), Vec3((&previous_[0])+idx), frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell(), limxyz); - fixedXYZ[0] += transVec[0]; - fixedXYZ[1] += transVec[1]; - fixedXYZ[2] += transVec[2]; - // Calculate the distance between this "fixed" coordinate - // and the reference (initial) frame. - delx = fixedXYZ[0] - iXYZ[0]; - dely = fixedXYZ[1] - iXYZ[1]; - delz = fixedXYZ[2] - iXYZ[2]; + if ( imageOpt_.ImagingEnabled() ) { + // Convert current Cartesian coords to fractional coords + Vec3 xyz_frac = frm.Frm().BoxCrd().FracCell() * Vec3(XYZ); + // Correct current frac coords + Vec3 ixyz = xyz_frac - previousFrac_[imask]; + ixyz[0] = ixyz[0] - round(ixyz[0]); + ixyz[1] = ixyz[1] - round(ixyz[1]); + ixyz[2] = ixyz[2] - round(ixyz[2]); + xyz_frac = previousFrac_[imask] + ixyz; + // Back to Cartesian + Vec3 xyz_cart1 = frm.Frm().BoxCrd().UnitCell().TransposeMult( xyz_frac ); + // Update reference frac coords + previousFrac_[imask] = xyz_frac; + // Calculate the distance between the fixed coordinates + // and reference (initial) frame coordinates. + delx = xyz_cart1[0] - iXYZ[0]; + dely = xyz_cart1[1] - iXYZ[1]; + delz = xyz_cart1[2] - iXYZ[2]; } else { // No imaging. Calculate distance from current position to initial position. delx = XYZ[0] - iXYZ[0]; @@ -415,10 +401,6 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) { fval = (float)dist2; atom_a_[at]->Add(frameNum, &fval); } - // Update the previous coordinate set to match the current coordinates - previous_[idx ] = fixedXYZ[0]; - previous_[idx+1] = fixedXYZ[1]; - previous_[idx+2] = fixedXYZ[2]; } // END loop over selected atoms # ifdef _OPENMP } // END pragma omp parallel diff --git a/src/Action_Diffusion.h b/src/Action_Diffusion.h index 176591c465..e0f50136c6 100644 --- a/src/Action_Diffusion.h +++ b/src/Action_Diffusion.h @@ -14,7 +14,7 @@ class Action_Diffusion : public Action { void Print(); typedef DataSetList::DataListType Dlist; - typedef std::vector Darray; + typedef std::vector Varray; inline void LoadInitial(Frame const&); void CalcDiffForSet(unsigned int&, Dlist const&, int, std::string const&) const; @@ -22,7 +22,7 @@ class Action_Diffusion : public Action { ImageOption imageOpt_; ///< Used to determine if imaging should be used. Frame initial_; ///< Initial frame (all atoms) - Darray previous_; ///< Previous coordinates (selected atoms) + Varray previousFrac_; ///< Previous fractional coordinates for imaging (selected atoms) DataSet* avg_x_; ///< Hold average diffusion in X direction each frame DataSet* avg_y_; ///< Hold average diffusion in Y direction each frame DataSet* avg_z_; ///< Hold average diffusion in Z direction each frame diff --git a/src/Action_Unwrap.cpp b/src/Action_Unwrap.cpp index 8ec969bb89..1e52a62732 100644 --- a/src/Action_Unwrap.cpp +++ b/src/Action_Unwrap.cpp @@ -11,7 +11,8 @@ Action_Unwrap::Action_Unwrap() : imageList_(0), imageMode_(Image::BYATOM), RefParm_(0), - center_(false) + center_(false), + refNeedsCalc_(false) { } /** DESTRUCTOR */ @@ -20,10 +21,11 @@ Action_Unwrap::~Action_Unwrap() { } void Action_Unwrap::Help() const { - mprintf("\t[center] [{bymol | byres | byatom}]\n" + mprintf("\t[center] [{byatom | byres | bymol}]\n" "\t[ %s ] []\n", DataSetList::RefArgs); mprintf(" Reverse of 'image'; unwrap coordinates in according\n" - " to a reference structure.\n"); + " to the first frame, or optionally a reference structure. Can\n" + " unwrap by atom (default), by residue, or by molecule.\n"); } // Action_Unwrap::Init() @@ -55,7 +57,9 @@ Action::RetType Action_Unwrap::Init(ArgList& actionArgs, ActionInit& init, int d RefFrame_ = REF.Coord(); // Get reference parm for frame RefParm_ = REF.ParmPtr(); - } + refNeedsCalc_ = true; + } else + refNeedsCalc_ = false; // Get mask string maskExpression_ = actionArgs.GetMaskNext(); @@ -132,18 +136,15 @@ Action::RetType Action_Unwrap::Setup(ActionSetup& setup) { // Action_Unwrap::DoAction() Action::RetType Action_Unwrap::DoAction(int frameNum, ActionFrame& frm) { - if (RefFrame_.empty()) { - // Set reference structure if not already set - RefFrame_ = frm.Frm(); - return Action::OK; + if (refNeedsCalc_) { + // Calculate initial fractional coords from reference frame. + Image::UnwrapFrac(fracCoords_, RefFrame_, *imageList_, + frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell()); + refNeedsCalc_ = false; } - Vec3 limxyz = frm.Frm().BoxCrd().Lengths() / 2.0; - if (frm.Frm().BoxCrd().Is_X_Aligned_Ortho()) - Image::UnwrapOrtho( frm.ModifyFrm(), RefFrame_, *imageList_, allEntities_, limxyz ); - else { - Image::UnwrapNonortho( frm.ModifyFrm(), RefFrame_, *imageList_, allEntities_, frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell(), limxyz ); - } + Image::UnwrapFrac(fracCoords_, frm.ModifyFrm(), *imageList_, + frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell()); return Action::MODIFY_COORDS; } diff --git a/src/Action_Unwrap.h b/src/Action_Unwrap.h index 59afec8b0e..01cccfd5f0 100644 --- a/src/Action_Unwrap.h +++ b/src/Action_Unwrap.h @@ -23,9 +23,12 @@ class Action_Unwrap : public Action { Frame RefFrame_; ///< Reference frame, updated each DoAction Topology* RefParm_; ///< Reference topology bool center_; ///< If true, determine distances to centers + bool refNeedsCalc_; ///< If true, need to calc frac. coords of ref Unit allEntities_; ///< Hold atoms to copy from target to reference # ifdef MPI Parallel::Comm trajComm_; # endif + typedef std::vector Varray; + Varray fracCoords_; ///< Hold fractional coords for last frame }; #endif diff --git a/src/Frame.h b/src/Frame.h index 6e3a8d8836..5f543fdb23 100644 --- a/src/Frame.h +++ b/src/Frame.h @@ -98,6 +98,8 @@ class Frame { int RepIdx() const { return repidx_; } /// \return overall coordinate index int CrdIdx() const { return crdidx_; } + /// Set XYZ for atom + inline void SetXYZ(int, Vec3 const&); /// Set box from another box void SetBox( Box const& b ) { box_ = b; } /// Modify box in place @@ -305,6 +307,12 @@ class Frame { bool setupFrame(unsigned int, CoordinateInfo const&); }; // ---------- INLINE FUNCTION DEFINITIONS -------------------------------------- +void Frame::SetXYZ(int atnum, Vec3 const& xyz) { + int idx = atnum*3; + X_[idx ] = xyz[0]; + X_[idx+1] = xyz[1]; + X_[idx+2] = xyz[2]; +} bool Frame::CheckCoordsInvalid() const { if (natom_ > 1) { diff --git a/src/ImageRoutines.cpp b/src/ImageRoutines.cpp index 18ea5966a1..3c7f42af51 100644 --- a/src/ImageRoutines.cpp +++ b/src/ImageRoutines.cpp @@ -1,8 +1,7 @@ -#include // floor +#include // floor, round #include "ImageRoutines.h" #include "DistRoutines.h" #include "CpptrajStdio.h" -#include "Unwrap.h" #include "Image_List.h" #include "Image_List_Pair_CoM.h" #include "Image_List_Pair_Geom.h" @@ -12,8 +11,6 @@ #include "Image_List_Unit_First.h" #include "Image_List_Mask.h" -using namespace Cpptraj; - /** \return Empty list for imaging. */ Image::List* Image::CreateImageList(Mode modeIn, bool useMass, bool center) { Image::List* listOut = 0; @@ -247,60 +244,66 @@ Vec3 Image::Ortho(Vec3 const& Coord, Vec3 const& bp, Vec3 const& bm, Box const& } // ----------------------------------------------------------------------------- -// Image::UnwrapNonortho() -void Image::UnwrapNonortho(Frame& tgtIn, Frame& refIn, List const& AtomPairs, - Unit const& allEntities, - Matrix_3x3 const& ucell, Matrix_3x3 const& frac, - Vec3 const& limxyz) -{ - int idx; - int maxidx = (int)AtomPairs.nEntities(); - // Loop over atom pairs -# ifdef _OPENMP -# pragma omp parallel private(idx) - { -# pragma omp for -# endif - for (idx = 0; idx < maxidx; idx++) - { - Vec3 vtgt = AtomPairs.GetCoord(idx, tgtIn); - Vec3 vref = AtomPairs.GetCoord(idx, refIn); - Vec3 boxTrans = Unwrap::UnwrapVec_Nonortho(vtgt, vref, ucell, frac, limxyz); - AtomPairs.DoTranslation( tgtIn, idx, boxTrans ); - } // END loop over atom pairs -# ifdef _OPENMP - } -# endif - // Save new ref positions - refIn.CopyFrom(tgtIn, allEntities); -} - -// Image::UnwrapOrtho() -void Image::UnwrapOrtho(Frame& tgtIn, Frame& refIn, List const& AtomPairs, - Unit const& allEntities, Vec3 const& limxyz) +void Image::UnwrapFrac(std::vector& previousFrac, + Frame& currentFrame, + List const& AtomPairs, + Matrix_3x3 const& ucell, Matrix_3x3 const& frac) { - Vec3 boxVec = tgtIn.BoxCrd().Lengths(); int idx; int maxidx = (int)AtomPairs.nEntities(); - // Loop over atom pairs -# ifdef _OPENMP -# pragma omp parallel private(idx) - { -# pragma omp for -# endif - for (idx = 0; idx < maxidx; idx++) - { - Vec3 vtgt = AtomPairs.GetCoord(idx, tgtIn); - Vec3 vref = AtomPairs.GetCoord(idx, refIn); - Vec3 boxTrans = Unwrap::UnwrapVec_Ortho(vtgt, vref, boxVec, limxyz); - // Translate atoms from first to last - AtomPairs.DoTranslation( tgtIn, idx, boxTrans ); - } // END loop over atom pairs -# ifdef _OPENMP + if (previousFrac.empty()) { + // Set initial frac coords + //mprintf("DEBUG: Initial set.\n"); + previousFrac.resize( maxidx ); +# ifdef _OPENMP +# pragma omp parallel private(idx) + { +# pragma omp for +# endif + for (idx = 0; idx < maxidx; idx++) + { + // Convert to fractional coords + //Vec3 xyz_cart( currentFrame.XYZ( idx ) ); + Vec3 xyz_cart = AtomPairs.GetCoord(idx, currentFrame); + //Vec3 xyz_frac = frac * xyz_cart; + previousFrac[idx] = ( frac * xyz_cart ); + //previousFrac.push_back( xyz_frac[1] ); + //previousFrac.push_back( xyz_frac[2] ); + } +# ifdef _OPENMP + } +# endif + } else { + //mprintf("DEBUG: Subsequent set.\n"); + // Update currentframe +# ifdef _OPENMP +# pragma omp parallel private(idx) + { +# pragma omp for +# endif + for (idx = 0; idx < maxidx; idx++) + { + // Convert to fractional coords + //Vec3 xyz_cart( currentFrame.XYZ( idx ) ); + Vec3 xyz_cart0 = AtomPairs.GetCoord(idx, currentFrame); + Vec3 xyz_frac = frac * xyz_cart0; + // Correct frac coords + Vec3 ixyz = xyz_frac - previousFrac[idx]; + ixyz[0] = ixyz[0] - round(ixyz[0]); + ixyz[1] = ixyz[1] - round(ixyz[1]); + ixyz[2] = ixyz[2] - round(ixyz[2]); + xyz_frac = previousFrac[idx] + ixyz; + // Back to Cartesian + Vec3 xyz_cart1 = ucell.TransposeMult( xyz_frac ); + //currentFrame.SetXYZ(idx, xyz_cart); + AtomPairs.DoTranslation( currentFrame, idx, xyz_cart1 - xyz_cart0 ); + // Update reference frac coords + previousFrac[idx] = xyz_frac; + } +# ifdef _OPENMP + } +# endif } -# endif - // Save new ref positions - refIn.CopyFrom(tgtIn, allEntities); } // ----------------------------------------------------------------------------- diff --git a/src/ImageRoutines.h b/src/ImageRoutines.h index bce75e4c9a..4feda053e3 100644 --- a/src/ImageRoutines.h +++ b/src/ImageRoutines.h @@ -34,11 +34,9 @@ namespace Image { /// Perform orthogonal imaging on given coordinates using given box boundaries Vec3 Ortho(Vec3 const&, Vec3 const&, Vec3 const&, Box const&); - /// Perform unwrap of non-orthogonal cell using given reference. - void UnwrapNonortho( Frame&, Frame&, List const&, Unit const&, - Matrix_3x3 const&, Matrix_3x3 const&, Vec3 const& ); - /// Perform unwrap of orthogonal cell using given reference. - void UnwrapOrtho( Frame&, Frame&, List const&, Unit const&, Vec3 const& ); + /// Perform unwrapping in fractional space + void UnwrapFrac(std::vector&, Frame&, List const&, + Matrix_3x3 const&, Matrix_3x3 const&); /// Wrap selected atom coords from given frame into primary cell, store in result. void WrapToCell0(std::vector&, Frame const&, AtomMask const&, diff --git a/src/Traj_CharmmDcd.cpp b/src/Traj_CharmmDcd.cpp index 10466f41ec..52f034985f 100644 --- a/src/Traj_CharmmDcd.cpp +++ b/src/Traj_CharmmDcd.cpp @@ -28,7 +28,10 @@ Traj_CharmmDcd::Traj_CharmmDcd() : freeat_(0), xcoord_(0), ycoord_(0), - zcoord_(0) + zcoord_(0), + timeStep_(0), + stepsBetweenFrames_(0), + initialStep_(0) {} // DESTRUCTOR @@ -296,8 +299,8 @@ int Traj_CharmmDcd::setupTrajin(FileName const& fname, Topology* trajParm) box.SetupFromXyzAbg( boxtmp ); } } - // Set traj info: No velocity, temperature, or time. - SetCoordInfo( CoordinateInfo( box, false, false, false ) ); + // Set traj info: No velocity or temperature, but has time. + SetCoordInfo( CoordinateInfo( box, false, false, true ) ); // If there are fixed atoms read the first frame now // TODO: Deal with fixed atoms closeTraj(); @@ -365,8 +368,14 @@ int Traj_CharmmDcd::readDcdHeader() { } else boxBytes_ = 0; // Timestep - convert from AKMA to ps - float timestep = buffer.f[9] / Constants::AMBERTIME_TO_PS; - if (debug_>0) mprintf("\tTimestep is %f\n",timestep); + timeStep_ = (double)(buffer.f[9]) / Constants::AMBERTIME_TO_PS; + stepsBetweenFrames_ = buffer.i[2]; + initialStep_= buffer.i[1]; + if (debug_>0) { + mprintf("\tTimestep is %f ps\n", timeStep_); + mprintf("\tSteps between frames is %i\n", stepsBetweenFrames_); + mprintf("\tStarting step is %i\n", initialStep_); + } // Read end size of first block, should also be 84 if (ReadBlock(84)<0) return 1; // ********** Step 2 - Read title block @@ -548,6 +557,9 @@ int Traj_CharmmDcd::readFrame(int set, Frame& frameIn) { else frameIn.ModifyBox().AssignFromXyzAbg( box ); } + // Calculate time + int step = initialStep_ + (set * stepsBetweenFrames_); + frameIn.SetTime( (double)step * timeStep_ ); return readXYZ(frameIn.xAddress()); } @@ -572,6 +584,9 @@ void Traj_CharmmDcd::WriteHelp() { "\tucell : Write older (v21) format trajectory that stores unit cell params\n" "\t instead of shape matrix.\n" "\tveltraj : Write velocity trajectory instead of coordinates.\n" + "\tdt : Set trajectory time step in ps.\n" + "\tnstep : # steps between frames.\n" + "\tstep0 : Initial step.\n" ); } @@ -593,6 +608,9 @@ int Traj_CharmmDcd::processWriteArgs(ArgList& argIn, DataSetList const& DSLin) { isVel_ = true; else isVel_ = false; + timeStep_ = argIn.getKeyDouble("dt", 0.001); + stepsBetweenFrames_ = argIn.getKeyInt("nstep", 1); + initialStep_ = argIn.getKeyInt("step0", 1); return 0; } @@ -689,16 +707,16 @@ int Traj_CharmmDcd::writeDcdHeader() { //buffer.i[0] = trajParm->parmFrames; buffer.i[0] = 0; // Starting timestep - buffer.i[1] = 1; + buffer.i[1] = initialStep_; // Number of steps between frames - buffer.i[2] = 1; + buffer.i[2] = stepsBetweenFrames_; // For velocity traj only if (isVel_) buffer.i[4] = 1; // Number of fixed atoms buffer.i[8] = 0; - // Timestep - buffer.f[9] = 0.001; + // Timestep - convert from ps to AKMA + buffer.f[9] = (float)(timeStep_ * Constants::AMBERTIME_TO_PS); // Default to SHAPE if (charmmCellType_ == UNKNOWN) charmmCellType_ = SHAPE; diff --git a/src/Traj_CharmmDcd.h b/src/Traj_CharmmDcd.h index 9bc46c587e..1ec8d7fada 100644 --- a/src/Traj_CharmmDcd.h +++ b/src/Traj_CharmmDcd.h @@ -32,6 +32,9 @@ class Traj_CharmmDcd : public TrajectoryIO { float* ycoord_; ///< Pointer to start of Y coords in master coord array float* zcoord_; ///< Pointer to start of Z coords in master coord array CpptrajFile file_; ///< Input/Output file + double timeStep_; ///< Time step for writing + int stepsBetweenFrames_; ///< # of steps between frames (write) + int initialStep_; ///< Initial step (write) union headerbyte { unsigned char c[80]; int i[20]; float f[20]; }; int ReadBlock(int); diff --git a/src/Unwrap.h b/src/Unwrap.h deleted file mode 100644 index 4e08bc46d5..0000000000 --- a/src/Unwrap.h +++ /dev/null @@ -1,73 +0,0 @@ -#ifndef INC_UNWRAP_H -#define INC_UNWRAP_H -#include "Matrix_3x3.h" -#include "Vec3.h" -namespace Cpptraj { -namespace Unwrap { - /// \return True if absolute value of any component of dxyz is greater than limxyz - template bool vGreaterThan(T const* dxyz, T const* limxyz) { - return ( fabs(dxyz[0]) > limxyz[0] || - fabs(dxyz[1]) > limxyz[1] || - fabs(dxyz[2]) > limxyz[2] ); - } - - /// \return Vector required for unwrapping in orthogonal box - template Vec3 UnwrapVec_Ortho(Vec3 const& vtgt, Vec3 const& vref, Vec3 const& boxVec, Vec3 const& limxyz) - { - T dxyz[3]; - dxyz[0] = (vtgt[0] - vref[0]); - dxyz[1] = (vtgt[1] - vref[1]); - dxyz[2] = (vtgt[2] - vref[2]); - if (vGreaterThan(dxyz, limxyz.Dptr())) - return Vec3( - -floor( dxyz[0] / boxVec[0] + 0.5 ) * boxVec[0], - -floor( dxyz[1] / boxVec[1] + 0.5 ) * boxVec[1], - -floor( dxyz[2] / boxVec[2] + 0.5 ) * boxVec[2] - ); - else - return Vec3(0.0); - } - - /// \return Vector required for unwrapping in non-orthogonal box - template Vec3 UnwrapVec_Nonortho(Vec3 const& vtgt, Vec3 const& vref, Matrix_3x3 const& ucell, Matrix_3x3 const& frac, Vec3 const& limxyz) - { - Vec3 boxTrans(0.0); - // Calculate original distance from the ref (previous) position. - Vec3 vd = vtgt - vref; // dx dy dz - if (!vGreaterThan(vd.Dptr(), limxyz.Dptr())) - return boxTrans; - T minDistanceSquare = vd.Magnitude2(); - // Reciprocal coordinates - vd = frac * vd ; // recip * dxyz - T cx = floor(vd[0]); - T cy = floor(vd[1]); - T cz = floor(vd[2]); - // Loop over all possible translations - for (int ix = -1; ix < 2; ++ix) { - for (int iy = -1; iy < 2; ++iy) { - for (int iz = -1; iz < 2; ++iz) { - // Calculate the translation. - Vec3 vcc = ucell.TransposeMult( Vec3( cx+(T)ix, - cy+(T)iy, - cz+(T)iz ) ); // ucell^T * ccxyz - // Calc. the potential new coordinate for tgt - Vec3 vnew = vtgt - vcc; - // Calc. the new distance from the ref (previous) position - Vec3 vr = vref - vnew; - T distanceSquare = vr.Magnitude2(); - // If the orig. distance is greater than the new distance, unwrap. - if ( minDistanceSquare > distanceSquare ) { - minDistanceSquare = distanceSquare; - boxTrans = vcc; - } - } - } - } - // Translate tgt atoms - boxTrans.Neg(); - return boxTrans; - } - -} // END namespace Unwrap -} // END namespace Cpptraj -#endif diff --git a/src/Version.h b/src/Version.h index a3a82267eb..5db18168fd 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 "V6.18.1" +#define CPPTRAJ_INTERNAL_VERSION "V6.19.0" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 471813c72a..b89b8356f8 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -23,7 +23,7 @@ Action_CreateReservoir.o : Action_CreateReservoir.cpp Action.h ActionState.h Act Action_DNAionTracker.o : Action_DNAionTracker.cpp Action.h ActionState.h Action_DNAionTracker.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h Action_DSSP.o : Action_DSSP.cpp Action.h ActionState.h Action_DSSP.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h Action_Density.o : Action_Density.cpp Action.h ActionState.h Action_Density.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OnlineVarT.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h -Action_Diffusion.o : Action_Diffusion.cpp Action.h ActionState.h Action_Diffusion.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Unwrap.h Vec3.h +Action_Diffusion.o : Action_Diffusion.cpp Action.h ActionState.h Action_Diffusion.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h Action_Dihedral.o : Action_Dihedral.cpp Action.h ActionState.h Action_Dihedral.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_double.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TorsionRoutines.h TypeNameHolder.h Unit.h Vec3.h Action_DihedralRMS.o : Action_DihedralRMS.cpp Action.h ActionState.h Action_DihedralRMS.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DihedralSearch.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceAction.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TorsionRoutines.h TypeNameHolder.h Unit.h Vec3.h Action_Dipole.o : Action_Dipole.cpp Action.h ActionState.h Action_Dipole.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h Grid.h GridAction.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h @@ -350,7 +350,7 @@ GridAction.o : GridAction.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h GridBin.o : GridBin.cpp Box.h CpptrajStdio.h GridBin.h Matrix_3x3.h Parallel.h Vec3.h HistBin.o : HistBin.cpp Constants.h CpptrajStdio.h Dimension.h HistBin.h Hungarian.o : Hungarian.cpp ArrayIterator.h Constants.h CpptrajStdio.h Hungarian.h Matrix.h -ImageRoutines.o : ImageRoutines.cpp Atom.h AtomMask.h Box.h CoordinateInfo.h CpptrajStdio.h DistRoutines.h Frame.h ImageOption.h ImageRoutines.h ImageTypes.h Image_List.h Image_List_Mask.h Image_List_Pair.h Image_List_Pair_CoM.h Image_List_Pair_First.h Image_List_Pair_Geom.h Image_List_Unit.h Image_List_Unit_CoM.h Image_List_Unit_First.h Image_List_Unit_Geom.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Unit.h Unwrap.h Vec3.h +ImageRoutines.o : ImageRoutines.cpp Atom.h AtomMask.h Box.h CoordinateInfo.h CpptrajStdio.h DistRoutines.h Frame.h ImageOption.h ImageRoutines.h ImageTypes.h Image_List.h Image_List_Mask.h Image_List_Pair.h Image_List_Pair_CoM.h Image_List_Pair_First.h Image_List_Pair_Geom.h Image_List_Unit.h Image_List_Unit_CoM.h Image_List_Unit_First.h Image_List_Unit_Geom.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Unit.h Vec3.h Image_List_Mask.o : Image_List_Mask.cpp Atom.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h Image_List.h Image_List_Mask.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h Image_List_Pair.o : Image_List_Pair.cpp Atom.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h Image_List.h Image_List_Pair.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h Image_List_Unit.o : Image_List_Unit.cpp Atom.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h Image_List.h Image_List_Unit.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h diff --git a/test/Test_Unwrap/RunTest.sh b/test/Test_Unwrap/RunTest.sh index 571a07f006..d6f8368b2f 100755 --- a/test/Test_Unwrap/RunTest.sh +++ b/test/Test_Unwrap/RunTest.sh @@ -10,7 +10,7 @@ INPUT="ptraj.in" TOP="../tz2.truncoct.parm7" cat > ptraj.in < ptraj.in <