From 8125772bb5ee3bd5fcccf8db675376911c7b8e90 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 9 May 2023 15:34:47 -0400 Subject: [PATCH 01/16] Do unwrapping entirely in fractional space to avoid problems unwrapping NPT trajectories --- src/Action_Unwrap.cpp | 4 ++- src/Action_Unwrap.h | 2 ++ src/Frame.h | 8 ++++++ src/ImageRoutines.cpp | 60 ++++++++++++++++++++++++++++++++++++++++++- src/ImageRoutines.h | 4 +++ 5 files changed, 76 insertions(+), 2 deletions(-) diff --git a/src/Action_Unwrap.cpp b/src/Action_Unwrap.cpp index 8ec969bb89..ff040de7b9 100644 --- a/src/Action_Unwrap.cpp +++ b/src/Action_Unwrap.cpp @@ -132,6 +132,8 @@ Action::RetType Action_Unwrap::Setup(ActionSetup& setup) { // Action_Unwrap::DoAction() Action::RetType Action_Unwrap::DoAction(int frameNum, ActionFrame& frm) { + Image::UnwrapFrac(fracCoords_, frm.ModifyFrm(), frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell()); +/* if (RefFrame_.empty()) { // Set reference structure if not already set RefFrame_ = frm.Frm(); @@ -144,6 +146,6 @@ Action::RetType Action_Unwrap::DoAction(int frameNum, ActionFrame& frm) { else { Image::UnwrapNonortho( frm.ModifyFrm(), RefFrame_, *imageList_, allEntities_, frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell(), limxyz ); } - +*/ return Action::MODIFY_COORDS; } diff --git a/src/Action_Unwrap.h b/src/Action_Unwrap.h index 59afec8b0e..09dac40c08 100644 --- a/src/Action_Unwrap.h +++ b/src/Action_Unwrap.h @@ -27,5 +27,7 @@ class Action_Unwrap : public Action { # 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..3c75b3c75b 100644 --- a/src/ImageRoutines.cpp +++ b/src/ImageRoutines.cpp @@ -1,4 +1,4 @@ -#include // floor +#include // floor, round #include "ImageRoutines.h" #include "DistRoutines.h" #include "CpptrajStdio.h" @@ -247,6 +247,64 @@ Vec3 Image::Ortho(Vec3 const& Coord, Vec3 const& bp, Vec3 const& bm, Box const& } // ----------------------------------------------------------------------------- +void Image::UnwrapFrac(std::vector& previousFrac, + Frame& currentFrame, + Matrix_3x3 const& ucell, Matrix_3x3 const& frac) +{ + int idx; + int natom = currentFrame.Natom(); + if (previousFrac.empty()) { + // Set initial frac coords + //mprintf("DEBUG: Initial set.\n"); + previousFrac.reserve( natom * 3 ); +# ifdef _OPENMP +# pragma omp parallel private(idx) + { +# pragma omp for +# endif + for (idx = 0; idx < natom; idx++) + { + // Convert to fractional coords + Vec3 xyz_cart( currentFrame.XYZ( idx ) ); + //Vec3 xyz_frac = frac * xyz_cart; + previousFrac.push_back( 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 < natom; idx++) + { + // Convert to fractional coords + Vec3 xyz_cart( currentFrame.XYZ( idx ) ); + Vec3 xyz_frac = frac * xyz_cart; + // 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 + xyz_cart = ucell.TransposeMult( xyz_frac ); + currentFrame.SetXYZ(idx, xyz_cart); + // Update reference frac coords + previousFrac[idx] = xyz_frac; + } +# ifdef _OPENMP + } +# endif + } +} + // Image::UnwrapNonortho() void Image::UnwrapNonortho(Frame& tgtIn, Frame& refIn, List const& AtomPairs, Unit const& allEntities, diff --git a/src/ImageRoutines.h b/src/ImageRoutines.h index bce75e4c9a..657a88da24 100644 --- a/src/ImageRoutines.h +++ b/src/ImageRoutines.h @@ -34,6 +34,10 @@ namespace Image { /// Perform orthogonal imaging on given coordinates using given box boundaries Vec3 Ortho(Vec3 const&, Vec3 const&, Vec3 const&, Box const&); + /// Perform unwrapping in fractional space + void UnwrapFrac(std::vector&, Frame&, + Matrix_3x3 const&, Matrix_3x3 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& ); From 592b96b779abb3790863aea4708f40738bd72ca9 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 9 May 2023 15:47:38 -0400 Subject: [PATCH 02/16] Add dt keyword for dcd trajout to set time step. Ensure time step is written in AKMA units. --- src/Traj_CharmmDcd.cpp | 13 ++++++++----- src/Traj_CharmmDcd.h | 1 + 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/Traj_CharmmDcd.cpp b/src/Traj_CharmmDcd.cpp index 10466f41ec..994334bc8d 100644 --- a/src/Traj_CharmmDcd.cpp +++ b/src/Traj_CharmmDcd.cpp @@ -28,7 +28,8 @@ Traj_CharmmDcd::Traj_CharmmDcd() : freeat_(0), xcoord_(0), ycoord_(0), - zcoord_(0) + zcoord_(0), + timeStep_(0) {} // DESTRUCTOR @@ -365,8 +366,8 @@ 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; + if (debug_>0) mprintf("\tTimestep is %f ps\n", timeStep_); // Read end size of first block, should also be 84 if (ReadBlock(84)<0) return 1; // ********** Step 2 - Read title block @@ -572,6 +573,7 @@ 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" ); } @@ -593,6 +595,7 @@ int Traj_CharmmDcd::processWriteArgs(ArgList& argIn, DataSetList const& DSLin) { isVel_ = true; else isVel_ = false; + timeStep_ = argIn.getKeyDouble("dt", 0.001); return 0; } @@ -697,8 +700,8 @@ int Traj_CharmmDcd::writeDcdHeader() { 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..f10d374241 100644 --- a/src/Traj_CharmmDcd.h +++ b/src/Traj_CharmmDcd.h @@ -32,6 +32,7 @@ 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 union headerbyte { unsigned char c[80]; int i[20]; float f[20]; }; int ReadBlock(int); From cd4bd8f7a97e661542b668f06452ccd2c1b4959a Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 10 May 2023 09:23:45 -0400 Subject: [PATCH 03/16] Use the image list so bymol and byres are honored --- src/Action_Unwrap.cpp | 2 +- src/ImageRoutines.cpp | 22 +++++++++++++--------- src/ImageRoutines.h | 2 +- 3 files changed, 15 insertions(+), 11 deletions(-) diff --git a/src/Action_Unwrap.cpp b/src/Action_Unwrap.cpp index ff040de7b9..eafc69a902 100644 --- a/src/Action_Unwrap.cpp +++ b/src/Action_Unwrap.cpp @@ -132,7 +132,7 @@ Action::RetType Action_Unwrap::Setup(ActionSetup& setup) { // Action_Unwrap::DoAction() Action::RetType Action_Unwrap::DoAction(int frameNum, ActionFrame& frm) { - Image::UnwrapFrac(fracCoords_, frm.ModifyFrm(), frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell()); + Image::UnwrapFrac(fracCoords_, frm.ModifyFrm(), *imageList_, frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell()); /* if (RefFrame_.empty()) { // Set reference structure if not already set diff --git a/src/ImageRoutines.cpp b/src/ImageRoutines.cpp index 3c75b3c75b..00596110b1 100644 --- a/src/ImageRoutines.cpp +++ b/src/ImageRoutines.cpp @@ -249,23 +249,25 @@ Vec3 Image::Ortho(Vec3 const& Coord, Vec3 const& bp, Vec3 const& bm, Box const& // ----------------------------------------------------------------------------- void Image::UnwrapFrac(std::vector& previousFrac, Frame& currentFrame, + List const& AtomPairs, Matrix_3x3 const& ucell, Matrix_3x3 const& frac) { int idx; - int natom = currentFrame.Natom(); + int maxidx = (int)AtomPairs.nEntities(); if (previousFrac.empty()) { // Set initial frac coords //mprintf("DEBUG: Initial set.\n"); - previousFrac.reserve( natom * 3 ); + previousFrac.reserve( maxidx ); # ifdef _OPENMP # pragma omp parallel private(idx) { # pragma omp for # endif - for (idx = 0; idx < natom; idx++) + for (idx = 0; idx < maxidx; idx++) { // Convert to fractional coords - Vec3 xyz_cart( currentFrame.XYZ( idx ) ); + //Vec3 xyz_cart( currentFrame.XYZ( idx ) ); + Vec3 xyz_cart = AtomPairs.GetCoord(idx, currentFrame); //Vec3 xyz_frac = frac * xyz_cart; previousFrac.push_back( frac * xyz_cart ); //previousFrac.push_back( xyz_frac[1] ); @@ -282,11 +284,12 @@ void Image::UnwrapFrac(std::vector& previousFrac, { # pragma omp for # endif - for (idx = 0; idx < natom; idx++) + for (idx = 0; idx < maxidx; idx++) { // Convert to fractional coords - Vec3 xyz_cart( currentFrame.XYZ( idx ) ); - Vec3 xyz_frac = frac * xyz_cart; + //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]); @@ -294,8 +297,9 @@ void Image::UnwrapFrac(std::vector& previousFrac, ixyz[2] = ixyz[2] - round(ixyz[2]); xyz_frac = previousFrac[idx] + ixyz; // Back to Cartesian - xyz_cart = ucell.TransposeMult( xyz_frac ); - currentFrame.SetXYZ(idx, xyz_cart); + 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; } diff --git a/src/ImageRoutines.h b/src/ImageRoutines.h index 657a88da24..ce7362463e 100644 --- a/src/ImageRoutines.h +++ b/src/ImageRoutines.h @@ -35,7 +35,7 @@ namespace Image { Vec3 Ortho(Vec3 const&, Vec3 const&, Vec3 const&, Box const&); /// Perform unwrapping in fractional space - void UnwrapFrac(std::vector&, Frame&, + void UnwrapFrac(std::vector&, Frame&, List const&, Matrix_3x3 const&, Matrix_3x3 const&); /// Perform unwrap of non-orthogonal cell using given reference. From 3e349384ca10108495bd6ed9f972ed4b23085056 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 10 May 2023 10:06:16 -0400 Subject: [PATCH 04/16] Cannot use push_back in an OpenMP loop --- src/ImageRoutines.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ImageRoutines.cpp b/src/ImageRoutines.cpp index 00596110b1..a9d9a1d397 100644 --- a/src/ImageRoutines.cpp +++ b/src/ImageRoutines.cpp @@ -257,7 +257,7 @@ void Image::UnwrapFrac(std::vector& previousFrac, if (previousFrac.empty()) { // Set initial frac coords //mprintf("DEBUG: Initial set.\n"); - previousFrac.reserve( maxidx ); + previousFrac.resize( maxidx ); # ifdef _OPENMP # pragma omp parallel private(idx) { @@ -269,7 +269,7 @@ void Image::UnwrapFrac(std::vector& previousFrac, //Vec3 xyz_cart( currentFrame.XYZ( idx ) ); Vec3 xyz_cart = AtomPairs.GetCoord(idx, currentFrame); //Vec3 xyz_frac = frac * xyz_cart; - previousFrac.push_back( frac * xyz_cart ); + previousFrac[idx] = ( frac * xyz_cart ); //previousFrac.push_back( xyz_frac[1] ); //previousFrac.push_back( xyz_frac[2] ); } From 8d48c8405315713e87f37dc9f64df87a74da5005 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 10 May 2023 10:25:20 -0400 Subject: [PATCH 05/16] Fix imaging in diffusion calc for NPT simulations --- src/Action_Diffusion.cpp | 89 +++++++++++++++++----------------------- src/Action_Diffusion.h | 4 +- 2 files changed, 39 insertions(+), 54 deletions(-) diff --git a/src/Action_Diffusion.cpp b/src/Action_Diffusion.cpp index 4363f04112..a9c3600665 100644 --- a/src/Action_Diffusion.cpp +++ b/src/Action_Diffusion.cpp @@ -1,4 +1,4 @@ -#include // sqrt +#include // sqrt, round #include "Action_Diffusion.h" #include "CpptrajStdio.h" #include "StringRoutines.h" // validDouble @@ -251,8 +251,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 +311,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 +327,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 +350,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 +404,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 From 042a2eef3b4623a678ef2625684b0b4d07c46eed Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 10 May 2023 10:51:32 -0400 Subject: [PATCH 06/16] If using a reference structure, ensure frac coords get calcd for it --- src/Action_Unwrap.cpp | 28 +++++++++++++--------------- src/Action_Unwrap.h | 1 + 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/src/Action_Unwrap.cpp b/src/Action_Unwrap.cpp index eafc69a902..a7fe5e8c18 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 */ @@ -55,7 +56,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,20 +135,15 @@ Action::RetType Action_Unwrap::Setup(ActionSetup& setup) { // Action_Unwrap::DoAction() Action::RetType Action_Unwrap::DoAction(int frameNum, ActionFrame& frm) { - Image::UnwrapFrac(fracCoords_, frm.ModifyFrm(), *imageList_, frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell()); -/* - 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 09dac40c08..01cccfd5f0 100644 --- a/src/Action_Unwrap.h +++ b/src/Action_Unwrap.h @@ -23,6 +23,7 @@ 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_; From 1d8073b7d1533c0432f4934f8d7d58793e6b8bef Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 10 May 2023 10:59:46 -0400 Subject: [PATCH 07/16] Explicitly make tests byatom --- test/Test_Unwrap/RunTest.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 < Date: Wed, 10 May 2023 11:06:05 -0400 Subject: [PATCH 08/16] Update help --- src/Action_Unwrap.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Action_Unwrap.cpp b/src/Action_Unwrap.cpp index a7fe5e8c18..1e52a62732 100644 --- a/src/Action_Unwrap.cpp +++ b/src/Action_Unwrap.cpp @@ -21,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() From 193ebb6b5e4f4c6fbbdf6d44530401e36d9b9dcd Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 10 May 2023 11:10:52 -0400 Subject: [PATCH 09/16] Remove the old unwrap routines --- src/Action_Diffusion.cpp | 3 -- src/ImageRoutines.cpp | 59 ---------------------------------------- src/cpptrajdepend | 4 +-- 3 files changed, 2 insertions(+), 64 deletions(-) diff --git a/src/Action_Diffusion.cpp b/src/Action_Diffusion.cpp index a9c3600665..e21605b7bf 100644 --- a/src/Action_Diffusion.cpp +++ b/src/Action_Diffusion.cpp @@ -2,7 +2,6 @@ #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), diff --git a/src/ImageRoutines.cpp b/src/ImageRoutines.cpp index a9d9a1d397..3c7f42af51 100644 --- a/src/ImageRoutines.cpp +++ b/src/ImageRoutines.cpp @@ -2,7 +2,6 @@ #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; @@ -309,62 +306,6 @@ void Image::UnwrapFrac(std::vector& previousFrac, } } -// 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) -{ - 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 - } -# endif - // Save new ref positions - refIn.CopyFrom(tgtIn, allEntities); -} - // ----------------------------------------------------------------------------- void Image::WrapToCell0(std::vector& CoordsIn, Frame const& frmIn, AtomMask const& maskIn, 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 From e79f2e37146ebc0df390a458e8363a150129cfa4 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 10 May 2023 11:11:34 -0400 Subject: [PATCH 10/16] Remove old Unwrap templates --- src/Unwrap.h | 73 ---------------------------------------------------- 1 file changed, 73 deletions(-) delete mode 100644 src/Unwrap.h 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 From de79febc6e44c0c1fa53457e82e938cd5df6fa8b Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 10 May 2023 11:14:55 -0400 Subject: [PATCH 11/16] Fix up unwrap manual entry; note that default is actually unwrap by atom --- doc/cpptraj.lyx | 33 +++++++++++++-------------------- 1 file changed, 13 insertions(+), 20 deletions(-) diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index fa81932b2a..398d8233ad 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -38950,7 +38950,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 +38963,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 +38976,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 +39049,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 From 4fd05d50cbc3cace091765b1297a2843e0413ba8 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 10 May 2023 11:16:25 -0400 Subject: [PATCH 12/16] Add dt option to dcd trajout --- doc/cpptraj.lyx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 398d8233ad..015ea76b3c 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -18088,6 +18088,10 @@ 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 + \end_deeper \begin_layout Standard Note that by default CPPTRAJ will try to write the symmetric shape matrix From 95a32f7dd25f20fd9037b73ba47947132f44e346 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 10 May 2023 11:40:53 -0400 Subject: [PATCH 13/16] 6.19.0. Minor version bump for fixed unwrapping and fixed diffusion imaging for NPT --- src/Version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 90b271e5e70f9b1c35d31f363f91391813b85d59 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 10 May 2023 12:19:08 -0400 Subject: [PATCH 14/16] Properly read time from charmm dcd. Add t0 and nstep options to charmm trajout --- src/Traj_CharmmDcd.cpp | 27 +++++++++++++++++++++------ src/Traj_CharmmDcd.h | 2 ++ 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/src/Traj_CharmmDcd.cpp b/src/Traj_CharmmDcd.cpp index 994334bc8d..b79eedf761 100644 --- a/src/Traj_CharmmDcd.cpp +++ b/src/Traj_CharmmDcd.cpp @@ -29,7 +29,9 @@ Traj_CharmmDcd::Traj_CharmmDcd() : xcoord_(0), ycoord_(0), zcoord_(0), - timeStep_(0) + timeStep_(0), + stepsBetweenFrames_(0), + initialStep_(0) {} // DESTRUCTOR @@ -297,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(); @@ -367,7 +369,13 @@ int Traj_CharmmDcd::readDcdHeader() { boxBytes_ = 0; // Timestep - convert from AKMA to ps timeStep_ = (double)(buffer.f[9]) / Constants::AMBERTIME_TO_PS; - if (debug_>0) mprintf("\tTimestep is %f ps\n", timeStep_); + 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 @@ -549,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()); } @@ -574,6 +585,8 @@ void Traj_CharmmDcd::WriteHelp() { "\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" + "\tt0 : Initial step.\n" ); } @@ -596,6 +609,8 @@ int Traj_CharmmDcd::processWriteArgs(ArgList& argIn, DataSetList const& DSLin) { else isVel_ = false; timeStep_ = argIn.getKeyDouble("dt", 0.001); + stepsBetweenFrames_ = argIn.getKeyInt("nstep", 1); + initialStep_ = argIn.getKeyInt("t0", 1); return 0; } @@ -692,9 +707,9 @@ 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; diff --git a/src/Traj_CharmmDcd.h b/src/Traj_CharmmDcd.h index f10d374241..1ec8d7fada 100644 --- a/src/Traj_CharmmDcd.h +++ b/src/Traj_CharmmDcd.h @@ -33,6 +33,8 @@ class Traj_CharmmDcd : public TrajectoryIO { 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); From f3b3f503409fa9ea362fabdd1cbce98acd73fed1 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 10 May 2023 12:22:38 -0400 Subject: [PATCH 15/16] Add options to manual. Change t0 to step0 --- doc/cpptraj.lyx | 8 ++++++++ src/Traj_CharmmDcd.cpp | 4 ++-- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 015ea76b3c..9d59717d98 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -18092,6 +18092,14 @@ veltraj Write velocity trajectory instead of coordinates. 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 diff --git a/src/Traj_CharmmDcd.cpp b/src/Traj_CharmmDcd.cpp index b79eedf761..52f034985f 100644 --- a/src/Traj_CharmmDcd.cpp +++ b/src/Traj_CharmmDcd.cpp @@ -586,7 +586,7 @@ void Traj_CharmmDcd::WriteHelp() { "\tveltraj : Write velocity trajectory instead of coordinates.\n" "\tdt : Set trajectory time step in ps.\n" "\tnstep : # steps between frames.\n" - "\tt0 : Initial step.\n" + "\tstep0 : Initial step.\n" ); } @@ -610,7 +610,7 @@ int Traj_CharmmDcd::processWriteArgs(ArgList& argIn, DataSetList const& DSLin) { isVel_ = false; timeStep_ = argIn.getKeyDouble("dt", 0.001); stepsBetweenFrames_ = argIn.getKeyInt("nstep", 1); - initialStep_ = argIn.getKeyInt("t0", 1); + initialStep_ = argIn.getKeyInt("step0", 1); return 0; } From eb3ab08a352a77e39001a8bf58111544a8bd7b5a Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 10 May 2023 12:24:59 -0400 Subject: [PATCH 16/16] Get rid of old unused function definitions --- src/ImageRoutines.h | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/ImageRoutines.h b/src/ImageRoutines.h index ce7362463e..4feda053e3 100644 --- a/src/ImageRoutines.h +++ b/src/ImageRoutines.h @@ -38,12 +38,6 @@ namespace Image { void UnwrapFrac(std::vector&, Frame&, List const&, Matrix_3x3 const&, Matrix_3x3 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& ); - /// Wrap selected atom coords from given frame into primary cell, store in result. void WrapToCell0(std::vector&, Frame const&, AtomMask const&, Matrix_3x3 const&, Matrix_3x3 const&);