From 953879da98bad09e2de57fac9cb4c552c3eddc49 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 7 Feb 2019 15:19:22 -0500 Subject: [PATCH 01/36] DRR - Cpptraj: Add code for calc C6 params for the LJ pme --- src/Ewald.cpp | 11 +++++++++++ src/Ewald.h | 3 +++ src/Ewald_ParticleMesh.cpp | 2 ++ 3 files changed, 16 insertions(+) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index 50537f115a..ce80feff70 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -177,6 +177,17 @@ void Ewald::CalculateCharges(Topology const& topIn, AtomMask const& maskIn) { Setup_VDW_Correction( topIn, maskIn ); } +void Ewald::CalculateC6params(Topology const& topIn, AtomMask const& maskIn) { + Cparam_.clear(); + for (AtomMask::const_iterator atom = maskIn.begin(); atom != maskIn.end(); ++atom) + { + double rmin = topIn.GetVDWradius( *atom ); + double eps = topIn.GetVDWdepth( *atom ); + Cparam_.push_back( 8.0 * (rmin*rmin*rmin) * sqrt(2 * eps) ); + mprintf("DEBUG: C6 param atom %8i = %16.8f\n", *atom+1, Cparam_.back()); + } +} + /** Set up exclusion lists for selected atoms. */ void Ewald::SetupExcluded(Topology const& topIn, AtomMask const& maskIn) { diff --git a/src/Ewald.h b/src/Ewald.h index f4458c7974..2261644e8f 100644 --- a/src/Ewald.h +++ b/src/Ewald.h @@ -45,6 +45,8 @@ class Ewald { int Setup_Pairlist(Box const&, Vec3 const&, double); /// Calculate sum q, sum q^2. Calls setup for vdw correction void CalculateCharges(Topology const&, AtomMask const&); + /// Calculate VDW C6 parameters for LJ PME + void CalculateC6params(Topology const&, AtomMask const&); /// Setup main excluded atom list void SetupExcluded(Topology const&, AtomMask const&); /// Setup VDW correction for selected atom types @@ -68,6 +70,7 @@ class Ewald { Varray Cells_; ///< Hold fractional translations to neighbor cells (non-pairlist only) # endif Darray Charge_; ///< Hold selected atomic charges converted to Amber units. + Darray Cparam_; ///< Hold selected atomic C6 coefficients for LJ PME PairList pairList_; ///< Atom pair list for direct sum. Darray erfc_table_; ///< Hold Erfc cubic spline Y values and coefficients (Y B C D). Iarray2D Excluded_; ///< Full exclusion list for each selected atom. diff --git a/src/Ewald_ParticleMesh.cpp b/src/Ewald_ParticleMesh.cpp index e5ce9eca5f..e0a501e5ef 100644 --- a/src/Ewald_ParticleMesh.cpp +++ b/src/Ewald_ParticleMesh.cpp @@ -127,6 +127,8 @@ int Ewald_ParticleMesh::Init(Box const& boxIn, double cutoffIn, double dsumTolIn /** Setup PME calculation. */ int Ewald_ParticleMesh::Setup(Topology const& topIn, AtomMask const& maskIn) { CalculateCharges(topIn, maskIn); + // TODO put branch here + CalculateC6params( topIn, maskIn ); coordsD_.clear(); coordsD_.reserve( maskIn.Nselected() * 3); SetupExcluded(topIn, maskIn); From 6dd4d6193fe44fd893907fa7c319af762b2dd001 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 7 Feb 2019 15:27:13 -0500 Subject: [PATCH 02/36] DRR - Cpptraj: Add test LJ pme calc. --- src/Action_Energy.cpp | 17 +++++++++-------- src/Ewald_ParticleMesh.cpp | 10 ++++++++++ 2 files changed, 19 insertions(+), 8 deletions(-) diff --git a/src/Action_Energy.cpp b/src/Action_Energy.cpp index 5b6411f434..38983e8044 100644 --- a/src/Action_Energy.cpp +++ b/src/Action_Energy.cpp @@ -332,6 +332,14 @@ Action::RetType Action_Energy::Setup(ActionSetup& setup) { return Action::ERR; EW_->Setup( setup.Top(), Imask_ ); } +# ifdef LIBPME + else if (elecType_ == PME) { + if (((Ewald_ParticleMesh*)EW_)->Init(setup.CoordInfo().TrajBox(), cutoff_, dsumtol_, + ewcoeff_, skinnb_, erfcDx_, npoints_, debug_, mlimits_)) + return Action::ERR; + EW_->Setup( setup.Top(), Imask_ ); + } +# endif // For KE, check for velocities/forces if (KEtype_ != KE_NONE) { if (!setup.CoordInfo().HasVel()) { @@ -348,14 +356,7 @@ Action::RetType Action_Energy::Setup(ActionSetup& setup) { "Warning: 'ketype vv' to estimate kinetic energy.\n"); } } -# ifdef LIBPME - else if (elecType_ == PME) { - if (((Ewald_ParticleMesh*)EW_)->Init(setup.CoordInfo().TrajBox(), cutoff_, dsumtol_, - ewcoeff_, skinnb_, erfcDx_, npoints_, debug_, mlimits_)) - return Action::ERR; - EW_->Setup( setup.Top(), Imask_ ); - } -# endif + currentParm_ = setup.TopAddress(); return Action::OK; } diff --git a/src/Ewald_ParticleMesh.cpp b/src/Ewald_ParticleMesh.cpp index e0a501e5ef..b660e9915d 100644 --- a/src/Ewald_ParticleMesh.cpp +++ b/src/Ewald_ParticleMesh.cpp @@ -183,6 +183,16 @@ double Ewald_ParticleMesh::Recip_ParticleMesh(Box const& boxIn) boxIn.Alpha(), boxIn.Beta(), boxIn.Gamma(), PMEInstanceD::LatticeType::XAligned); double erecip = pme_object->computeERec(0, chargesD, coordsD); + + // LJ PME + Mat cparamD(&Cparam_[0], Cparam_.size(), 1); + auto pme_vdw = std::unique_ptr(new PMEInstanceD()); + pme_vdw->setup(6, ew_coeff_, order_, nfft1, nfft2, nfft3, -1.0, 0); + pme_vdw->setLatticeVectors(boxIn.BoxX(), boxIn.BoxY(), boxIn.BoxZ(), + boxIn.Alpha(), boxIn.Beta(), boxIn.Gamma(), + PMEInstanceD::LatticeType::XAligned); + double evdwrecip = pme_vdw->computeERec(0, cparamD, coordsD); + mprintf("DEBUG: Evdwrecip = %16.8f\n", evdwrecip); t_recip_.Stop(); return erecip; } From 4f4796f84d32732deded1a5e2a12e9ab58b9cc12 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 7 Feb 2019 15:27:56 -0500 Subject: [PATCH 03/36] DRR - Cpptraj: Start LJ PME test. --- test/Test_LJ_PME/RunTest.sh | 22 ++++ test/Test_LJ_PME/water_2.crd | 9 ++ test/Test_LJ_PME/water_2.parm7 | 195 +++++++++++++++++++++++++++++++++ 3 files changed, 226 insertions(+) create mode 100755 test/Test_LJ_PME/RunTest.sh create mode 100644 test/Test_LJ_PME/water_2.crd create mode 100644 test/Test_LJ_PME/water_2.parm7 diff --git a/test/Test_LJ_PME/RunTest.sh b/test/Test_LJ_PME/RunTest.sh new file mode 100755 index 0000000000..c82dae1514 --- /dev/null +++ b/test/Test_LJ_PME/RunTest.sh @@ -0,0 +1,22 @@ +#!/bin/bash + +. ../MasterTest.sh + +CleanFiles lj.in + +INPUT='-i lj.in' +TESTNAME='LJ PME tests.' +Requires libpme maxthreads 1 + +cat > lj.in < Date: Thu, 7 Feb 2019 15:38:35 -0500 Subject: [PATCH 04/36] DRR - Cpptraj: add lj pme self correction --- src/Ewald.cpp | 11 +++++++++++ src/Ewald.h | 4 +++- src/Ewald_ParticleMesh.cpp | 3 +++ 3 files changed, 17 insertions(+), 1 deletion(-) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index ce80feff70..3edbd38d99 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -304,6 +304,17 @@ double Ewald::Self(double volume) { return ene; } +double Ewald::Self6() { + t_self_.Start(); // TODO precalc + double ew2 = ew_coeff_ * ew_coeff_; + double ew6 = ew2 * ew2 * ew2; + double c6sum = 0.0; + for (Darray::const_iterator it = Cparam_.begin(); it != Cparam_.end(); ++it) + c6sum += ew6 * (*it * *it); + t_self_.Stop(); + return c6sum / 12.0; +} + // Ewald::Adjust() # ifdef _OPENMP double Ewald::Adjust(double q0, double q1, double rij) const { diff --git a/src/Ewald.h b/src/Ewald.h index 2261644e8f..2652cdacff 100644 --- a/src/Ewald.h +++ b/src/Ewald.h @@ -37,6 +37,8 @@ class Ewald { inline double ERFC(double) const; /// Ewald "self" energy double Self(double); + /// Ewald "self" energy for C6 term + double Self6(); /// Get analytical estimate of energy due to dispersion interactions > cutoff double Vdw_Correction(double); /// Box, debug, cutoff, dsum tol, ew coeff, erfc dx, nb skin @@ -80,7 +82,7 @@ class Ewald { static const double INVSQRTPI_; double sumq_; ///< Sum of charges double sumq2_; ///< Sum of charges squared - double ew_coeff_; ///< Ewald coefficient + double ew_coeff_; ///< Ewald coefficient TODO separate one for dispersion. double cutoff_; ///< Direct space cutoff double dsumTol_; ///< Direct space sum tolerance. double erfcTableDx_; ///< Spacing of X values in Erfc table. diff --git a/src/Ewald_ParticleMesh.cpp b/src/Ewald_ParticleMesh.cpp index b660e9915d..ea9ef012b1 100644 --- a/src/Ewald_ParticleMesh.cpp +++ b/src/Ewald_ParticleMesh.cpp @@ -205,6 +205,9 @@ double Ewald_ParticleMesh::CalcEnergy(Frame const& frameIn, AtomMask const& mask double volume = frameIn.BoxCrd().ToRecip(ucell, recip); double e_self = Self( volume ); double e_vdwr = Vdw_Correction( volume ); + // TODO branch + double e_vdw6self = Self6(); + mprintf("DEBUG: e_vdw6self = %16.8f\n", e_vdw6self); pairList_.CreatePairList(frameIn, ucell, recip, maskIn); From 93bf88e7534fe481a3498155c7c240be499af03b Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 7 Feb 2019 16:04:41 -0500 Subject: [PATCH 05/36] DRR - Cpptraj: Attempt correction. Not there yet. --- src/Ewald.cpp | 42 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 40 insertions(+), 2 deletions(-) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index 3edbd38d99..cbd829a555 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -347,6 +347,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) double Eelec = 0.0; double e_adjust = 0.0; double Evdw = 0.0; + double Eljpme_correction = 0.0; int cidx; # ifdef _OPENMP # pragma omp parallel private(cidx) reduction(+: Eelec, Evdw, e_adjust) @@ -419,10 +420,28 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) double e_vdw = f12 - f6; // (A/r^12)-(B/r^6) Evdw += e_vdw; //mprintf("PVDW %8i%8i%20.6f%20.6f\n", ta0+1, ta1+1, e_vdw, r2); + // LJ PME direct space correction + double kr2 = ew_coeff_ * ew_coeff_ * rij2; + double kr4 = kr2 * kr2; + //double kr6 = kr2 * kr4; + double expterm = exp(-kr2); + double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; + Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) * r6 * Cij; } } - } else + } else { e_adjust += Adjust(q0, q1, sqrt(rij2)); + // LJ PME direct space correction + // NOTE: Assuming excluded pair is within cutoff + double kr2 = ew_coeff_ * ew_coeff_ * rij2; + double kr4 = kr2 * kr2; + //double kr6 = kr2 * kr4; + double expterm = exp(-kr2); + double r4 = rij2 * rij2; + double r6 = rij2 * r4; + double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; + Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) / r6 * Cij; + } } // END loop over other atoms in thisCell // Loop over all neighbor cells for (unsigned int nidx = 1; nidx != cellList.size(); nidx++) @@ -484,10 +503,28 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) double e_vdw = f12 - f6; // (A/r^12)-(B/r^6) Evdw += e_vdw; //mprintf("PVDW %8i%8i%20.6f%20.6f\n", ta0+1, ta1+1, e_vdw, r2); + // LJ PME direct space correction + double kr2 = ew_coeff_ * ew_coeff_ * rij2; + double kr4 = kr2 * kr2; + //double kr6 = kr2 * kr4; + double expterm = exp(-kr2); + double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; + Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) * r6 * Cij; } } - } else + } else { e_adjust += Adjust(q0, q1, sqrt(rij2)); + // LJ PME direct space correction + // NOTE: Assuming excluded pair is within cutoff + double kr2 = ew_coeff_ * ew_coeff_ * rij2; + double kr4 = kr2 * kr2; + //double kr6 = kr2 * kr4; + double expterm = exp(-kr2); + double r4 = rij2 * rij2; + double r6 = rij2 * r4; + double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; + Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) / r6 * Cij; + } } // END loop over neighbor cell atoms } // END Loop over neighbor cells } // Loop over thisCell atoms @@ -499,6 +536,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) t_direct_.Stop(); e_adjust_out = e_adjust; evdw_out = Evdw; + mprintf("DEBUG: LJ vdw correction = %16.8f\n", Eljpme_correction); return Eelec; } From da7a4212b6b23e255d1824a2981f574deec2093d Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 08:40:04 -0500 Subject: [PATCH 06/36] DRR - Cpptraj: Split up LJ pme correction into excluded and non-excluded parts --- src/Ewald.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index cbd829a555..ec8b6968ca 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -348,6 +348,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) double e_adjust = 0.0; double Evdw = 0.0; double Eljpme_correction = 0.0; + double Eljpme_correction_excl = 0.0; int cidx; # ifdef _OPENMP # pragma omp parallel private(cidx) reduction(+: Eelec, Evdw, e_adjust) @@ -440,7 +441,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) double r4 = rij2 * rij2; double r6 = rij2 * r4; double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; - Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) / r6 * Cij; + Eljpme_correction_excl += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) / r6 * Cij; } } // END loop over other atoms in thisCell // Loop over all neighbor cells @@ -523,7 +524,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) double r4 = rij2 * rij2; double r6 = rij2 * r4; double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; - Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) / r6 * Cij; + Eljpme_correction_excl += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) / r6 * Cij; } } // END loop over neighbor cell atoms } // END Loop over neighbor cells @@ -536,7 +537,8 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) t_direct_.Stop(); e_adjust_out = e_adjust; evdw_out = Evdw; - mprintf("DEBUG: LJ vdw correction = %16.8f\n", Eljpme_correction); + mprintf("DEBUG: LJ vdw correction = %16.8f\n", Eljpme_correction); + mprintf("DEBUG: LJ vdw correction (excluded) = %16.8f\n", Eljpme_correction_excl); return Eelec; } From 3b2458f9feb1d7482d2ca1d2dd02ba022eeacffc Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 08:54:28 -0500 Subject: [PATCH 07/36] DRR - Cpptraj: Add code for switch function, not yet enabled. --- src/Ewald.cpp | 32 +++++++++++++++++++++++++++++--- 1 file changed, 29 insertions(+), 3 deletions(-) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index ec8b6968ca..77ba3dfc64 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -335,6 +335,19 @@ double Ewald::Adjust(double q0, double q1, double rij) { } # endif +static inline double switch_fn(double rij2, double cut2_0, double cut2_1) +{ + if (rij2 <= cut2_0) + return 1.0; + else if (rij2 > cut2_1) + return 0.0; + else { + double xoff_m_x = cut2_1 - rij2; + double fac = 1.0 / (cut2_1 - cut2_0); + return (xoff_m_x*xoff_m_x) * (cut2_1 + 2.0*rij2 - 3.0*cut2_0) * (fac*fac*fac); + } +} + // Ewald::Direct() /** Calculate direct space energy. This is the faster version that uses * a pair list. Also calculate the energy adjustment for excluded @@ -350,8 +363,19 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) double Eljpme_correction = 0.0; double Eljpme_correction_excl = 0.0; int cidx; + + double vswitch = 1.0; + bool use_switch = false; // TODO enable? maybe + double window_width; + if (use_switch) + window_width = 2.0; + else + window_width = 0.0; + double cut0 = cutoff_ - window_width; // TODO check for negative + double cut2_0 = cut0 * cut0; + double cut2_1 = cut2; # ifdef _OPENMP -# pragma omp parallel private(cidx) reduction(+: Eelec, Evdw, e_adjust) +# pragma omp parallel private(cidx, vswitch) reduction(+: Eelec, Evdw, e_adjust,Eljpme_correction,Eljpme_correction_excl ) { # pragma omp for # endif @@ -427,7 +451,8 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) //double kr6 = kr2 * kr4; double expterm = exp(-kr2); double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; - Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) * r6 * Cij; + vswitch = switch_fn(rij2, cut2_0, cut2_1); + Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) * r6 * vswitch * Cij; } } } else { @@ -510,7 +535,8 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) //double kr6 = kr2 * kr4; double expterm = exp(-kr2); double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; - Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) * r6 * Cij; + vswitch = switch_fn(rij2, cut2_0, cut2_1); + Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) * r6 * vswitch * Cij; } } } else { From dd4f0fdca17f13eeea8589f61b132cb394142688 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 09:45:27 -0500 Subject: [PATCH 08/36] DRR - Cpptraj: Add separate terms for LJ PME. Ensure LJ pme is off when LJ Ewald parameter is negative. --- src/Ewald.cpp | 49 +++++++++++++++++++++++++------------- src/Ewald.h | 8 ++++--- src/Ewald_ParticleMesh.cpp | 44 +++++++++++++++++++++++++++------- src/Ewald_ParticleMesh.h | 2 ++ src/Ewald_Regular.cpp | 2 +- 5 files changed, 76 insertions(+), 29 deletions(-) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index 77ba3dfc64..b0272166ca 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -13,6 +13,8 @@ Ewald::Ewald() : sumq_(0.0), sumq2_(0.0), ew_coeff_(0.0), + lw_coeff_(0.0), + switch_width_(0.0), cutoff_(0.0), dsumTol_(0.0), erfcTableDx_(0.0), @@ -230,12 +232,15 @@ void Ewald::SetupExcluded(Topology const& topIn, AtomMask const& maskIn) /** Check some common input. */ int Ewald::CheckInput(Box const& boxIn, int debugIn, double cutoffIn, double dsumTolIn, - double ew_coeffIn, double erfcTableDxIn, double skinnbIn) + double ew_coeffIn, double lw_coeffIn, double switch_widthIn, + double erfcTableDxIn, double skinnbIn) { debug_ = debugIn; cutoff_ = cutoffIn; dsumTol_ = dsumTolIn; ew_coeff_ = ew_coeffIn; + lw_coeff_ = lw_coeffIn; + switch_width_ = switch_widthIn; erfcTableDx_ = erfcTableDxIn; // Check input if (cutoff_ < Constants::SMALL) { @@ -254,6 +259,11 @@ int Ewald::CheckInput(Box const& boxIn, int debugIn, double cutoffIn, double dsu mprinterr("Error: skinnb is less than 0.0\n"); return 1; } + if (switch_width_ < 0.0) switch_width_ = 0.0; + if (switch_width_ > cutoff_) { + mprinterr("Error: Switch width must be less than the cutoff.\n"); + return 1; + } // Set defaults if necessary if (dsumTol_ < Constants::SMALL) @@ -262,7 +272,13 @@ int Ewald::CheckInput(Box const& boxIn, int debugIn, double cutoffIn, double dsu ew_coeff_ = FindEwaldCoefficient( cutoff_, dsumTol_ ); if (erfcTableDx_ <= 0.0) erfcTableDx_ = 1.0 / 5000; // TODO make this optional - FillErfcTable( cutoff_, ew_coeff_ ); + FillErfcTable( cutoff_, ew_coeff_ ); + // TODO do for C6 as well + // TODO for C6 correction term + if (lw_coeff_ < 0.0) + lw_coeff_ = 0.0; + else if (DABS(lw_coeff_) < Constants::SMALL) + lw_coeff_ = ew_coeff_; return 0; } @@ -291,7 +307,7 @@ int Ewald::Setup_Pairlist(Box const& boxIn, Vec3 const& recipLengths, double ski return 0; } -/** Self energy. This is the cancelling Gaussian plus the "neutralizing plasma". */ +/** Electrostatic self energy. This is the cancelling Gaussian plus the "neutralizing plasma". */ double Ewald::Self(double volume) { t_self_.Start(); double d0 = -ew_coeff_ * INVSQRTPI_; @@ -304,9 +320,10 @@ double Ewald::Self(double volume) { return ene; } +/** Lennard-Jones self energy. */ double Ewald::Self6() { t_self_.Start(); // TODO precalc - double ew2 = ew_coeff_ * ew_coeff_; + double ew2 = lw_coeff_ * lw_coeff_; double ew6 = ew2 * ew2 * ew2; double c6sum = 0.0; for (Darray::const_iterator it = Cparam_.begin(); it != Cparam_.end(); ++it) @@ -335,6 +352,7 @@ double Ewald::Adjust(double q0, double q1, double rij) { } # endif +/** Switching function for Lennard-Jones. */ static inline double switch_fn(double rij2, double cut2_0, double cut2_1) { if (rij2 <= cut2_0) @@ -352,6 +370,9 @@ static inline double switch_fn(double rij2, double cut2_0, double cut2_1) /** Calculate direct space energy. This is the faster version that uses * a pair list. Also calculate the energy adjustment for excluded * atoms. + * \param PL The pairlist used to calculate energy. + * \param e_adjust_out The electrostatic adjust energy for excluded atoms. + * \param evdw_out The direct space van der Waals term */ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) { @@ -364,14 +385,8 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) double Eljpme_correction_excl = 0.0; int cidx; - double vswitch = 1.0; - bool use_switch = false; // TODO enable? maybe - double window_width; - if (use_switch) - window_width = 2.0; - else - window_width = 0.0; - double cut0 = cutoff_ - window_width; // TODO check for negative + double vswitch; + double cut0 = cutoff_ - switch_width_; double cut2_0 = cut0 * cut0; double cut2_1 = cut2; # ifdef _OPENMP @@ -446,7 +461,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) Evdw += e_vdw; //mprintf("PVDW %8i%8i%20.6f%20.6f\n", ta0+1, ta1+1, e_vdw, r2); // LJ PME direct space correction - double kr2 = ew_coeff_ * ew_coeff_ * rij2; + double kr2 = lw_coeff_ * lw_coeff_ * rij2; double kr4 = kr2 * kr2; //double kr6 = kr2 * kr4; double expterm = exp(-kr2); @@ -459,7 +474,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) e_adjust += Adjust(q0, q1, sqrt(rij2)); // LJ PME direct space correction // NOTE: Assuming excluded pair is within cutoff - double kr2 = ew_coeff_ * ew_coeff_ * rij2; + double kr2 = lw_coeff_ * lw_coeff_ * rij2; double kr4 = kr2 * kr2; //double kr6 = kr2 * kr4; double expterm = exp(-kr2); @@ -530,7 +545,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) Evdw += e_vdw; //mprintf("PVDW %8i%8i%20.6f%20.6f\n", ta0+1, ta1+1, e_vdw, r2); // LJ PME direct space correction - double kr2 = ew_coeff_ * ew_coeff_ * rij2; + double kr2 = lw_coeff_ * lw_coeff_ * rij2; double kr4 = kr2 * kr2; //double kr6 = kr2 * kr4; double expterm = exp(-kr2); @@ -543,7 +558,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) e_adjust += Adjust(q0, q1, sqrt(rij2)); // LJ PME direct space correction // NOTE: Assuming excluded pair is within cutoff - double kr2 = ew_coeff_ * ew_coeff_ * rij2; + double kr2 = lw_coeff_ * lw_coeff_ * rij2; double kr4 = kr2 * kr2; //double kr6 = kr2 * kr4; double expterm = exp(-kr2); @@ -562,7 +577,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) # endif t_direct_.Stop(); e_adjust_out = e_adjust; - evdw_out = Evdw; + evdw_out = Evdw + Eljpme_correction + Eljpme_correction_excl; mprintf("DEBUG: LJ vdw correction = %16.8f\n", Eljpme_correction); mprintf("DEBUG: LJ vdw correction (excluded) = %16.8f\n", Eljpme_correction_excl); return Eelec; diff --git a/src/Ewald.h b/src/Ewald.h index 2652cdacff..ca6575ded5 100644 --- a/src/Ewald.h +++ b/src/Ewald.h @@ -41,8 +41,8 @@ class Ewald { double Self6(); /// Get analytical estimate of energy due to dispersion interactions > cutoff double Vdw_Correction(double); - /// Box, debug, cutoff, dsum tol, ew coeff, erfc dx, nb skin - int CheckInput(Box const&, int, double, double, double, double, double); + /// Box, debug, cutoff, dsum tol, ew coeff, lj coeff, switch window, erfc dx, nb skin + int CheckInput(Box const&, int, double, double, double, double, double, double, double); /// Set up pair list int Setup_Pairlist(Box const&, Vec3 const&, double); /// Calculate sum q, sum q^2. Calls setup for vdw correction @@ -82,7 +82,9 @@ class Ewald { static const double INVSQRTPI_; double sumq_; ///< Sum of charges double sumq2_; ///< Sum of charges squared - double ew_coeff_; ///< Ewald coefficient TODO separate one for dispersion. + double ew_coeff_; ///< Ewald coefficient for electrostatics + double lw_coeff_; ///< Ewald coefficient for LJ + double switch_width_; ///< Switching window size for LJ switch if active double cutoff_; ///< Direct space cutoff double dsumTol_; ///< Direct space sum tolerance. double erfcTableDx_; ///< Spacing of X values in Erfc table. diff --git a/src/Ewald_ParticleMesh.cpp b/src/Ewald_ParticleMesh.cpp index ea9ef012b1..d7e9799233 100644 --- a/src/Ewald_ParticleMesh.cpp +++ b/src/Ewald_ParticleMesh.cpp @@ -91,7 +91,7 @@ int Ewald_ParticleMesh::Init(Box const& boxIn, double cutoffIn, double dsumTolIn double ew_coeffIn, double skinnbIn, double erfcTableDxIn, int orderIn, int debugIn, const int* nfftIn) { - if (CheckInput(boxIn, debugIn, cutoffIn, dsumTolIn, ew_coeffIn, erfcTableDxIn, skinnbIn)) + if (CheckInput(boxIn, debugIn, cutoffIn, dsumTolIn, ew_coeffIn, -1.0, 0.0, erfcTableDxIn, skinnbIn)) return 1; if (nfftIn != 0) std::copy(nfftIn, nfftIn+3, nfft_); @@ -184,8 +184,25 @@ double Ewald_ParticleMesh::Recip_ParticleMesh(Box const& boxIn) PMEInstanceD::LatticeType::XAligned); double erecip = pme_object->computeERec(0, chargesD, coordsD); - // LJ PME + t_recip_.Stop(); + return erecip; +} + +/** The LJ PME reciprocal term. */ +double Ewald_ParticleMesh::LJ_Recip_ParticleMesh(Box const& boxIn) +{ + t_recip_.Start(); + int nfft1 = nfft_[0]; + int nfft2 = nfft_[1]; + int nfft3 = nfft_[2]; + if ( DetermineNfft(nfft1, nfft2, nfft3, boxIn) ) { + mprinterr("Error: Could not determine grid spacing.\n"); + return 0.0; + } + + Mat coordsD(&coordsD_[0], Charge_.size(), 3); Mat cparamD(&Cparam_[0], Cparam_.size(), 1); + auto pme_vdw = std::unique_ptr(new PMEInstanceD()); pme_vdw->setup(6, ew_coeff_, order_, nfft1, nfft2, nfft3, -1.0, 0); pme_vdw->setLatticeVectors(boxIn.BoxX(), boxIn.BoxY(), boxIn.BoxZ(), @@ -194,7 +211,7 @@ double Ewald_ParticleMesh::Recip_ParticleMesh(Box const& boxIn) double evdwrecip = pme_vdw->computeERec(0, cparamD, coordsD); mprintf("DEBUG: Evdwrecip = %16.8f\n", evdwrecip); t_recip_.Stop(); - return erecip; + return evdwrecip; } /** Calculate full nonbonded energy with PME */ @@ -204,10 +221,7 @@ double Ewald_ParticleMesh::CalcEnergy(Frame const& frameIn, AtomMask const& mask Matrix_3x3 ucell, recip; double volume = frameIn.BoxCrd().ToRecip(ucell, recip); double e_self = Self( volume ); - double e_vdwr = Vdw_Correction( volume ); - // TODO branch - double e_vdw6self = Self6(); - mprintf("DEBUG: e_vdw6self = %16.8f\n", e_vdw6self); + double e_vdw_lr_correction; pairList_.CreatePairList(frameIn, ucell, recip, maskIn); @@ -223,13 +237,27 @@ double Ewald_ParticleMesh::CalcEnergy(Frame const& frameIn, AtomMask const& mask // MapCoords(frameIn, ucell, recip, maskIn); double e_recip = Recip_ParticleMesh( frameIn.BoxCrd() ); + + // TODO branch + double e_vdw6self, e_vdw6recip; + if (lw_coeff_ > 0.0) { + e_vdw6self = Self6(); + mprintf("DEBUG: e_vdw6self = %16.8f\n", e_vdw6self); + e_vdw6recip = LJ_Recip_ParticleMesh( frameIn.BoxCrd() ); + e_vdw_lr_correction = 0.0; + } else { + e_vdw6self = 0.0; + e_vdw6recip = 0.0; + e_vdw_lr_correction = Vdw_Correction( volume ); + } + double e_adjust = 0.0; e_vdw = 0.0; double e_direct = Direct( pairList_, e_adjust, e_vdw ); if (debug_ > 0) mprintf("DEBUG: Eself= %20.10f Erecip= %20.10f Edirect= %20.10f Eadjust= %20.10f Evdw= %20.10f\n", e_self, e_recip, e_direct, e_adjust, e_vdw); - e_vdw += e_vdwr; + e_vdw += (e_vdw_lr_correction + e_vdw6self + e_vdw6recip); t_total_.Stop(); return e_self + e_recip + e_direct + e_adjust; } diff --git a/src/Ewald_ParticleMesh.h b/src/Ewald_ParticleMesh.h index 73c625bbf7..2c56102d54 100644 --- a/src/Ewald_ParticleMesh.h +++ b/src/Ewald_ParticleMesh.h @@ -19,6 +19,8 @@ class Ewald_ParticleMesh : public Ewald { int DetermineNfft(int&, int&, int&, Box const&) const; /// Particle mesh Ewald reciprocal energy double Recip_ParticleMesh(Box const&); + /// Particle mesh Ewald LJ recip energy + double LJ_Recip_ParticleMesh(Box const&); Darray coordsD_; ///< Hold coordinates for selected atoms diff --git a/src/Ewald_Regular.cpp b/src/Ewald_Regular.cpp index 0a713040a6..5073c25207 100644 --- a/src/Ewald_Regular.cpp +++ b/src/Ewald_Regular.cpp @@ -100,7 +100,7 @@ int Ewald_Regular::Init(Box const& boxIn, double cutoffIn, double dsumTolIn, dou double ew_coeffIn, double maxexpIn, double skinnbIn, double erfcTableDxIn, int debugIn, const int* mlimitsIn) { - if (CheckInput(boxIn, debugIn, cutoffIn, dsumTolIn, ew_coeffIn, erfcTableDxIn, skinnbIn)) + if (CheckInput(boxIn, debugIn, cutoffIn, dsumTolIn, ew_coeffIn, -1.0, 0.0, erfcTableDxIn, skinnbIn)) return 1; rsumTol_ = rsumTolIn; maxexp_ = maxexpIn; From cafd687d208efeaf9c8b72bb27657f43720082d7 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 10:06:46 -0500 Subject: [PATCH 09/36] DRR - Cpptraj: Add function for calculating sigma --- src/Topology.cpp | 12 ++++++++++++ src/Topology.h | 1 + 2 files changed, 13 insertions(+) diff --git a/src/Topology.cpp b/src/Topology.cpp index 58da874c03..5ee8816d70 100644 --- a/src/Topology.cpp +++ b/src/Topology.cpp @@ -491,6 +491,7 @@ void Topology::Resize(Pointers const& pIn) { dihedralparm_.resize( pIn.nDihParm_ ); } +/** \return Rmin for given atom. */ double Topology::GetVDWradius(int a1) const { //TODO: return zero when no params? NonbondType const& LJ = GetLJparam(a1, a1); @@ -500,6 +501,17 @@ double Topology::GetVDWradius(int a1) const { return 0.0; } +/** \return sigma for given atom. */ +double Topology::GetVDWsigma(int a1) const { + //TODO: return zero when no params? + NonbondType const& LJ = GetLJparam(a1, a1); + if (LJ.B() > 0.0) + return ( 0.5 * pow(LJ.A() / LJ.B(), (1.0/6.0)) ); + else + return 0.0; +} + +/** \return epsilon for given atom. */ double Topology::GetVDWdepth(int a1) const { NonbondType const& LJ = GetLJparam(a1, a1); if (LJ.A() > 0.0) diff --git a/src/Topology.h b/src/Topology.h index e05cb94058..a60350ec79 100644 --- a/src/Topology.h +++ b/src/Topology.h @@ -101,6 +101,7 @@ class Topology { NonbondParmType const& Nonbond() const { return nonbond_; } NonbondParmType& SetNonbond() { return nonbond_; } double GetVDWradius(int) const; + double GetVDWsigma(int) const; double GetVDWdepth(int) const; /// \return Lennard-Jones 6-12 parameters for given pair of atoms inline NonbondType const& GetLJparam(int, int) const; From 0d9d8cfd57cde741a8e7ea755120a37df8e278ff Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 10:09:52 -0500 Subject: [PATCH 10/36] DRR - Cpptraj: Add keywords for lennard jones pme --- src/Action_Energy.cpp | 24 ++++++++++++++++++++---- src/Action_Energy.h | 1 + src/Ewald_ParticleMesh.cpp | 4 ++-- src/Ewald_ParticleMesh.h | 4 ++-- 4 files changed, 25 insertions(+), 8 deletions(-) diff --git a/src/Action_Energy.cpp b/src/Action_Energy.cpp index 38983e8044..805479a169 100644 --- a/src/Action_Energy.cpp +++ b/src/Action_Energy.cpp @@ -24,7 +24,7 @@ void Action_Energy::Help() const { "\t [ewcoeff ] [maxexp ] [skinnb ]\n" "\t [mlimits ,,] [erfcdx ]\n" "\t pme [cut ] [dsumtol ] [order ]\n" - "\t [ewcoeff ] [skinnb ]\n" + "\t [ewcoeff ] [ljpme] [ewcoefflj] [skinnb ]\n" "\t [nfft ,,] [erfcdx ]\n" "\t } ]\n" " Calculate energy for atoms in mask.\n"); @@ -140,6 +140,11 @@ Action::RetType Action_Energy::Init(ArgList& actionArgs, ActionInit& init, int d cutoff_ = actionArgs.getKeyDouble("cut", 8.0); dsumtol_ = actionArgs.getKeyDouble("dsumtol", 1E-5); ewcoeff_ = actionArgs.getKeyDouble("ewcoeff", 0.0); + lwcoeff_ = -1.0; + if (actionArgs.hasKey("ljpme")) + lwcoeff_ = 0.4; + lwcoeff_ = actionArgs.getKeyDouble("ewcoefflj", lwcoeff_); + // TODO arguments for switching skinnb_ = actionArgs.getKeyDouble("skinnb", 2.0); erfcDx_ = actionArgs.getKeyDouble("erfcdx", 0.0); npoints_ = actionArgs.getKeyInt("order", 6); @@ -221,6 +226,8 @@ Action::RetType Action_Energy::Init(ArgList& actionArgs, ActionInit& init, int d need_lj_params_ = true; } } + if (lj_longrange_correction && lwcoeff_ >= 0.0) + lj_longrange_correction = false; // Get Masks Mask1_.SetMaskString( actionArgs.GetMaskNext() ); @@ -290,8 +297,16 @@ Action::RetType Action_Energy::Init(ArgList& actionArgs, ActionInit& init, int d if (erfcDx_ > 0.0) mprintf("\tERFC table dx= %g\n", erfcDx_); } - if (termEnabled[VDW] && lj_longrange_correction) - mprintf("\tUsing long range correction for nonbond VDW calc.\n"); + if (termEnabled[VDW]) { + if (lj_longrange_correction) + mprintf("\tUsing long range correction for nonbond VDW calc.\n"); + else if (lwcoeff_ >= 0.0) { + if (lwcoeff_ > 0.0) + mprintf("\tUsing Lennard-Jones PME with Ewald coefficient %.4f\n", lwcoeff_); + else + mprintf("\tLennard-Jones PME Ewald coefficient will be set to elec. Ewald coefficient.\n"); + } + } if (KEtype_ != KE_NONE) { if (KEtype_ == KE_AUTO) mprintf("\tIf forces and velocities present KE will be calculated assuming\n" @@ -335,7 +350,8 @@ Action::RetType Action_Energy::Setup(ActionSetup& setup) { # ifdef LIBPME else if (elecType_ == PME) { if (((Ewald_ParticleMesh*)EW_)->Init(setup.CoordInfo().TrajBox(), cutoff_, dsumtol_, - ewcoeff_, skinnb_, erfcDx_, npoints_, debug_, mlimits_)) + ewcoeff_, lwcoeff_, skinnb_, erfcDx_, npoints_, + debug_, mlimits_)) return Action::ERR; EW_->Setup( setup.Top(), Imask_ ); } diff --git a/src/Action_Energy.h b/src/Action_Energy.h index 668a7acdd3..d5520d68f5 100644 --- a/src/Action_Energy.h +++ b/src/Action_Energy.h @@ -45,6 +45,7 @@ class Action_Energy: public Action { double dsumtol_; ///< Ewald direct sum tolerance. double rsumtol_; ///< Regular Ewald reciprocal sum tolerance. double ewcoeff_; ///< Ewald coefficient. + double lwcoeff_; ///< LJ Ewald coefficient. double maxexp_; double skinnb_; ///< Size of non-bonded "skin" double erfcDx_; ///< Spacing for ERFC table (default 1/5000) diff --git a/src/Ewald_ParticleMesh.cpp b/src/Ewald_ParticleMesh.cpp index d7e9799233..282ce3e18a 100644 --- a/src/Ewald_ParticleMesh.cpp +++ b/src/Ewald_ParticleMesh.cpp @@ -88,10 +88,10 @@ int Ewald_ParticleMesh::DetermineNfft(int& nfft1, int& nfft2, int& nfft3, Box co /** Set up PME parameters. */ int Ewald_ParticleMesh::Init(Box const& boxIn, double cutoffIn, double dsumTolIn, - double ew_coeffIn, double skinnbIn, double erfcTableDxIn, + double ew_coeffIn, double lw_coeffIn,double skinnbIn, double erfcTableDxIn, int orderIn, int debugIn, const int* nfftIn) { - if (CheckInput(boxIn, debugIn, cutoffIn, dsumTolIn, ew_coeffIn, -1.0, 0.0, erfcTableDxIn, skinnbIn)) + if (CheckInput(boxIn, debugIn, cutoffIn, dsumTolIn, ew_coeffIn, lw_coeffIn, 0.0, erfcTableDxIn, skinnbIn)) return 1; if (nfftIn != 0) std::copy(nfftIn, nfftIn+3, nfft_); diff --git a/src/Ewald_ParticleMesh.h b/src/Ewald_ParticleMesh.h index 2c56102d54..a8864d9cb5 100644 --- a/src/Ewald_ParticleMesh.h +++ b/src/Ewald_ParticleMesh.h @@ -6,8 +6,8 @@ class Ewald_ParticleMesh : public Ewald { public: Ewald_ParticleMesh(); - /// Box, cut, dsum tol, ew coeff, NB skin, erfc dx, order, debug, nfft - int Init(Box const&, double, double, double, double, double, int, int, const int*); + /// Box, cut, dsum tol, ew coeff, lj ew coeff, NB skin, erfc dx, order, debug, nfft + int Init(Box const&, double, double, double, double, double, double, int, int, const int*); // ----- Inherited --------------------------- int Setup(Topology const&, AtomMask const&); double CalcEnergy(Frame const&, AtomMask const&, double&); // TODO const? From c40fd0d9af08d7385ec06a3a5176c9676c216bfa Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 10:13:19 -0500 Subject: [PATCH 11/36] DRR - Cpptraj: Code doc --- src/Ewald_ParticleMesh.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Ewald_ParticleMesh.cpp b/src/Ewald_ParticleMesh.cpp index 282ce3e18a..0ae1487c26 100644 --- a/src/Ewald_ParticleMesh.cpp +++ b/src/Ewald_ParticleMesh.cpp @@ -127,7 +127,8 @@ int Ewald_ParticleMesh::Init(Box const& boxIn, double cutoffIn, double dsumTolIn /** Setup PME calculation. */ int Ewald_ParticleMesh::Setup(Topology const& topIn, AtomMask const& maskIn) { CalculateCharges(topIn, maskIn); - // TODO put branch here + // NOTE: These dont need to actually be calculated if the lj ewald coeff + // is 0.0, but do it here anyway to avoid segfaults. CalculateC6params( topIn, maskIn ); coordsD_.clear(); coordsD_.reserve( maskIn.Nselected() * 3); From f611ee880b956697191bd1731209fbc0afcb7aec Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 10:13:32 -0500 Subject: [PATCH 12/36] DRR - Cpptraj: Simple LJ pme test --- test/Test_LJ_PME/RunTest.sh | 5 +++-- test/Test_LJ_PME/ene.dat.save | 2 ++ 2 files changed, 5 insertions(+), 2 deletions(-) create mode 100644 test/Test_LJ_PME/ene.dat.save diff --git a/test/Test_LJ_PME/RunTest.sh b/test/Test_LJ_PME/RunTest.sh index c82dae1514..7fb6648ed6 100755 --- a/test/Test_LJ_PME/RunTest.sh +++ b/test/Test_LJ_PME/RunTest.sh @@ -2,7 +2,7 @@ . ../MasterTest.sh -CleanFiles lj.in +CleanFiles lj.in ene.dat INPUT='-i lj.in' TESTNAME='LJ PME tests.' @@ -14,9 +14,10 @@ trajin water_2.crd box x 20 y 20 z 20 alpha 90 beta 90 gamma 90 debug 10 -energy out ene.dat etype pme cut 8.0 dsumtol 0.0000001 skinnb 0.01 ewcoeff 0.4 +energy out ene.dat prec 16.8 etype pme cut 8.0 dsumtol 0.0000001 skinnb 0.01 ewcoeff 0.3 ewcoefflj 0.3 EOF RunCpptraj "LJ PME test." +DoTest ene.dat.save ene.dat EndTest exit 0 diff --git a/test/Test_LJ_PME/ene.dat.save b/test/Test_LJ_PME/ene.dat.save new file mode 100644 index 0000000000..71fe4d54ae --- /dev/null +++ b/test/Test_LJ_PME/ene.dat.save @@ -0,0 +1,2 @@ +#Frame ENE_00001[bond] ENE_00001[angle] ENE_00001[dih] ENE_00001[vdw14] ENE_00001[elec14] ENE_00001[vdw] ENE_00001[elec] ENE_00001[total] + 1 46.56162949 88.49177893 0.00000000 0.00000000 0.00000000 -0.18643240 -0.19203566 134.67494036 From 68487d0eb152785ab458bcaec8c3043ea14dc24d Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 10:23:34 -0500 Subject: [PATCH 13/36] DRR - Cpptraj: Make sure the correct Ewald param passed in for LJ --- src/Ewald_ParticleMesh.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Ewald_ParticleMesh.cpp b/src/Ewald_ParticleMesh.cpp index 0ae1487c26..707027549a 100644 --- a/src/Ewald_ParticleMesh.cpp +++ b/src/Ewald_ParticleMesh.cpp @@ -105,6 +105,8 @@ int Ewald_ParticleMesh::Init(Box const& boxIn, double cutoffIn, double dsumTolIn mprintf("\tParticle Mesh Ewald params:\n"); mprintf("\t Cutoff= %g Direct Sum Tol= %g Ewald coeff.= %g NB skin= %g\n", cutoff_, dsumTol_, ew_coeff_, skinnbIn); + if (lw_coeff_ > 0.0) + mprintf("\t LJ Ewald coeff.= %g\n", lw_coeff_); mprintf("\t Bspline order= %i\n", order_); mprintf("\t Erfc table dx= %g, size= %zu\n", erfcTableDx_, erfc_table_.size()/4); mprintf("\t "); @@ -205,7 +207,7 @@ double Ewald_ParticleMesh::LJ_Recip_ParticleMesh(Box const& boxIn) Mat cparamD(&Cparam_[0], Cparam_.size(), 1); auto pme_vdw = std::unique_ptr(new PMEInstanceD()); - pme_vdw->setup(6, ew_coeff_, order_, nfft1, nfft2, nfft3, -1.0, 0); + pme_vdw->setup(6, lw_coeff_, order_, nfft1, nfft2, nfft3, -1.0, 0); pme_vdw->setLatticeVectors(boxIn.BoxX(), boxIn.BoxY(), boxIn.BoxZ(), boxIn.Alpha(), boxIn.Beta(), boxIn.Gamma(), PMEInstanceD::LatticeType::XAligned); From 7ea3c70b702dae480a169a9fda7f516501acc800 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 10:37:18 -0500 Subject: [PATCH 14/36] DRR - Cpptraj: Fix regular Ewald by allocating C6 arrays even though they arent used --- src/Ewald_Regular.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Ewald_Regular.cpp b/src/Ewald_Regular.cpp index 5073c25207..0f99aee0a7 100644 --- a/src/Ewald_Regular.cpp +++ b/src/Ewald_Regular.cpp @@ -158,6 +158,8 @@ int Ewald_Regular::Init(Box const& boxIn, double cutoffIn, double dsumTolIn, dou /** Setup regular Ewald calculation. */ int Ewald_Regular::Setup(Topology const& topIn, AtomMask const& maskIn) { CalculateCharges(topIn, maskIn); + // Blank C6 Arrays. TODO actually blank them + CalculateC6params( topIn, maskIn ); // Build exponential factors for use in structure factors. // These arrays are laid out in 1D; value for each atom at each m, i.e. From fa46d6c09265dae8a44a99dcb8206accfd855bd1 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 10:48:28 -0500 Subject: [PATCH 15/36] DRR - Cpptraj: Add lj PME kappa sweep test. --- test/Test_LJ_PME/RunTest.sh | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/test/Test_LJ_PME/RunTest.sh b/test/Test_LJ_PME/RunTest.sh index 7fb6648ed6..e94ffbdb35 100755 --- a/test/Test_LJ_PME/RunTest.sh +++ b/test/Test_LJ_PME/RunTest.sh @@ -2,7 +2,7 @@ . ../MasterTest.sh -CleanFiles lj.in ene.dat +CleanFiles lj.in ene.dat ene.dat.? INPUT='-i lj.in' TESTNAME='LJ PME tests.' @@ -13,11 +13,30 @@ parm water_2.parm7 trajin water_2.crd box x 20 y 20 z 20 alpha 90 beta 90 gamma 90 -debug 10 +#debug 10 energy out ene.dat prec 16.8 etype pme cut 8.0 dsumtol 0.0000001 skinnb 0.01 ewcoeff 0.3 ewcoefflj 0.3 EOF RunCpptraj "LJ PME test." DoTest ene.dat.save ene.dat +# Kappa sweep +i=0 +for kappa in '0.25' '0.35' '0.45' '0.5' ; do + cat > lj.in < Date: Thu, 14 Feb 2019 10:54:15 -0500 Subject: [PATCH 16/36] DRR - Cpptraj: Hide debug info. If LJ pme will not be active set C6 params to 0.0 --- src/Ewald.cpp | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index b0272166ca..beea06c355 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -181,13 +181,17 @@ void Ewald::CalculateCharges(Topology const& topIn, AtomMask const& maskIn) { void Ewald::CalculateC6params(Topology const& topIn, AtomMask const& maskIn) { Cparam_.clear(); - for (AtomMask::const_iterator atom = maskIn.begin(); atom != maskIn.end(); ++atom) - { - double rmin = topIn.GetVDWradius( *atom ); - double eps = topIn.GetVDWdepth( *atom ); - Cparam_.push_back( 8.0 * (rmin*rmin*rmin) * sqrt(2 * eps) ); - mprintf("DEBUG: C6 param atom %8i = %16.8f\n", *atom+1, Cparam_.back()); - } + if (lw_coeff_ > 0.0) { + for (AtomMask::const_iterator atom = maskIn.begin(); atom != maskIn.end(); ++atom) + { + double rmin = topIn.GetVDWradius( *atom ); + double eps = topIn.GetVDWdepth( *atom ); + Cparam_.push_back( 8.0 * (rmin*rmin*rmin) * sqrt(2 * eps) ); + if (debug_ > 0) + mprintf("DEBUG: C6 param atom %8i = %16.8f\n", *atom+1, Cparam_.back()); + } + } else + Cparam_.assign(maskIn.Nselected(), 0.0); } /** Set up exclusion lists for selected atoms. */ From 4f5c76144f88fe96d5ee5b85e6770bf3ae686bc4 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 10:59:50 -0500 Subject: [PATCH 17/36] DRR - Cpptraj: Hide some debug info --- src/Ewald.cpp | 7 +++++-- src/Ewald_ParticleMesh.cpp | 6 ++++-- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index beea06c355..196ea1be69 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -581,9 +581,12 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) # endif t_direct_.Stop(); e_adjust_out = e_adjust; +# ifdef DEBUG_PAIRLIST + mprintf("DEBUG: LJ vdw = %16.8f\n", Evdw); + mprintf("DEBUG: LJ vdw PME correction = %16.8f\n", Eljpme_correction); + mprintf("DEBUG: LJ vdw PME correction (excluded) = %16.8f\n", Eljpme_correction_excl); +# endif evdw_out = Evdw + Eljpme_correction + Eljpme_correction_excl; - mprintf("DEBUG: LJ vdw correction = %16.8f\n", Eljpme_correction); - mprintf("DEBUG: LJ vdw correction (excluded) = %16.8f\n", Eljpme_correction_excl); return Eelec; } diff --git a/src/Ewald_ParticleMesh.cpp b/src/Ewald_ParticleMesh.cpp index 707027549a..bf79aacec5 100644 --- a/src/Ewald_ParticleMesh.cpp +++ b/src/Ewald_ParticleMesh.cpp @@ -212,7 +212,6 @@ double Ewald_ParticleMesh::LJ_Recip_ParticleMesh(Box const& boxIn) boxIn.Alpha(), boxIn.Beta(), boxIn.Gamma(), PMEInstanceD::LatticeType::XAligned); double evdwrecip = pme_vdw->computeERec(0, cparamD, coordsD); - mprintf("DEBUG: Evdwrecip = %16.8f\n", evdwrecip); t_recip_.Stop(); return evdwrecip; } @@ -245,8 +244,11 @@ double Ewald_ParticleMesh::CalcEnergy(Frame const& frameIn, AtomMask const& mask double e_vdw6self, e_vdw6recip; if (lw_coeff_ > 0.0) { e_vdw6self = Self6(); - mprintf("DEBUG: e_vdw6self = %16.8f\n", e_vdw6self); e_vdw6recip = LJ_Recip_ParticleMesh( frameIn.BoxCrd() ); + if (debug_ > 0) { + mprintf("DEBUG: e_vdw6self = %16.8f\n", e_vdw6self); + mprintf("DEBUG: Evdwrecip = %16.8f\n", e_vdw6recip); + } e_vdw_lr_correction = 0.0; } else { e_vdw6self = 0.0; From e7e74fd0968ed740ec8c87188f269956a9167ed9 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 11:00:01 -0500 Subject: [PATCH 18/36] DRR - Cpptraj: Enable lj pme test --- test/Makefile | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Makefile b/test/Makefile index f2e7ef3939..88060d3656 100644 --- a/test/Makefile +++ b/test/Makefile @@ -352,6 +352,7 @@ test.gromacs: test.ewald: @-cd Test_Ewald && ./RunTest.sh $(OPT) @-cd Test_LibPME && ./RunTest.sh $(OPT) + @-cd Test_LJ_PME && ./RunTest.sh $(OPT) test.ti: @-cd Test_TI && ./RunTest.sh $(OPT) From 31c0df36d835dbb4f8e8c0afc9b7a6c545ef7bd4 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 11:31:23 -0500 Subject: [PATCH 19/36] DRR - Cpptraj: Box action modifies topology, return value to indicate this. --- src/Action_Box.cpp | 2 +- test/Test_LJ_PME/RunTest.sh | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Action_Box.cpp b/src/Action_Box.cpp index 559c500cc9..5439165c27 100644 --- a/src/Action_Box.cpp +++ b/src/Action_Box.cpp @@ -54,7 +54,7 @@ Action::RetType Action_Box::Setup(ActionSetup& setup) { cInfo_.SetBox( pbox ); } setup.SetCoordInfo( &cInfo_ ); - return Action::OK; + return Action::MODIFY_TOPOLOGY; } Action::RetType Action_Box::DoAction(int frameNum, ActionFrame& frm) { diff --git a/test/Test_LJ_PME/RunTest.sh b/test/Test_LJ_PME/RunTest.sh index e94ffbdb35..41eac4d613 100755 --- a/test/Test_LJ_PME/RunTest.sh +++ b/test/Test_LJ_PME/RunTest.sh @@ -8,6 +8,7 @@ INPUT='-i lj.in' TESTNAME='LJ PME tests.' Requires libpme maxthreads 1 +# Basic test cat > lj.in < Date: Thu, 14 Feb 2019 11:34:41 -0500 Subject: [PATCH 20/36] DRR - Cpptraj: Add switch to inside the direct space calc --- src/Ewald.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index 196ea1be69..8cb8ccb9ef 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -455,6 +455,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) int nbindex = NB_->GetLJindex(TypeIndices_[it0->Idx()], TypeIndices_[it1->Idx()]); if (nbindex > -1) { + vswitch = switch_fn(rij2, cut2_0, cut2_1); NonbondType const& LJ = NB_->NBarray()[ nbindex ]; double r2 = 1.0 / rij2; double r6 = r2 * r2 * r2; @@ -462,7 +463,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) double f12 = LJ.A() * r12; // A/r^12 double f6 = LJ.B() * r6; // B/r^6 double e_vdw = f12 - f6; // (A/r^12)-(B/r^6) - Evdw += e_vdw; + Evdw += (e_vdw * vswitch); //mprintf("PVDW %8i%8i%20.6f%20.6f\n", ta0+1, ta1+1, e_vdw, r2); // LJ PME direct space correction double kr2 = lw_coeff_ * lw_coeff_ * rij2; @@ -470,7 +471,6 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) //double kr6 = kr2 * kr4; double expterm = exp(-kr2); double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; - vswitch = switch_fn(rij2, cut2_0, cut2_1); Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) * r6 * vswitch * Cij; } } @@ -539,6 +539,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) int nbindex = NB_->GetLJindex(TypeIndices_[it0->Idx()], TypeIndices_[it1->Idx()]); if (nbindex > -1) { + vswitch = switch_fn(rij2, cut2_0, cut2_1); NonbondType const& LJ = NB_->NBarray()[ nbindex ]; double r2 = 1.0 / rij2; double r6 = r2 * r2 * r2; @@ -546,7 +547,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) double f12 = LJ.A() * r12; // A/r^12 double f6 = LJ.B() * r6; // B/r^6 double e_vdw = f12 - f6; // (A/r^12)-(B/r^6) - Evdw += e_vdw; + Evdw += (e_vdw * vswitch); //mprintf("PVDW %8i%8i%20.6f%20.6f\n", ta0+1, ta1+1, e_vdw, r2); // LJ PME direct space correction double kr2 = lw_coeff_ * lw_coeff_ * rij2; @@ -554,7 +555,6 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) //double kr6 = kr2 * kr4; double expterm = exp(-kr2); double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; - vswitch = switch_fn(rij2, cut2_0, cut2_1); Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) * r6 * vswitch * Cij; } } From e387fd5c46a19494ce9ffcefbc1817959fa012de Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 11:40:50 -0500 Subject: [PATCH 21/36] DRR - Cpptraj: Add LJ switching argument --- src/Action_Energy.cpp | 8 +++++--- src/Action_Energy.h | 1 + src/Ewald_ParticleMesh.cpp | 6 ++++-- src/Ewald_ParticleMesh.h | 5 +++-- 4 files changed, 13 insertions(+), 7 deletions(-) diff --git a/src/Action_Energy.cpp b/src/Action_Energy.cpp index 805479a169..3f10507117 100644 --- a/src/Action_Energy.cpp +++ b/src/Action_Energy.cpp @@ -23,7 +23,7 @@ void Action_Energy::Help() const { "\t ewald [cut ] [dsumtol ] [rsumtol ]\n" "\t [ewcoeff ] [maxexp ] [skinnb ]\n" "\t [mlimits ,,] [erfcdx ]\n" - "\t pme [cut ] [dsumtol ] [order ]\n" + "\t pme [cut ] [dsumtol ] [order ] [ljswidth ]\n" "\t [ewcoeff ] [ljpme] [ewcoefflj] [skinnb ]\n" "\t [nfft ,,] [erfcdx ]\n" "\t } ]\n" @@ -144,7 +144,7 @@ Action::RetType Action_Energy::Init(ArgList& actionArgs, ActionInit& init, int d if (actionArgs.hasKey("ljpme")) lwcoeff_ = 0.4; lwcoeff_ = actionArgs.getKeyDouble("ewcoefflj", lwcoeff_); - // TODO arguments for switching + ljswidth_ = actionArgs.getKeyDouble("ljswidth", 0.0); skinnb_ = actionArgs.getKeyDouble("skinnb", 2.0); erfcDx_ = actionArgs.getKeyDouble("erfcdx", 0.0); npoints_ = actionArgs.getKeyInt("order", 6); @@ -306,6 +306,8 @@ Action::RetType Action_Energy::Init(ArgList& actionArgs, ActionInit& init, int d else mprintf("\tLennard-Jones PME Ewald coefficient will be set to elec. Ewald coefficient.\n"); } + if (ljswidth_ > 0.0) + mprintf("\tWidth of LJ switch region: %.4f Ang.\n", ljswidth_); } if (KEtype_ != KE_NONE) { if (KEtype_ == KE_AUTO) @@ -350,7 +352,7 @@ Action::RetType Action_Energy::Setup(ActionSetup& setup) { # ifdef LIBPME else if (elecType_ == PME) { if (((Ewald_ParticleMesh*)EW_)->Init(setup.CoordInfo().TrajBox(), cutoff_, dsumtol_, - ewcoeff_, lwcoeff_, skinnb_, erfcDx_, npoints_, + ewcoeff_, lwcoeff_, ljswidth_, skinnb_, erfcDx_, npoints_, debug_, mlimits_)) return Action::ERR; EW_->Setup( setup.Top(), Imask_ ); diff --git a/src/Action_Energy.h b/src/Action_Energy.h index d5520d68f5..979e637cba 100644 --- a/src/Action_Energy.h +++ b/src/Action_Energy.h @@ -46,6 +46,7 @@ class Action_Energy: public Action { double rsumtol_; ///< Regular Ewald reciprocal sum tolerance. double ewcoeff_; ///< Ewald coefficient. double lwcoeff_; ///< LJ Ewald coefficient. + double ljswidth_; ///< Size of LJ switch region double maxexp_; double skinnb_; ///< Size of non-bonded "skin" double erfcDx_; ///< Spacing for ERFC table (default 1/5000) diff --git a/src/Ewald_ParticleMesh.cpp b/src/Ewald_ParticleMesh.cpp index bf79aacec5..9f372ee949 100644 --- a/src/Ewald_ParticleMesh.cpp +++ b/src/Ewald_ParticleMesh.cpp @@ -88,10 +88,12 @@ int Ewald_ParticleMesh::DetermineNfft(int& nfft1, int& nfft2, int& nfft3, Box co /** Set up PME parameters. */ int Ewald_ParticleMesh::Init(Box const& boxIn, double cutoffIn, double dsumTolIn, - double ew_coeffIn, double lw_coeffIn,double skinnbIn, double erfcTableDxIn, + double ew_coeffIn, double lw_coeffIn, double switch_widthIn, + double skinnbIn, double erfcTableDxIn, int orderIn, int debugIn, const int* nfftIn) { - if (CheckInput(boxIn, debugIn, cutoffIn, dsumTolIn, ew_coeffIn, lw_coeffIn, 0.0, erfcTableDxIn, skinnbIn)) + if (CheckInput(boxIn, debugIn, cutoffIn, dsumTolIn, ew_coeffIn, lw_coeffIn, switch_widthIn, + erfcTableDxIn, skinnbIn)) return 1; if (nfftIn != 0) std::copy(nfftIn, nfftIn+3, nfft_); diff --git a/src/Ewald_ParticleMesh.h b/src/Ewald_ParticleMesh.h index a8864d9cb5..f7acbff3ec 100644 --- a/src/Ewald_ParticleMesh.h +++ b/src/Ewald_ParticleMesh.h @@ -6,8 +6,9 @@ class Ewald_ParticleMesh : public Ewald { public: Ewald_ParticleMesh(); - /// Box, cut, dsum tol, ew coeff, lj ew coeff, NB skin, erfc dx, order, debug, nfft - int Init(Box const&, double, double, double, double, double, double, int, int, const int*); + /// Box, cut, dsum tol, ew coeff, lj ew coeff, switch width, NB skin, erfc dx, order, dbg, nfft + int Init(Box const&, double, double, double, double, double, double, double, + int, int, const int*); // ----- Inherited --------------------------- int Setup(Topology const&, AtomMask const&); double CalcEnergy(Frame const&, AtomMask const&, double&); // TODO const? From e8741d417bcab0b6e49627c78ffac67a319d6566 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 13:26:54 -0500 Subject: [PATCH 22/36] DRR - Cpptraj: First incarnation of the catcrd command --- src/Exec_CatCrd.cpp | 92 +++++++++++++++++++++++++++++++++++++++++++++ src/Exec_CatCrd.h | 12 ++++++ 2 files changed, 104 insertions(+) create mode 100644 src/Exec_CatCrd.cpp create mode 100644 src/Exec_CatCrd.h diff --git a/src/Exec_CatCrd.cpp b/src/Exec_CatCrd.cpp new file mode 100644 index 0000000000..8dc62a46e5 --- /dev/null +++ b/src/Exec_CatCrd.cpp @@ -0,0 +1,92 @@ +#include "Exec_CatCrd.h" +#include "CpptrajStdio.h" + +// Exec_CatCrd::Help() +void Exec_CatCrd::Help() const +{ + mprintf("\t [ ...] name \n" + " Combine 2 or more COORDS data sets into a single one.\n"); +} + +// Exec_CatCrd::Execute() +Exec::RetType Exec_CatCrd::Execute(CpptrajState& State, ArgList& argIn) +{ + // Set up metadata with file name and output set name + std::string setname = argIn.GetStringKey("name"); + if (setname.empty()) { + mprinterr("Error: No output COORDS set name specified: 'name '\n"); + return CpptrajState::ERR; + } + Topology* parm = 0; + // Get COORDS data sets to concatenate. + typedef std::vector DCarray; + DCarray inputSets; + std::string dsarg = argIn.GetStringNext(); + while (!dsarg.empty()) { + DataSetList dsl = State.DSL().GetMultipleSets( dsarg ); + for (DataSetList::const_iterator ds = dsl.begin(); ds != dsl.end(); ++ds) + { + if ( (*ds)->Group() != DataSet::COORDINATES ) + mprintf("Warning: Set '%s' is not COORDS, skipping.\n", (*ds)->legend()); + else { + DataSet_Coords* coordsIn = (DataSet_Coords*)*ds; + if (parm == 0) + parm = coordsIn->TopPtr(); + else { + // Check that topology matches. For now just check # atoms. + if (parm->Natom() != coordsIn->Top().Natom()) { + mprinterr("Error: Set '%s' # atoms (%i) differs from first (%i)\n", + coordsIn->legend(), coordsIn->Top().Natom(), parm->Natom()); + return CpptrajState::ERR; + } + } + inputSets.push_back( coordsIn ); + } + } + } + if (inputSets.empty()) { + mprinterr("Error: No input COORDS sets.\n"); + return CpptrajState::ERR; + } + + // Check if output already present + DataSet_Coords* coordsOut = 0; + DataSet* ds = State.DSL().FindSetOfType( setname, DataSet::COORDS ); + if (ds == 0) { + // Create Set + MetaData md( setname ); + coordsOut = (DataSet_Coords*)State.DSL().AddSet(DataSet::COORDS, md); + if (coordsOut == 0) { + mprinterr("Error: Could not allocate COORDS data set.\n"); + return CpptrajState::ERR; + } + coordsOut->CoordsSetup( *parm, inputSets.front()->CoordsInfo() ); + mprintf("\tNew COORDS set '%s'\n", coordsOut->legend()); + } else { + // Check that set is actually coords. + if (ds->Type() != DataSet::COORDS) { + mprinterr("Error: Set %s present but is not of type COORDS.\n", ds->legend()); + return CpptrajState::ERR; + } + coordsOut = (DataSet_Coords*)ds; + // Check that topology matches. For now just check # atoms. + if (parm->Natom() != coordsOut->Top().Natom()) { + mprinterr("Error: # atoms %i does not match COORDS data set '%s' (%i)\n", + parm->Natom(), coordsOut->legend(), coordsOut->Top().Natom()); + return CpptrajState::ERR; + } + mprintf("\tAppending to COORDS data set '%s'\n", coordsOut->legend()); + } + + for (DCarray::const_iterator in = inputSets.begin(); in != inputSets.end(); ++in) + { + mprintf("\t'%s'\n", (*in)->legend()); + Frame frameIn = (*in)->AllocateFrame(); + for (unsigned int frm = 0; frm != (*in)->Size(); frm++) { + (*in)->GetFrame(frm, frameIn); + coordsOut->AddFrame( frameIn ); + } + } + + return CpptrajState::OK; +} diff --git a/src/Exec_CatCrd.h b/src/Exec_CatCrd.h new file mode 100644 index 0000000000..32f6861b60 --- /dev/null +++ b/src/Exec_CatCrd.h @@ -0,0 +1,12 @@ +#ifndef INC_EXEC_CATCRD_H +#define INC_EXEC_CATCRD_H +#include "Exec.h" +/// Concatenate two or more COORDS data sets +class Exec_CatCrd : public Exec { + public: + Exec_CatCrd() : Exec(COORDS) {} + void Help() const; + DispatchObject* Alloc() const { return (DispatchObject*)new Exec_CatCrd(); } + RetType Execute(CpptrajState&, ArgList&); +}; +#endif From 85c2ff11c42e0236ada0a9d1688c4eafd1c00d25 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 14:22:33 -0500 Subject: [PATCH 23/36] DRR - Cpptraj: Output switch width when active --- src/Ewald_ParticleMesh.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Ewald_ParticleMesh.cpp b/src/Ewald_ParticleMesh.cpp index 9f372ee949..0a3aaa90dd 100644 --- a/src/Ewald_ParticleMesh.cpp +++ b/src/Ewald_ParticleMesh.cpp @@ -109,6 +109,8 @@ int Ewald_ParticleMesh::Init(Box const& boxIn, double cutoffIn, double dsumTolIn cutoff_, dsumTol_, ew_coeff_, skinnbIn); if (lw_coeff_ > 0.0) mprintf("\t LJ Ewald coeff.= %g\n", lw_coeff_); + if (switch_width_ > 0.0) + mprintf("\t LJ switch width= %g\n", switch_width_); mprintf("\t Bspline order= %i\n", order_); mprintf("\t Erfc table dx= %g, size= %zu\n", erfcTableDx_, erfc_table_.size()/4); mprintf("\t "); From 82cb87e9da8277860ea0d19f58bb8f471e6c803c Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 14:22:54 -0500 Subject: [PATCH 24/36] DRR - Cpptraj: Enable catcrd --- src/Command.cpp | 2 ++ src/Exec_CatCrd.cpp | 1 + src/cpptrajdepend | 3 ++- src/cpptrajfiles | 1 + 4 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/Command.cpp b/src/Command.cpp index 7ecf86bfe4..77a492dd65 100644 --- a/src/Command.cpp +++ b/src/Command.cpp @@ -40,6 +40,7 @@ #include "Exec_PermuteDihedrals.h" #include "Exec_RotateDihedral.h" #include "Exec_SplitCoords.h" +#include "Exec_CatCrd.h" // ----- TRAJECTORY ------------------------------------------------------------ #include "Exec_Traj.h" // ----- TOPOLOGY -------------------------------------------------------------- @@ -233,6 +234,7 @@ void Command::Init() { // SYSTEM Command::AddCmd( new Exec_System(), Cmd::EXE, 6, "gnuplot", "head", "less", "ls", "pwd", "xmgrace" ); // COORDS + Command::AddCmd( new Exec_CatCrd(), Cmd::EXE, 1, "catcrd" ); Command::AddCmd( new Exec_CombineCoords(), Cmd::EXE, 1, "combinecrd" ); Command::AddCmd( new Exec_CrdAction(), Cmd::EXE, 1, "crdaction" ); Command::AddCmd( new Exec_CrdOut(), Cmd::EXE, 1, "crdout" ); diff --git a/src/Exec_CatCrd.cpp b/src/Exec_CatCrd.cpp index 8dc62a46e5..623fecd280 100644 --- a/src/Exec_CatCrd.cpp +++ b/src/Exec_CatCrd.cpp @@ -43,6 +43,7 @@ Exec::RetType Exec_CatCrd::Execute(CpptrajState& State, ArgList& argIn) inputSets.push_back( coordsIn ); } } + dsarg = argIn.GetStringNext(); } if (inputSets.empty()) { mprinterr("Error: No input COORDS sets.\n"); diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 178553e614..858c2754b4 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -154,7 +154,7 @@ Cluster_ReadInfo.o : Cluster_ReadInfo.cpp ArgList.h ArrayIterator.h AssociatedDa Cmd.o : Cmd.cpp Cmd.h DispatchObject.h CmdInput.o : CmdInput.cpp CmdInput.h StringRoutines.h CmdList.o : CmdList.cpp Cmd.h CmdList.h DispatchObject.h -Command.o : Command.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Action_Align.h Action_Angle.h Action_AreaPerMol.h Action_AtomMap.h Action_AtomicCorr.h Action_AtomicFluct.h Action_AutoImage.h Action_Average.h Action_Bounds.h Action_Box.h Action_Center.h Action_Channel.h Action_CheckChirality.h Action_CheckStructure.h Action_Closest.h Action_ClusterDihedral.h Action_Contacts.h Action_CreateCrd.h Action_CreateReservoir.h Action_DNAionTracker.h Action_DSSP.h Action_Density.h Action_Diffusion.h Action_Dihedral.h Action_Dipole.h Action_DistRmsd.h Action_Distance.h Action_Energy.h Action_Esander.h Action_FilterByData.h Action_FixAtomOrder.h Action_FixImagedBonds.h Action_GIST.h Action_Grid.h Action_GridFreeEnergy.h Action_HydrogenBond.h Action_Image.h Action_InfraredSpectrum.h Action_Jcoupling.h Action_LESsplit.h Action_LIE.h Action_LipidOrder.h Action_MakeStructure.h Action_Mask.h Action_Matrix.h Action_MinImage.h Action_Molsurf.h Action_MultiDihedral.h Action_MultiVector.h Action_NAstruct.h Action_NMRrst.h Action_NativeContacts.h Action_OrderParameter.h Action_Outtraj.h Action_PairDist.h Action_Pairwise.h Action_Principal.h Action_Projection.h Action_Pucker.h Action_Radgyr.h Action_Radial.h Action_RandomizeIons.h Action_Remap.h Action_ReplicateCell.h Action_Rmsd.h Action_Rotate.h Action_RunningAvg.h Action_STFC_Diffusion.h Action_Scale.h Action_SetVelocity.h Action_Spam.h Action_Strip.h Action_Surf.h Action_SymmetricRmsd.h Action_Temperature.h Action_Translate.h Action_Unstrip.h Action_Unwrap.h Action_Vector.h Action_VelocityAutoCorr.h Action_Volmap.h Action_Volume.h Action_Watershell.h Analysis.h AnalysisList.h AnalysisState.h Analysis_AmdBias.h Analysis_AutoCorr.h Analysis_Average.h Analysis_Clustering.h Analysis_ConstantPHStats.h Analysis_Corr.h Analysis_CrankShaft.h Analysis_CrdFluct.h Analysis_CrossCorr.h Analysis_CurveFit.h Analysis_Divergence.h Analysis_FFT.h Analysis_HausdorffDistance.h Analysis_Hist.h Analysis_IRED.h Analysis_Integrate.h Analysis_KDE.h Analysis_Lifetime.h Analysis_LowestCurve.h Analysis_Matrix.h Analysis_MeltCurve.h Analysis_Modes.h Analysis_MultiHist.h Analysis_Multicurve.h Analysis_Overlap.h Analysis_PhiPsi.h Analysis_Regression.h Analysis_RemLog.h Analysis_Rms2d.h Analysis_RmsAvgCorr.h Analysis_Rotdif.h Analysis_RunningAvg.h Analysis_Spline.h Analysis_State.h Analysis_Statistics.h Analysis_TI.h Analysis_Timecorr.h Analysis_VectorMath.h Analysis_Wavelet.h ArgList.h Array1D.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AxisType.h BaseIOtype.h Box.h BufferedLine.h CharMask.h ClusterDist.h ClusterList.h ClusterMap.h ClusterNode.h ClusterSieve.h Cmd.h CmdInput.h CmdList.h Command.h ComplexArray.h Constraints.h Control.h CoordinateInfo.h Corr.h Cph.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Cmatrix.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_GridFlt.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_RemLog.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_pH.h DataSet_string.h Deprecated.h DihedralSearch.h Dimension.h DispatchObject.h DistRoutines.h Energy.h Energy_Sander.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h Ewald.h Exec.h Exec_Analyze.h Exec_Calc.h Exec_Change.h Exec_ClusterMap.h Exec_CombineCoords.h Exec_Commands.h Exec_CompareTop.h Exec_CrdAction.h Exec_CrdOut.h Exec_CreateSet.h Exec_DataFile.h Exec_DataFilter.h Exec_DataSetCmd.h Exec_GenerateAmberRst.h Exec_Help.h Exec_LoadCrd.h Exec_LoadTraj.h Exec_ParallelAnalysis.h Exec_ParmBox.h Exec_ParmSolvent.h Exec_ParmStrip.h Exec_ParmWrite.h Exec_PermuteDihedrals.h Exec_Precision.h Exec_PrintData.h Exec_ReadData.h Exec_ReadEnsembleData.h Exec_ReadInput.h Exec_RotateDihedral.h Exec_RunAnalysis.h Exec_ScaleDihedralK.h Exec_SequenceAlign.h Exec_SortEnsembleData.h Exec_SplitCoords.h Exec_System.h Exec_Top.h Exec_Traj.h Exec_UpdateParameters.h Exec_ViewRst.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h Grid.h GridAction.h GridBin.h HistBin.h Hungarian.h ImageTypes.h ImagedAction.h InputTrajCommon.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h NetcdfFile.h OnlineVarT.h OutputTrajCommon.h PDBfile.h PairList.h Parallel.h ParameterHolders.h ParameterTypes.h PubFFT.h RPNcalc.h Random.h Range.h ReferenceAction.h ReferenceFrame.h RemdReservoirNC.h ReplicaDimArray.h ReplicaInfo.h Residue.h Spline.h StructureCheck.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h VariableArray.h Vec3.h molsurf.h +Command.o : Command.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Action_Align.h Action_Angle.h Action_AreaPerMol.h Action_AtomMap.h Action_AtomicCorr.h Action_AtomicFluct.h Action_AutoImage.h Action_Average.h Action_Bounds.h Action_Box.h Action_Center.h Action_Channel.h Action_CheckChirality.h Action_CheckStructure.h Action_Closest.h Action_ClusterDihedral.h Action_Contacts.h Action_CreateCrd.h Action_CreateReservoir.h Action_DNAionTracker.h Action_DSSP.h Action_Density.h Action_Diffusion.h Action_Dihedral.h Action_Dipole.h Action_DistRmsd.h Action_Distance.h Action_Energy.h Action_Esander.h Action_FilterByData.h Action_FixAtomOrder.h Action_FixImagedBonds.h Action_GIST.h Action_Grid.h Action_GridFreeEnergy.h Action_HydrogenBond.h Action_Image.h Action_InfraredSpectrum.h Action_Jcoupling.h Action_LESsplit.h Action_LIE.h Action_LipidOrder.h Action_MakeStructure.h Action_Mask.h Action_Matrix.h Action_MinImage.h Action_Molsurf.h Action_MultiDihedral.h Action_MultiVector.h Action_NAstruct.h Action_NMRrst.h Action_NativeContacts.h Action_OrderParameter.h Action_Outtraj.h Action_PairDist.h Action_Pairwise.h Action_Principal.h Action_Projection.h Action_Pucker.h Action_Radgyr.h Action_Radial.h Action_RandomizeIons.h Action_Remap.h Action_ReplicateCell.h Action_Rmsd.h Action_Rotate.h Action_RunningAvg.h Action_STFC_Diffusion.h Action_Scale.h Action_SetVelocity.h Action_Spam.h Action_Strip.h Action_Surf.h Action_SymmetricRmsd.h Action_Temperature.h Action_Translate.h Action_Unstrip.h Action_Unwrap.h Action_Vector.h Action_VelocityAutoCorr.h Action_Volmap.h Action_Volume.h Action_Watershell.h Analysis.h AnalysisList.h AnalysisState.h Analysis_AmdBias.h Analysis_AutoCorr.h Analysis_Average.h Analysis_Clustering.h Analysis_ConstantPHStats.h Analysis_Corr.h Analysis_CrankShaft.h Analysis_CrdFluct.h Analysis_CrossCorr.h Analysis_CurveFit.h Analysis_Divergence.h Analysis_FFT.h Analysis_HausdorffDistance.h Analysis_Hist.h Analysis_IRED.h Analysis_Integrate.h Analysis_KDE.h Analysis_Lifetime.h Analysis_LowestCurve.h Analysis_Matrix.h Analysis_MeltCurve.h Analysis_Modes.h Analysis_MultiHist.h Analysis_Multicurve.h Analysis_Overlap.h Analysis_PhiPsi.h Analysis_Regression.h Analysis_RemLog.h Analysis_Rms2d.h Analysis_RmsAvgCorr.h Analysis_Rotdif.h Analysis_RunningAvg.h Analysis_Spline.h Analysis_State.h Analysis_Statistics.h Analysis_TI.h Analysis_Timecorr.h Analysis_VectorMath.h Analysis_Wavelet.h ArgList.h Array1D.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AxisType.h BaseIOtype.h Box.h BufferedLine.h CharMask.h ClusterDist.h ClusterList.h ClusterMap.h ClusterNode.h ClusterSieve.h Cmd.h CmdInput.h CmdList.h Command.h ComplexArray.h Constraints.h Control.h CoordinateInfo.h Corr.h Cph.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Cmatrix.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_GridFlt.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_RemLog.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_pH.h DataSet_string.h Deprecated.h DihedralSearch.h Dimension.h DispatchObject.h DistRoutines.h Energy.h Energy_Sander.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h Ewald.h Exec.h Exec_Analyze.h Exec_Calc.h Exec_CatCrd.h Exec_Change.h Exec_ClusterMap.h Exec_CombineCoords.h Exec_Commands.h Exec_CompareTop.h Exec_CrdAction.h Exec_CrdOut.h Exec_CreateSet.h Exec_DataFile.h Exec_DataFilter.h Exec_DataSetCmd.h Exec_GenerateAmberRst.h Exec_Help.h Exec_LoadCrd.h Exec_LoadTraj.h Exec_ParallelAnalysis.h Exec_ParmBox.h Exec_ParmSolvent.h Exec_ParmStrip.h Exec_ParmWrite.h Exec_PermuteDihedrals.h Exec_Precision.h Exec_PrintData.h Exec_ReadData.h Exec_ReadEnsembleData.h Exec_ReadInput.h Exec_RotateDihedral.h Exec_RunAnalysis.h Exec_ScaleDihedralK.h Exec_SequenceAlign.h Exec_SortEnsembleData.h Exec_SplitCoords.h Exec_System.h Exec_Top.h Exec_Traj.h Exec_UpdateParameters.h Exec_ViewRst.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h Grid.h GridAction.h GridBin.h HistBin.h Hungarian.h ImageTypes.h ImagedAction.h InputTrajCommon.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h NetcdfFile.h OnlineVarT.h OutputTrajCommon.h PDBfile.h PairList.h Parallel.h ParameterHolders.h ParameterTypes.h PubFFT.h RPNcalc.h Random.h Range.h ReferenceAction.h ReferenceFrame.h RemdReservoirNC.h ReplicaDimArray.h ReplicaInfo.h Residue.h Spline.h StructureCheck.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h VariableArray.h Vec3.h molsurf.h ComplexArray.o : ComplexArray.cpp ArrayIterator.h ComplexArray.h Constraints.o : Constraints.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h Box.h CharMask.h Constants.h Constraints.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Topology.h Vec3.h Control.o : Control.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Control.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h VariableArray.h Vec3.h @@ -235,6 +235,7 @@ Ewald_ParticleMesh.o : Ewald_ParticleMesh.cpp Atom.h AtomExtra.h AtomMask.h Box. Ewald_Regular.o : Ewald_Regular.cpp Atom.h AtomExtra.h AtomMask.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajStdio.h Ewald.h Ewald_Regular.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h PairList.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h Timer.h Topology.h Vec3.h Exec_Analyze.o : Exec_Analyze.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Cmd.h CmdList.h Command.h Control.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h Exec.h Exec_Analyze.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h VariableArray.h Vec3.h Exec_Calc.o : Exec_Calc.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h Exec.h Exec_Calc.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h RPNcalc.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h Vec3.h +Exec_CatCrd.o : Exec_CatCrd.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h Exec.h Exec_CatCrd.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h Vec3.h Exec_Change.o : Exec_Change.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h Exec.h Exec_Change.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterHolders.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h Vec3.h Exec_ClusterMap.o : Exec_ClusterMap.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h ClusterMap.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_2D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_MatrixFlt.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h Exec.h Exec_ClusterMap.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h Vec3.h Exec_CombineCoords.o : Exec_CombineCoords.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h Exec.h Exec_CombineCoords.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h Vec3.h diff --git a/src/cpptrajfiles b/src/cpptrajfiles index f9395fd816..633a1da86f 100644 --- a/src/cpptrajfiles +++ b/src/cpptrajfiles @@ -232,6 +232,7 @@ COMMON_SOURCES=ActionFrameCounter.cpp \ Ewald_Regular.cpp \ Exec_Analyze.cpp \ Exec_Calc.cpp \ + Exec_CatCrd.cpp \ Exec_Change.cpp \ Exec_ClusterMap.cpp \ Exec_CombineCoords.cpp \ From 02217e0a38809cb49848c85bcc42c5107c83a996 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 15:10:32 -0500 Subject: [PATCH 25/36] DRR - Cpptraj: Add test for catcrd command. --- test/Makefile | 3 +++ test/Test_CatCrd/RunTest.sh | 28 ++++++++++++++++++++++++++++ 2 files changed, 31 insertions(+) create mode 100755 test/Test_CatCrd/RunTest.sh diff --git a/test/Makefile b/test/Makefile index 88060d3656..598355e1f4 100644 --- a/test/Makefile +++ b/test/Makefile @@ -444,6 +444,9 @@ test.xyzfmt: test.hausdorff: @-cd Test_Hausdorff && ./RunTest.sh $(OPT) +test.catcrd: + @-cd Test_CatCrd && ./RunTest.sh $(OPT) + # Every test target should go here. COMPLETETESTS=test.general \ test.strip \ diff --git a/test/Test_CatCrd/RunTest.sh b/test/Test_CatCrd/RunTest.sh new file mode 100755 index 0000000000..02a50845b5 --- /dev/null +++ b/test/Test_CatCrd/RunTest.sh @@ -0,0 +1,28 @@ +#!/bin/bash + +. ../MasterTest.sh + +CleanFiles cat.in cat.crd.save cat.crd + +TESTNAME='Concatenate COORDS data set test' +Requires maxthreads 1 + +INPUT='-i cat.in' + +cat > cat.in < Date: Thu, 14 Feb 2019 15:12:26 -0500 Subject: [PATCH 26/36] DRR - Cpptraj: Enable catcrd test --- test/Makefile | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/Makefile b/test/Makefile index 598355e1f4..8406e62b0c 100644 --- a/test/Makefile +++ b/test/Makefile @@ -589,7 +589,8 @@ COMPLETETESTS=test.general \ test.cphstats \ test.cmdline \ test.xyzfmt \ - test.hausdorff + test.hausdorff \ + test.catcrd test.all: $(MAKE) test.complete summary From e7445dcdec509ba81f78c680837b6f00cc1082c7 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 15:12:53 -0500 Subject: [PATCH 27/36] DRR - Cpptraj: Revision bump; LJ PME support and catcrd command. --- src/Version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Version.h b/src/Version.h index c711599bcd..25b535bb24 100644 --- a/src/Version.h +++ b/src/Version.h @@ -20,5 +20,5 @@ * Whenever a number that precedes is incremented, all subsequent * numbers should be reset to 0. */ -#define CPPTRAJ_INTERNAL_VERSION "V4.11.1" +#define CPPTRAJ_INTERNAL_VERSION "V4.11.2" #endif From 0dafd9d2ce0021ae2bb65a0f6f2092408c9b3219 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 14 Feb 2019 15:19:46 -0500 Subject: [PATCH 28/36] DRR - Cpptraj: Add LJ with switch test --- test/Test_LJ_PME/RunTest.sh | 17 ++++++++++++++++- test/Test_LJ_PME/switch.dat.save | 2 ++ 2 files changed, 18 insertions(+), 1 deletion(-) create mode 100644 test/Test_LJ_PME/switch.dat.save diff --git a/test/Test_LJ_PME/RunTest.sh b/test/Test_LJ_PME/RunTest.sh index 41eac4d613..da749c08d8 100755 --- a/test/Test_LJ_PME/RunTest.sh +++ b/test/Test_LJ_PME/RunTest.sh @@ -2,7 +2,7 @@ . ../MasterTest.sh -CleanFiles lj.in ene.dat ene.dat.? +CleanFiles lj.in ene.dat ene.dat.? switch.dat INPUT='-i lj.in' TESTNAME='LJ PME tests.' @@ -39,6 +39,21 @@ EOF ((i++)) done +# Test with switching +cat > lj.in < Date: Fri, 15 Feb 2019 09:34:39 -0500 Subject: [PATCH 29/36] DRR - Cpptraj: Attempt to simplify code maintenance while maintaining performance by putting common energy kernels in include files. --- src/EnergyKernel_Adjust.h | 11 ++++ src/EnergyKernel_Nonbond.h | 41 ++++++++++++++ src/Ewald.cpp | 107 ++----------------------------------- src/cpptrajdepend | 2 +- 4 files changed, 57 insertions(+), 104 deletions(-) create mode 100644 src/EnergyKernel_Adjust.h create mode 100644 src/EnergyKernel_Nonbond.h diff --git a/src/EnergyKernel_Adjust.h b/src/EnergyKernel_Adjust.h new file mode 100644 index 0000000000..c5d7423138 --- /dev/null +++ b/src/EnergyKernel_Adjust.h @@ -0,0 +1,11 @@ + e_adjust += Adjust(q0, q1, sqrt(rij2)); + // LJ PME direct space correction + // NOTE: Assuming excluded pair is within cutoff + double kr2 = lw_coeff_ * lw_coeff_ * rij2; + double kr4 = kr2 * kr2; + //double kr6 = kr2 * kr4; + double expterm = exp(-kr2); + double r4 = rij2 * rij2; + double r6 = rij2 * r4; + double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; + Eljpme_correction_excl += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) / r6 * Cij; diff --git a/src/EnergyKernel_Nonbond.h b/src/EnergyKernel_Nonbond.h new file mode 100644 index 0000000000..4362dba5cb --- /dev/null +++ b/src/EnergyKernel_Nonbond.h @@ -0,0 +1,41 @@ + double rij = sqrt( rij2 ); + double qiqj = q0 * q1; +# ifndef _OPENMP + t_erfc_.Start(); +# endif + //double erfc = erfc_func(ew_coeff_ * rij); + double erfc = ERFC(ew_coeff_ * rij); +# ifndef _OPENMP + t_erfc_.Stop(); +# endif + double e_elec = qiqj * erfc / rij; + Eelec += e_elec; + //mprintf("EELEC %4i%4i%12.5f%12.5f%12.5f%3.0f%3.0f%3.0f\n", + //int ta0, ta1; + //if (it0->Idx() < it1->Idx()) { + // ta0=it0->Idx(); ta1=it1->Idx(); + //} else { + // ta1=it0->Idx(); ta0=it1->Idx(); + //} + //mprintf("PELEC %6i%6i%12.5f%12.5f%12.5f\n", ta0, ta1, rij, erfc, e_elec); + int nbindex = NB_->GetLJindex(TypeIndices_[it0->Idx()], + TypeIndices_[it1->Idx()]); + if (nbindex > -1) { + vswitch = switch_fn(rij2, cut2_0, cut2_1); + NonbondType const& LJ = NB_->NBarray()[ nbindex ]; + double r2 = 1.0 / rij2; + double r6 = r2 * r2 * r2; + double r12 = r6 * r6; + double f12 = LJ.A() * r12; // A/r^12 + double f6 = LJ.B() * r6; // B/r^6 + double e_vdw = f12 - f6; // (A/r^12)-(B/r^6) + Evdw += (e_vdw * vswitch); + //mprintf("PVDW %8i%8i%20.6f%20.6f\n", ta0+1, ta1+1, e_vdw, r2); + // LJ PME direct space correction + double kr2 = lw_coeff_ * lw_coeff_ * rij2; + double kr4 = kr2 * kr2; + //double kr6 = kr2 * kr4; + double expterm = exp(-kr2); + double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; + Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) * r6 * vswitch * Cij; + } diff --git a/src/Ewald.cpp b/src/Ewald.cpp index 8cb8ccb9ef..af25477de0 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -433,59 +433,10 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) if (excluded.find( it1->Idx() ) == excluded.end()) { if ( rij2 < cut2 ) { - double rij = sqrt( rij2 ); - double qiqj = q0 * q1; -# ifndef _OPENMP - t_erfc_.Start(); -# endif - //double erfc = erfc_func(ew_coeff_ * rij); - double erfc = ERFC(ew_coeff_ * rij); -# ifndef _OPENMP - t_erfc_.Stop(); -# endif - double e_elec = qiqj * erfc / rij; - Eelec += e_elec; - //int ta0, ta1; - //if (it0->Idx() < it1->Idx()) { - // ta0=it0->Idx(); ta1=it1->Idx(); - //} else { - // ta1=it0->Idx(); ta0=it1->Idx(); - //} - //mprintf("PELEC %6i%6i%12.5f%12.5f%12.5f\n", ta0, ta1, rij, erfc, e_elec); - int nbindex = NB_->GetLJindex(TypeIndices_[it0->Idx()], - TypeIndices_[it1->Idx()]); - if (nbindex > -1) { - vswitch = switch_fn(rij2, cut2_0, cut2_1); - NonbondType const& LJ = NB_->NBarray()[ nbindex ]; - double r2 = 1.0 / rij2; - double r6 = r2 * r2 * r2; - double r12 = r6 * r6; - double f12 = LJ.A() * r12; // A/r^12 - double f6 = LJ.B() * r6; // B/r^6 - double e_vdw = f12 - f6; // (A/r^12)-(B/r^6) - Evdw += (e_vdw * vswitch); - //mprintf("PVDW %8i%8i%20.6f%20.6f\n", ta0+1, ta1+1, e_vdw, r2); - // LJ PME direct space correction - double kr2 = lw_coeff_ * lw_coeff_ * rij2; - double kr4 = kr2 * kr2; - //double kr6 = kr2 * kr4; - double expterm = exp(-kr2); - double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; - Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) * r6 * vswitch * Cij; - } +# include "EnergyKernel_Nonbond.h" } } else { - e_adjust += Adjust(q0, q1, sqrt(rij2)); - // LJ PME direct space correction - // NOTE: Assuming excluded pair is within cutoff - double kr2 = lw_coeff_ * lw_coeff_ * rij2; - double kr4 = kr2 * kr2; - //double kr6 = kr2 * kr4; - double expterm = exp(-kr2); - double r4 = rij2 * rij2; - double r6 = rij2 * r4; - double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; - Eljpme_correction_excl += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) / r6 * Cij; +# include "EnergyKernel_Adjust.h" } } // END loop over other atoms in thisCell // Loop over all neighbor cells @@ -516,60 +467,10 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) { //mprintf("\t\t\tdist= %f\n", sqrt(rij2)); if ( rij2 < cut2 ) { - double rij = sqrt( rij2 ); - double qiqj = q0 * q1; -# ifndef _OPENMP - t_erfc_.Start(); -# endif - //double erfc = erfc_func(ew_coeff_ * rij); - double erfc = ERFC(ew_coeff_ * rij); -# ifndef _OPENMP - t_erfc_.Stop(); -# endif - double e_elec = qiqj * erfc / rij; - Eelec += e_elec; - //mprintf("EELEC %4i%4i%12.5f%12.5f%12.5f%3.0f%3.0f%3.0f\n", - //int ta0, ta1; - //if (it0->Idx() < it1->Idx()) { - // ta0=it0->Idx(); ta1=it1->Idx(); - //} else { - // ta1=it0->Idx(); ta0=it1->Idx(); - //} - //mprintf("PELEC %6i%6i%12.5f%12.5f%12.5f\n", ta0, ta1, rij, erfc, e_elec); - int nbindex = NB_->GetLJindex(TypeIndices_[it0->Idx()], - TypeIndices_[it1->Idx()]); - if (nbindex > -1) { - vswitch = switch_fn(rij2, cut2_0, cut2_1); - NonbondType const& LJ = NB_->NBarray()[ nbindex ]; - double r2 = 1.0 / rij2; - double r6 = r2 * r2 * r2; - double r12 = r6 * r6; - double f12 = LJ.A() * r12; // A/r^12 - double f6 = LJ.B() * r6; // B/r^6 - double e_vdw = f12 - f6; // (A/r^12)-(B/r^6) - Evdw += (e_vdw * vswitch); - //mprintf("PVDW %8i%8i%20.6f%20.6f\n", ta0+1, ta1+1, e_vdw, r2); - // LJ PME direct space correction - double kr2 = lw_coeff_ * lw_coeff_ * rij2; - double kr4 = kr2 * kr2; - //double kr6 = kr2 * kr4; - double expterm = exp(-kr2); - double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; - Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) * r6 * vswitch * Cij; - } +# include "EnergyKernel_Nonbond.h" } } else { - e_adjust += Adjust(q0, q1, sqrt(rij2)); - // LJ PME direct space correction - // NOTE: Assuming excluded pair is within cutoff - double kr2 = lw_coeff_ * lw_coeff_ * rij2; - double kr4 = kr2 * kr2; - //double kr6 = kr2 * kr4; - double expterm = exp(-kr2); - double r4 = rij2 * rij2; - double r6 = rij2 * r4; - double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; - Eljpme_correction_excl += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) / r6 * Cij; +# include "EnergyKernel_Adjust.h" } } // END loop over neighbor cell atoms } // END Loop over neighbor cells diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 858c2754b4..61e567cb76 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -230,7 +230,7 @@ EnsembleOut.o : EnsembleOut.cpp ActionFrameCounter.h ArgList.h Atom.h AtomExtra. EnsembleOutList.o : EnsembleOutList.cpp ActionFrameCounter.h ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h EnsembleOut.h EnsembleOutList.h EnsembleOut_Multi.h EnsembleOut_Single.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Topology.h TrajectoryFile.h TrajectoryIO.h Vec3.h EnsembleOut_Multi.o : EnsembleOut_Multi.cpp ActionFrameCounter.h ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h EnsembleOut.h EnsembleOut_Multi.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h Topology.h TrajectoryFile.h TrajectoryIO.h Vec3.h EnsembleOut_Single.o : EnsembleOut_Single.cpp ActionFrameCounter.h ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h EnsembleOut.h EnsembleOut_Single.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Topology.h TrajectoryFile.h TrajectoryIO.h Vec3.h -Ewald.o : Ewald.cpp Atom.h AtomExtra.h AtomMask.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajStdio.h Ewald.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h PairList.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h SymbolExporting.h Timer.h Topology.h Vec3.h +Ewald.o : Ewald.cpp Atom.h AtomExtra.h AtomMask.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajStdio.h EnergyKernel_Adjust.h EnergyKernel_Nonbond.h Ewald.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h PairList.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h SymbolExporting.h Timer.h Topology.h Vec3.h Ewald_ParticleMesh.o : Ewald_ParticleMesh.cpp Atom.h AtomExtra.h AtomMask.h Box.h CharMask.h CoordinateInfo.h CpptrajStdio.h Ewald.h Ewald_ParticleMesh.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h PairList.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Timer.h Topology.h Vec3.h helpme_standalone.h Ewald_Regular.o : Ewald_Regular.cpp Atom.h AtomExtra.h AtomMask.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajStdio.h Ewald.h Ewald_Regular.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h PairList.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h Timer.h Topology.h Vec3.h Exec_Analyze.o : Exec_Analyze.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Cmd.h CmdList.h Command.h Control.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h Exec.h Exec_Analyze.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h VariableArray.h Vec3.h From fd52ce46f7d386d0d559534592a3625533a8bddd Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 15 Feb 2019 10:03:36 -0500 Subject: [PATCH 30/36] DRR - Cpptraj: Put pair list loop as an include file; this somewhat improves the speed of the direct space calc for cases where e.g. the lj pme part is not needed. --- src/EnergyKernel_Adjust.h | 2 + src/EnergyKernel_Nonbond.h | 4 +- src/Ewald.cpp | 71 ++++++++++++++++++++++++++++++++++ src/Ewald.h | 5 ++- src/PairListLoop.h | 79 ++++++++++++++++++++++++++++++++++++++ src/cpptrajdepend | 2 +- 6 files changed, 160 insertions(+), 3 deletions(-) create mode 100644 src/PairListLoop.h diff --git a/src/EnergyKernel_Adjust.h b/src/EnergyKernel_Adjust.h index c5d7423138..8ea81e953a 100644 --- a/src/EnergyKernel_Adjust.h +++ b/src/EnergyKernel_Adjust.h @@ -1,4 +1,5 @@ e_adjust += Adjust(q0, q1, sqrt(rij2)); +# ifdef CPPTRAJ_EKERNEL_LJPME // LJ PME direct space correction // NOTE: Assuming excluded pair is within cutoff double kr2 = lw_coeff_ * lw_coeff_ * rij2; @@ -9,3 +10,4 @@ double r6 = rij2 * r4; double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; Eljpme_correction_excl += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) / r6 * Cij; +# endif diff --git a/src/EnergyKernel_Nonbond.h b/src/EnergyKernel_Nonbond.h index 4362dba5cb..b44a426599 100644 --- a/src/EnergyKernel_Nonbond.h +++ b/src/EnergyKernel_Nonbond.h @@ -21,7 +21,7 @@ int nbindex = NB_->GetLJindex(TypeIndices_[it0->Idx()], TypeIndices_[it1->Idx()]); if (nbindex > -1) { - vswitch = switch_fn(rij2, cut2_0, cut2_1); + double vswitch = switch_fn(rij2, cut2_0, cut2_1); NonbondType const& LJ = NB_->NBarray()[ nbindex ]; double r2 = 1.0 / rij2; double r6 = r2 * r2 * r2; @@ -31,6 +31,7 @@ double e_vdw = f12 - f6; // (A/r^12)-(B/r^6) Evdw += (e_vdw * vswitch); //mprintf("PVDW %8i%8i%20.6f%20.6f\n", ta0+1, ta1+1, e_vdw, r2); +# ifdef CPPTRAJ_EKERNEL_LJPME // LJ PME direct space correction double kr2 = lw_coeff_ * lw_coeff_ * rij2; double kr4 = kr2 * kr2; @@ -38,4 +39,5 @@ double expterm = exp(-kr2); double Cij = Cparam_[it0->Idx()] * Cparam_[it1->Idx()]; Eljpme_correction += (1.0 - (1.0 + kr2 + kr4/2.0)*expterm) * r6 * vswitch * Cij; +# endif } diff --git a/src/Ewald.cpp b/src/Ewald.cpp index af25477de0..04b43f849c 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -370,6 +370,71 @@ static inline double switch_fn(double rij2, double cut2_0, double cut2_1) } } +double Ewald::Direct_VDW_LongRangeCorrection(PairList const& PL, double& e_adjust_out, double& evdw_out) +{ + t_direct_.Start(); + double cut2 = cutoff_ * cutoff_; + double Eelec = 0.0; + double e_adjust = 0.0; + double Evdw = 0.0; + int cidx; + double cut0 = cutoff_ - switch_width_; + double cut2_0 = cut0 * cut0; + double cut2_1 = cut2; +# ifdef _OPENMP +# pragma omp parallel private(cidx) reduction(+: Eelec, Evdw, e_adjust) + { +# pragma omp for +# endif +# include "PairListLoop.h" +# ifdef _OPENMP + } // END pragma omp parallel +# endif + t_direct_.Stop(); + e_adjust_out = e_adjust; +# ifdef DEBUG_PAIRLIST + mprintf("DEBUG: LJ vdw = %16.8f\n", Evdw); +# endif + evdw_out = Evdw; + return Eelec; +} + +double Ewald::Direct_VDW_LJPME(PairList const& PL, double& e_adjust_out, double& evdw_out) +{ + t_direct_.Start(); + double cut2 = cutoff_ * cutoff_; + double Eelec = 0.0; + double e_adjust = 0.0; + double Evdw = 0.0; + double Eljpme_correction = 0.0; + double Eljpme_correction_excl = 0.0; + int cidx; + double cut0 = cutoff_ - switch_width_; + double cut2_0 = cut0 * cut0; + double cut2_1 = cut2; +# define CPPTRAJ_EKERNEL_LJPME +# ifdef _OPENMP +# pragma omp parallel private(cidx) reduction(+: Eelec, Evdw, e_adjust, Eljpme_correction,Eljpme_correction_excl) + { +# pragma omp for +# endif +# include "PairListLoop.h" +# ifdef _OPENMP + } // END pragma omp parallel +# endif +# undef CPPTRAJ_EKERNEL_LJPME + t_direct_.Stop(); + e_adjust_out = e_adjust; +# ifdef DEBUG_PAIRLIST + mprintf("DEBUG: LJ vdw = %16.8f\n", Evdw); + mprintf("DEBUG: LJ vdw PME correction = %16.8f\n", Eljpme_correction); + mprintf("DEBUG: LJ vdw PME correction (excluded) = %16.8f\n", Eljpme_correction_excl); +# endif + evdw_out = Evdw + Eljpme_correction + Eljpme_correction_excl; + return Eelec; +} + + // Ewald::Direct() /** Calculate direct space energy. This is the faster version that uses * a pair list. Also calculate the energy adjustment for excluded @@ -380,6 +445,11 @@ static inline double switch_fn(double rij2, double cut2_0, double cut2_1) */ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) { + if (lw_coeff_ > 0.0) + return Direct_VDW_LJPME(PL, e_adjust_out, evdw_out); + else + return Direct_VDW_LongRangeCorrection(PL, e_adjust_out, evdw_out); +/* t_direct_.Start(); double cut2 = cutoff_ * cutoff_; double Eelec = 0.0; @@ -489,6 +559,7 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) # endif evdw_out = Evdw + Eljpme_correction + Eljpme_correction_excl; return Eelec; +*/ } /** Determine VDW long range correction prefactor. */ diff --git a/src/Ewald.h b/src/Ewald.h index ca6575ded5..fac2b33619 100644 --- a/src/Ewald.h +++ b/src/Ewald.h @@ -91,12 +91,15 @@ class Ewald { double one_over_Dx_; ///< One over erfcTableDx_. double Vdw_Recip_term_; ///< VDW recip correction term from # types and B parameters int debug_; - Timer t_total_; + Timer t_total_; // TODO make timing external Timer t_self_; Timer t_recip_; Timer t_trig_tables_; Timer t_direct_; Timer t_erfc_; Timer t_adjust_; + private: + double Direct_VDW_LongRangeCorrection(PairList const&, double&, double&); + double Direct_VDW_LJPME(PairList const&, double&, double&); }; #endif diff --git a/src/PairListLoop.h b/src/PairListLoop.h new file mode 100644 index 0000000000..0dad69c27a --- /dev/null +++ b/src/PairListLoop.h @@ -0,0 +1,79 @@ + for (cidx = 0; cidx < PL.NGridMax(); cidx++) + { + PairList::CellType const& thisCell = PL.Cell( cidx ); + if (thisCell.NatomsInGrid() > 0) + { + // cellList contains this cell index and all neighbors. + PairList::Iarray const& cellList = thisCell.CellList(); + // transList contains index to translation for the neighbor. + PairList::Iarray const& transList = thisCell.TransList(); + // Loop over all atoms of thisCell. + for (PairList::CellType::const_iterator it0 = thisCell.begin(); + it0 != thisCell.end(); ++it0) + { + Vec3 const& xyz0 = it0->ImageCoords(); + double q0 = Charge_[it0->Idx()]; +# ifdef DEBUG_PAIRLIST + mprintf("DBG: Cell %6i (%6i atoms):\n", cidx+1, thisCell.NatomsInGrid()); +# endif + // Exclusion list for this atom + Iset const& excluded = Excluded_[it0->Idx()]; + // Calc interaction of atom to all other atoms in thisCell. + for (PairList::CellType::const_iterator it1 = it0 + 1; + it1 != thisCell.end(); ++it1) + { + Vec3 const& xyz1 = it1->ImageCoords(); + double q1 = Charge_[it1->Idx()]; + Vec3 dxyz = xyz1 - xyz0; + double rij2 = dxyz.Magnitude2(); +# ifdef DEBUG_PAIRLIST + mprintf("\tAtom %6i to atom %6i (%f)\n", it0->Idx()+1, it1->Idx()+1, sqrt(rij2)); +# endif + // If atom excluded, calc adjustment, otherwise calc elec. energy. + if (excluded.find( it1->Idx() ) == excluded.end()) + { + if ( rij2 < cut2 ) { +# include "EnergyKernel_Nonbond.h" + } + } else { +# include "EnergyKernel_Adjust.h" + } + } // END loop over other atoms in thisCell + // Loop over all neighbor cells + for (unsigned int nidx = 1; nidx != cellList.size(); nidx++) + { + PairList::CellType const& nbrCell = PL.Cell( cellList[nidx] ); +# ifdef DEBUG_PAIRLIST + if (nbrCell.NatomsInGrid()>0) mprintf("\tto neighbor cell %6i\n", cellList[nidx]+1); +# endif + // Translate vector for neighbor cell + Vec3 const& tVec = PL.TransVec( transList[nidx] ); + //mprintf("\tNEIGHBOR %i (idxs %i - %i)\n", nbrCell, beg1, end1); + // Loop over every atom in nbrCell + for (PairList::CellType::const_iterator it1 = nbrCell.begin(); + it1 != nbrCell.end(); ++it1) + { + Vec3 const& xyz1 = it1->ImageCoords(); + double q1 = Charge_[it1->Idx()]; + Vec3 dxyz = xyz1 + tVec - xyz0; + double rij2 = dxyz.Magnitude2(); +# ifdef DEBUG_PAIRLIST + mprintf("\t\tAtom %6i to atom %6i (%f)\n", it0->Idx()+1, it1->Idx()+1, sqrt(rij2)); +# endif + //mprintf("\t\tNbrAtom %06i\n",atnum1); + // If atom excluded, calc adjustment, otherwise calc elec. energy. + // TODO Is there better way of checking this? + if (excluded.find( it1->Idx() ) == excluded.end()) + { + //mprintf("\t\t\tdist= %f\n", sqrt(rij2)); + if ( rij2 < cut2 ) { +# include "EnergyKernel_Nonbond.h" + } + } else { +# include "EnergyKernel_Adjust.h" + } + } // END loop over neighbor cell atoms + } // END Loop over neighbor cells + } // Loop over thisCell atoms + } // END if thisCell is not empty + } // Loop over cells diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 61e567cb76..a5e29e4bea 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -230,7 +230,7 @@ EnsembleOut.o : EnsembleOut.cpp ActionFrameCounter.h ArgList.h Atom.h AtomExtra. EnsembleOutList.o : EnsembleOutList.cpp ActionFrameCounter.h ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h EnsembleOut.h EnsembleOutList.h EnsembleOut_Multi.h EnsembleOut_Single.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Topology.h TrajectoryFile.h TrajectoryIO.h Vec3.h EnsembleOut_Multi.o : EnsembleOut_Multi.cpp ActionFrameCounter.h ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h EnsembleOut.h EnsembleOut_Multi.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h Topology.h TrajectoryFile.h TrajectoryIO.h Vec3.h EnsembleOut_Single.o : EnsembleOut_Single.cpp ActionFrameCounter.h ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h EnsembleOut.h EnsembleOut_Single.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Topology.h TrajectoryFile.h TrajectoryIO.h Vec3.h -Ewald.o : Ewald.cpp Atom.h AtomExtra.h AtomMask.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajStdio.h EnergyKernel_Adjust.h EnergyKernel_Nonbond.h Ewald.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h PairList.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h SymbolExporting.h Timer.h Topology.h Vec3.h +Ewald.o : Ewald.cpp Atom.h AtomExtra.h AtomMask.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajStdio.h EnergyKernel_Adjust.h EnergyKernel_Nonbond.h Ewald.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h PairList.h PairListLoop.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h SymbolExporting.h Timer.h Topology.h Vec3.h Ewald_ParticleMesh.o : Ewald_ParticleMesh.cpp Atom.h AtomExtra.h AtomMask.h Box.h CharMask.h CoordinateInfo.h CpptrajStdio.h Ewald.h Ewald_ParticleMesh.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h PairList.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Timer.h Topology.h Vec3.h helpme_standalone.h Ewald_Regular.o : Ewald_Regular.cpp Atom.h AtomExtra.h AtomMask.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajStdio.h Ewald.h Ewald_Regular.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h PairList.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h Timer.h Topology.h Vec3.h Exec_Analyze.o : Exec_Analyze.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Cmd.h CmdList.h Command.h Control.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h Exec.h Exec_Analyze.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h VariableArray.h Vec3.h From 421cf005cc8a072f558408ba858f851d597b4c2c Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 15 Feb 2019 10:47:54 -0500 Subject: [PATCH 31/36] DRR - Cpptraj: Caculate some common prefactors. --- src/EnergyKernel_Adjust.h | 18 ++++++++++++++++++ src/EnergyKernel_Nonbond.h | 2 +- src/Ewald.cpp | 15 +++++++-------- src/Ewald.h | 2 ++ src/PairListLoop.h | 4 ++-- 5 files changed, 30 insertions(+), 11 deletions(-) diff --git a/src/EnergyKernel_Adjust.h b/src/EnergyKernel_Adjust.h index 8ea81e953a..ae7e18383b 100644 --- a/src/EnergyKernel_Adjust.h +++ b/src/EnergyKernel_Adjust.h @@ -1,4 +1,22 @@ e_adjust += Adjust(q0, q1, sqrt(rij2)); +/* + // Electrostatic exclusion adjustment +# ifndef _OPENMP + t_adjust_.Start(); + t_erfc_.Start(); +# endif + double rij = sqrt(rij2); + //double erfc = erfc_func(ew_coeff_ * rij); + double erfc = ERFC(ew_coeff_ * rij); +# ifndef _OPENMP + t_erfc_.Stop(); +# endif + double d0 = (erfc - 1.0) / rij; +# ifndef _OPENMP + t_adjust_.Stop(); +# endif + e_adjust += q0 * q1 * d0; +*/ # ifdef CPPTRAJ_EKERNEL_LJPME // LJ PME direct space correction // NOTE: Assuming excluded pair is within cutoff diff --git a/src/EnergyKernel_Nonbond.h b/src/EnergyKernel_Nonbond.h index b44a426599..d8c58bd91b 100644 --- a/src/EnergyKernel_Nonbond.h +++ b/src/EnergyKernel_Nonbond.h @@ -21,7 +21,7 @@ int nbindex = NB_->GetLJindex(TypeIndices_[it0->Idx()], TypeIndices_[it1->Idx()]); if (nbindex > -1) { - double vswitch = switch_fn(rij2, cut2_0, cut2_1); + double vswitch = switch_fn(rij2, cut2_0_, cut2_); NonbondType const& LJ = NB_->NBarray()[ nbindex ]; double r2 = 1.0 / rij2; double r6 = r2 * r2 * r2; diff --git a/src/Ewald.cpp b/src/Ewald.cpp index 04b43f849c..812c16844c 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -16,6 +16,8 @@ Ewald::Ewald() : lw_coeff_(0.0), switch_width_(0.0), cutoff_(0.0), + cut2_(0.0), + cut2_0_(0.0), dsumTol_(0.0), erfcTableDx_(0.0), one_over_Dx_(0.0), @@ -284,6 +286,11 @@ int Ewald::CheckInput(Box const& boxIn, int debugIn, double cutoffIn, double dsu else if (DABS(lw_coeff_) < Constants::SMALL) lw_coeff_ = ew_coeff_; + // Calculate some common factors. + cut2_ = cutoff_ * cutoff_; + double cut0 = cutoff_ - switch_width_; + cut2_0_ = cut0 * cut0; + return 0; } @@ -373,14 +380,10 @@ static inline double switch_fn(double rij2, double cut2_0, double cut2_1) double Ewald::Direct_VDW_LongRangeCorrection(PairList const& PL, double& e_adjust_out, double& evdw_out) { t_direct_.Start(); - double cut2 = cutoff_ * cutoff_; double Eelec = 0.0; double e_adjust = 0.0; double Evdw = 0.0; int cidx; - double cut0 = cutoff_ - switch_width_; - double cut2_0 = cut0 * cut0; - double cut2_1 = cut2; # ifdef _OPENMP # pragma omp parallel private(cidx) reduction(+: Eelec, Evdw, e_adjust) { @@ -402,16 +405,12 @@ double Ewald::Direct_VDW_LongRangeCorrection(PairList const& PL, double& e_adjus double Ewald::Direct_VDW_LJPME(PairList const& PL, double& e_adjust_out, double& evdw_out) { t_direct_.Start(); - double cut2 = cutoff_ * cutoff_; double Eelec = 0.0; double e_adjust = 0.0; double Evdw = 0.0; double Eljpme_correction = 0.0; double Eljpme_correction_excl = 0.0; int cidx; - double cut0 = cutoff_ - switch_width_; - double cut2_0 = cut0 * cut0; - double cut2_1 = cut2; # define CPPTRAJ_EKERNEL_LJPME # ifdef _OPENMP # pragma omp parallel private(cidx) reduction(+: Eelec, Evdw, e_adjust, Eljpme_correction,Eljpme_correction_excl) diff --git a/src/Ewald.h b/src/Ewald.h index fac2b33619..825f48f770 100644 --- a/src/Ewald.h +++ b/src/Ewald.h @@ -86,6 +86,8 @@ class Ewald { double lw_coeff_; ///< Ewald coefficient for LJ double switch_width_; ///< Switching window size for LJ switch if active double cutoff_; ///< Direct space cutoff + double cut2_; ///< Direct space cutoff squared. + double cut2_0_; ///< Direct space cutoff minus switch width, squared. double dsumTol_; ///< Direct space sum tolerance. double erfcTableDx_; ///< Spacing of X values in Erfc table. double one_over_Dx_; ///< One over erfcTableDx_. diff --git a/src/PairListLoop.h b/src/PairListLoop.h index 0dad69c27a..b799cc0e5c 100644 --- a/src/PairListLoop.h +++ b/src/PairListLoop.h @@ -32,7 +32,7 @@ // If atom excluded, calc adjustment, otherwise calc elec. energy. if (excluded.find( it1->Idx() ) == excluded.end()) { - if ( rij2 < cut2 ) { + if ( rij2 < cut2_ ) { # include "EnergyKernel_Nonbond.h" } } else { @@ -66,7 +66,7 @@ if (excluded.find( it1->Idx() ) == excluded.end()) { //mprintf("\t\t\tdist= %f\n", sqrt(rij2)); - if ( rij2 < cut2 ) { + if ( rij2 < cut2_ ) { # include "EnergyKernel_Nonbond.h" } } else { From 51a151ef61bdbdaf5283b2bbedbcc4b432e23ff6 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 15 Feb 2019 11:29:13 -0500 Subject: [PATCH 32/36] DRR - Cpptraj: Remove old code. --- src/EnergyKernel_Adjust.h | 18 ------- src/Ewald.cpp | 111 -------------------------------------- 2 files changed, 129 deletions(-) diff --git a/src/EnergyKernel_Adjust.h b/src/EnergyKernel_Adjust.h index ae7e18383b..8ea81e953a 100644 --- a/src/EnergyKernel_Adjust.h +++ b/src/EnergyKernel_Adjust.h @@ -1,22 +1,4 @@ e_adjust += Adjust(q0, q1, sqrt(rij2)); -/* - // Electrostatic exclusion adjustment -# ifndef _OPENMP - t_adjust_.Start(); - t_erfc_.Start(); -# endif - double rij = sqrt(rij2); - //double erfc = erfc_func(ew_coeff_ * rij); - double erfc = ERFC(ew_coeff_ * rij); -# ifndef _OPENMP - t_erfc_.Stop(); -# endif - double d0 = (erfc - 1.0) / rij; -# ifndef _OPENMP - t_adjust_.Stop(); -# endif - e_adjust += q0 * q1 * d0; -*/ # ifdef CPPTRAJ_EKERNEL_LJPME // LJ PME direct space correction // NOTE: Assuming excluded pair is within cutoff diff --git a/src/Ewald.cpp b/src/Ewald.cpp index 812c16844c..3d5f170ea5 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -448,117 +448,6 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) return Direct_VDW_LJPME(PL, e_adjust_out, evdw_out); else return Direct_VDW_LongRangeCorrection(PL, e_adjust_out, evdw_out); -/* - t_direct_.Start(); - double cut2 = cutoff_ * cutoff_; - double Eelec = 0.0; - double e_adjust = 0.0; - double Evdw = 0.0; - double Eljpme_correction = 0.0; - double Eljpme_correction_excl = 0.0; - int cidx; - - double vswitch; - double cut0 = cutoff_ - switch_width_; - double cut2_0 = cut0 * cut0; - double cut2_1 = cut2; -# ifdef _OPENMP -# pragma omp parallel private(cidx, vswitch) reduction(+: Eelec, Evdw, e_adjust,Eljpme_correction,Eljpme_correction_excl ) - { -# pragma omp for -# endif - for (cidx = 0; cidx < PL.NGridMax(); cidx++) - { - PairList::CellType const& thisCell = PL.Cell( cidx ); - if (thisCell.NatomsInGrid() > 0) - { - // cellList contains this cell index and all neighbors. - PairList::Iarray const& cellList = thisCell.CellList(); - // transList contains index to translation for the neighbor. - PairList::Iarray const& transList = thisCell.TransList(); - // Loop over all atoms of thisCell. - for (PairList::CellType::const_iterator it0 = thisCell.begin(); - it0 != thisCell.end(); ++it0) - { - Vec3 const& xyz0 = it0->ImageCoords(); - double q0 = Charge_[it0->Idx()]; -# ifdef DEBUG_PAIRLIST - mprintf("DBG: Cell %6i (%6i atoms):\n", cidx+1, thisCell.NatomsInGrid()); -# endif - // Exclusion list for this atom - Iset const& excluded = Excluded_[it0->Idx()]; - // Calc interaction of atom to all other atoms in thisCell. - for (PairList::CellType::const_iterator it1 = it0 + 1; - it1 != thisCell.end(); ++it1) - { - Vec3 const& xyz1 = it1->ImageCoords(); - double q1 = Charge_[it1->Idx()]; - Vec3 dxyz = xyz1 - xyz0; - double rij2 = dxyz.Magnitude2(); -# ifdef DEBUG_PAIRLIST - mprintf("\tAtom %6i to atom %6i (%f)\n", it0->Idx()+1, it1->Idx()+1, sqrt(rij2)); -# endif - // If atom excluded, calc adjustment, otherwise calc elec. energy. - if (excluded.find( it1->Idx() ) == excluded.end()) - { - if ( rij2 < cut2 ) { -# include "EnergyKernel_Nonbond.h" - } - } else { -# include "EnergyKernel_Adjust.h" - } - } // END loop over other atoms in thisCell - // Loop over all neighbor cells - for (unsigned int nidx = 1; nidx != cellList.size(); nidx++) - { - PairList::CellType const& nbrCell = PL.Cell( cellList[nidx] ); -# ifdef DEBUG_PAIRLIST - if (nbrCell.NatomsInGrid()>0) mprintf("\tto neighbor cell %6i\n", cellList[nidx]+1); -# endif - // Translate vector for neighbor cell - Vec3 const& tVec = PL.TransVec( transList[nidx] ); - //mprintf("\tNEIGHBOR %i (idxs %i - %i)\n", nbrCell, beg1, end1); - // Loop over every atom in nbrCell - for (PairList::CellType::const_iterator it1 = nbrCell.begin(); - it1 != nbrCell.end(); ++it1) - { - Vec3 const& xyz1 = it1->ImageCoords(); - double q1 = Charge_[it1->Idx()]; - Vec3 dxyz = xyz1 + tVec - xyz0; - double rij2 = dxyz.Magnitude2(); -# ifdef DEBUG_PAIRLIST - mprintf("\t\tAtom %6i to atom %6i (%f)\n", it0->Idx()+1, it1->Idx()+1, sqrt(rij2)); -# endif - //mprintf("\t\tNbrAtom %06i\n",atnum1); - // If atom excluded, calc adjustment, otherwise calc elec. energy. - // TODO Is there better way of checking this? - if (excluded.find( it1->Idx() ) == excluded.end()) - { - //mprintf("\t\t\tdist= %f\n", sqrt(rij2)); - if ( rij2 < cut2 ) { -# include "EnergyKernel_Nonbond.h" - } - } else { -# include "EnergyKernel_Adjust.h" - } - } // END loop over neighbor cell atoms - } // END Loop over neighbor cells - } // Loop over thisCell atoms - } // END if thisCell is not empty - } // Loop over cells -# ifdef _OPENMP - } // END pragma omp parallel -# endif - t_direct_.Stop(); - e_adjust_out = e_adjust; -# ifdef DEBUG_PAIRLIST - mprintf("DEBUG: LJ vdw = %16.8f\n", Evdw); - mprintf("DEBUG: LJ vdw PME correction = %16.8f\n", Eljpme_correction); - mprintf("DEBUG: LJ vdw PME correction (excluded) = %16.8f\n", Eljpme_correction_excl); -# endif - evdw_out = Evdw + Eljpme_correction + Eljpme_correction_excl; - return Eelec; -*/ } /** Determine VDW long range correction prefactor. */ From 8536e76e6026cbe9310e73b909de14de3fdbc152 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 15 Feb 2019 12:05:49 -0500 Subject: [PATCH 33/36] DRR - Cpptraj: Avoid allocating PME instance every call --- src/Ewald_ParticleMesh.cpp | 19 +++++++++---------- src/Ewald_ParticleMesh.h | 4 ++++ src/cpptrajdepend | 2 +- 3 files changed, 14 insertions(+), 11 deletions(-) diff --git a/src/Ewald_ParticleMesh.cpp b/src/Ewald_ParticleMesh.cpp index 0a3aaa90dd..3c6f5a3195 100644 --- a/src/Ewald_ParticleMesh.cpp +++ b/src/Ewald_ParticleMesh.cpp @@ -1,8 +1,7 @@ #ifdef LIBPME #include // copy/fill -#include // unique_ptr +//#incl ude // unique_ptr #include "Ewald_ParticleMesh.h" -#include "helpme_standalone.h" #include "CpptrajStdio.h" typedef helpme::Matrix Mat; @@ -176,8 +175,8 @@ double Ewald_ParticleMesh::Recip_ParticleMesh(Box const& boxIn) // 8 = max # threads to use for each MPI instance; 0 = all available threads used. // NOTE: Scale factor for Charmm is 332.0716 // NOTE: The electrostatic constant has been baked into the Charge_ array already. - auto pme_object = std::unique_ptr(new PMEInstanceD()); - pme_object->setup(1, ew_coeff_, order_, nfft1, nfft2, nfft3, 1.0, 0); + //auto pme_object = std::unique_ptr(new PMEInstanceD()); + pme_object_.setup(1, ew_coeff_, order_, nfft1, nfft2, nfft3, 1.0, 0); // Sets the unit cell lattice vectors, with units consistent with those used to specify coordinates. // Args: 1 = the A lattice parameter in units consistent with the coordinates. // 2 = the B lattice parameter in units consistent with the coordinates. @@ -186,10 +185,10 @@ double Ewald_ParticleMesh::Recip_ParticleMesh(Box const& boxIn) // 5 = the beta lattice parameter in degrees. // 6 = the gamma lattice parameter in degrees. // 7 = lattice type - pme_object->setLatticeVectors(boxIn.BoxX(), boxIn.BoxY(), boxIn.BoxZ(), + pme_object_.setLatticeVectors(boxIn.BoxX(), boxIn.BoxY(), boxIn.BoxZ(), boxIn.Alpha(), boxIn.Beta(), boxIn.Gamma(), PMEInstanceD::LatticeType::XAligned); - double erecip = pme_object->computeERec(0, chargesD, coordsD); + double erecip = pme_object_.computeERec(0, chargesD, coordsD); t_recip_.Stop(); return erecip; @@ -210,12 +209,12 @@ double Ewald_ParticleMesh::LJ_Recip_ParticleMesh(Box const& boxIn) Mat coordsD(&coordsD_[0], Charge_.size(), 3); Mat cparamD(&Cparam_[0], Cparam_.size(), 1); - auto pme_vdw = std::unique_ptr(new PMEInstanceD()); - pme_vdw->setup(6, lw_coeff_, order_, nfft1, nfft2, nfft3, -1.0, 0); - pme_vdw->setLatticeVectors(boxIn.BoxX(), boxIn.BoxY(), boxIn.BoxZ(), + //auto pme_vdw = std::unique_ptr(new PMEInstanceD()); + pme_vdw_.setup(6, lw_coeff_, order_, nfft1, nfft2, nfft3, -1.0, 0); + pme_vdw_.setLatticeVectors(boxIn.BoxX(), boxIn.BoxY(), boxIn.BoxZ(), boxIn.Alpha(), boxIn.Beta(), boxIn.Gamma(), PMEInstanceD::LatticeType::XAligned); - double evdwrecip = pme_vdw->computeERec(0, cparamD, coordsD); + double evdwrecip = pme_vdw_.computeERec(0, cparamD, coordsD); t_recip_.Stop(); return evdwrecip; } diff --git a/src/Ewald_ParticleMesh.h b/src/Ewald_ParticleMesh.h index f7acbff3ec..1e3b4acea8 100644 --- a/src/Ewald_ParticleMesh.h +++ b/src/Ewald_ParticleMesh.h @@ -2,6 +2,7 @@ #define INC_EWALD_PARTICLEMESH_H #ifdef LIBPME #include "Ewald.h" +#include "helpme_standalone.h" /// Class for calculating electrostatics with particle mesh Ewald. class Ewald_ParticleMesh : public Ewald { public: @@ -27,6 +28,9 @@ class Ewald_ParticleMesh : public Ewald { int nfft_[3]; ///< Number of FFT grid points in each direction int order_; ///< PME B spline order + + PMEInstanceD pme_object_; + PMEInstanceD pme_vdw_; }; #endif /* LIBPME */ #endif diff --git a/src/cpptrajdepend b/src/cpptrajdepend index a5e29e4bea..e4baf5c089 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -27,7 +27,7 @@ Action_Dihedral.o : Action_Dihedral.cpp Action.h ActionState.h Action_Dihedral.h Action_Dipole.o : Action_Dipole.cpp Action.h ActionState.h Action_Dipole.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h Action_DistRmsd.o : Action_DistRmsd.cpp Action.h ActionState.h Action_DistRmsd.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.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 Parallel.h ParameterTypes.h Range.h ReferenceAction.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h Action_Distance.o : Action_Distance.cpp Action.h ActionState.h Action_Distance.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ImagedAction.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h -Action_Energy.o : Action_Energy.cpp Action.h ActionState.h Action_Energy.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h Energy.h Ewald.h Ewald_ParticleMesh.h Ewald_Regular.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h PairList.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h +Action_Energy.o : Action_Energy.cpp Action.h ActionState.h Action_Energy.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h Energy.h Ewald.h Ewald_ParticleMesh.h Ewald_Regular.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h PairList.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h helpme_standalone.h Action_Esander.o : Action_Esander.cpp Action.h ActionState.h Action_Esander.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h Energy_Sander.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h Action_FilterByData.o : Action_FilterByData.cpp Action.h ActionState.h Action_FilterByData.h ArgList.h Array1D.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h Action_FixAtomOrder.o : Action_FixAtomOrder.cpp Action.h ActionState.h Action_FixAtomOrder.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h BondSearch.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.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 Parallel.h ParameterTypes.h ParmFile.h ParmIO.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h From fb5f1e7428794cf418a6d25cb7dbd958cfe9ed4b Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 15 Feb 2019 12:06:43 -0500 Subject: [PATCH 34/36] DRR - Cpptraj: Protect openmp pragmas --- src/helpme_standalone.h | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/helpme_standalone.h b/src/helpme_standalone.h index 3a3ea56520..1fc1e919b9 100644 --- a/src/helpme_standalone.h +++ b/src/helpme_standalone.h @@ -2881,7 +2881,9 @@ class PMEInstance { // Exclude m=0 cell. int start = (nodeZero ? 1 : 0); // Writing the three nested loops in one allows for better load balancing in parallel. +#ifdef _OPENMP #pragma omp parallel for reduction(+ : energy, Vxx, Vxy, Vyy, Vxz, Vyz, Vzz) num_threads(nThreads) +#endif for (size_t yxz = start; yxz < nyxz; ++yxz) { size_t xz = yxz % nxz; short ky = yxz / nxz; @@ -2976,7 +2978,9 @@ class PMEInstance { // Exclude m=0 cell. int start = (nodeZero ? 1 : 0); // Writing the three nested loops in one allows for better load balancing in parallel. +#ifdef _OPENMP #pragma omp parallel for num_threads(nThreads) +#endif for (size_t yxz = start; yxz < nyxz; ++yxz) { size_t xz = yxz % nxz; short ky = yxz / nxz; @@ -3587,7 +3591,9 @@ class PMEInstance { } transformedGrid[0] = Complex(0, 0); +#ifdef _OPENMP #pragma omp parallel for reduction(+ : energy) num_threads(nThreads_) +#endif for (size_t yxz = 0; yxz < nyxz; ++yxz) { size_t xz = yxz % nxz; int kx = startX + xz / nz; @@ -3765,7 +3771,9 @@ class PMEInstance { size_t rowSize = std::ceil(nForceComponents / cacheLineSizeInReals_) * cacheLineSizeInReals_; RealMat fractionalPhis(nThreads_, rowSize); size_t nAtoms = atomList_.size(); +#ifdef _OPENMP #pragma omp parallel for num_threads(nThreads_) +#endif for (size_t relativeAtomNumber = 0; relativeAtomNumber < nAtoms; ++relativeAtomNumber) { const auto &entry = splineCache_[relativeAtomNumber]; const int &atom = entry.absoluteAtomNumber; From b3c1dc44d5a8f86e009e049a278adceca7350449 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 15 Feb 2019 12:24:02 -0500 Subject: [PATCH 35/36] DRR - Cpptraj: Squash a bunch of compiler warnings. --- src/helpme_standalone.h | 37 ++++++++++++++++++++++--------------- 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/src/helpme_standalone.h b/src/helpme_standalone.h index 1fc1e919b9..1beb9ab868 100644 --- a/src/helpme_standalone.h +++ b/src/helpme_standalone.h @@ -207,7 +207,8 @@ void JacobiCyclicDiagonalization(Real *eigenvalues, Real *eigenvectors, const Re Real threshold_norm; Real threshold; Real tan_phi, sin_phi, cos_phi, tan2_phi, sin2_phi, cos2_phi; - Real sin_2phi, cos_2phi, cot_2phi; + Real sin_2phi, cot_2phi; + //Real cos_2phi Real dum1; Real dum2; Real dum3; @@ -260,7 +261,7 @@ void JacobiCyclicDiagonalization(Real *eigenvalues, Real *eigenvectors, const Re if (tan_phi < 0) sin_phi = -sin_phi; cos_phi = sqrt(cos2_phi); sin_2phi = 2 * sin_phi * cos_phi; - cos_2phi = cos2_phi - sin2_phi; + //cos_2phi = cos2_phi - sin2_phi; // Rotate columns k and m for both the matrix A // and the matrix of eigenvectors. @@ -787,8 +788,8 @@ class Matrix { */ void assertSymmetric(const Real& threshold = 1e-10f) const { assertSquare(); - for (int row = 0; row < nRows_; ++row) { - for (int col = 0; col < row; ++col) { + for (unsigned int row = 0; row < nRows_; ++row) { + for (unsigned int col = 0; col < row; ++col) { if (std::abs(data_[row * nCols_ + col] - data_[col * nCols_ + row]) > threshold) throw std::runtime_error("Unexpected non-symmetric matrix found."); } @@ -817,7 +818,7 @@ class Matrix { Matrix evecs = std::get<1>(eigenPairs); evalsReal.applyOperationToEachElement(function); Matrix evecsT = evecs.transpose(); - for (int row = 0; row < nRows_; ++row) { + for (unsigned int row = 0; row < nRows_; ++row) { Real transformedEigenvalue = evalsReal[row][0]; std::for_each(evecsT.data_ + row * nCols_, evecsT.data_ + (row + 1) * nCols_, [&](Real& val) { val *= transformedEigenvalue; }); @@ -853,10 +854,10 @@ class Matrix { throw std::runtime_error("Attempting to multiply matrices with incompatible dimensions."); Matrix product(nRows_, other.nCols_); Real* output = product.data_; - for (int row = 0; row < nRows_; ++row) { + for (unsigned int row = 0; row < nRows_; ++row) { const Real* rowPtr = data_ + row * nCols_; - for (int col = 0; col < other.nCols_; ++col) { - for (int link = 0; link < nCols_; ++link) { + for (unsigned int col = 0; col < other.nCols_; ++col) { + for (unsigned int link = 0; link < nCols_; ++link) { *output += rowPtr[link] * other.data_[link * other.nCols_ + col]; } ++output; @@ -987,10 +988,10 @@ class Matrix { unsortedEigenVectors.transposeInPlace(); std::vector> eigenPairs; - for (int val = 0; val < nRows_; ++val) eigenPairs.push_back({eigenValues[val][0], unsortedEigenVectors[val]}); + for (unsigned int val = 0; val < nRows_; ++val) eigenPairs.push_back({eigenValues[val][0], unsortedEigenVectors[val]}); std::sort(eigenPairs.begin(), eigenPairs.end()); if (order == SortOrder::Descending) std::reverse(eigenPairs.begin(), eigenPairs.end()); - for (int val = 0; val < nRows_; ++val) { + for (unsigned int val = 0; val < nRows_; ++val) { const auto& e = eigenPairs[val]; eigenValues.data_[val] = std::get<0>(e); std::copy(std::get<1>(e), std::get<1>(e) + nCols_, sortedEigenVectors[val]); @@ -2497,7 +2498,7 @@ class PMEInstance { atomList_.clear(); size_t nAtoms = coords.nRows(); - for (int atom = 0; atom < nAtoms; ++atom) { + for (unsigned int atom = 0; atom < nAtoms; ++atom) { const Real *atomCoords = coords[atom]; constexpr float EPS = 1e-6; Real aCoord = @@ -2527,11 +2528,11 @@ class PMEInstance { nAtoms = atomList_.size(); if (splineCache_.size() < nAtoms) { size_t newSize = static_cast(1.2 * nAtoms); - for (int atom = splineCache_.size(); atom < newSize; ++atom) + for (unsigned int atom = splineCache_.size(); atom < newSize; ++atom) splineCache_.emplace_back(splineOrder_, splineDerivativeLevel); } - for (int atomListNum = 0; atomListNum < nAtoms; ++atomListNum) { + for (unsigned int atomListNum = 0; atomListNum < nAtoms; ++atomListNum) { const auto &entry = atomList_[atomListNum]; const int absoluteAtomNumber = std::get<0>(entry); const Real aCoord = std::get<1>(entry); @@ -2814,7 +2815,7 @@ class PMEInstance { if (coordinates.nRows() != parameters.nRows()) throw std::runtime_error( "Inconsistent number of coordinates and parameters; there should be nAtoms of each."); - if (parameters.nCols() != nCartesian(parameterAngMom)) + if ((int)parameters.nCols() != nCartesian(parameterAngMom)) throw std::runtime_error( "Mismatch in the number of parameters provided and the parameter angular momentum"); } @@ -3284,14 +3285,20 @@ class PMEInstance { * \return Pointer to the transformed grid, which is stored in one of the buffers in BAC order. */ Complex *forwardTransform(Real *realGrid) { +# if HAVE_MPI == 1 Real *realCBA; +# endif Complex *buffer1, *buffer2; if (realGrid == reinterpret_cast(workSpace1_.data())) { +# if HAVE_MPI == 1 realCBA = reinterpret_cast(workSpace2_.data()); +# endif buffer1 = workSpace2_.data(); buffer2 = workSpace1_.data(); } else { +# if HAVE_MPI == 1 realCBA = reinterpret_cast(workSpace2_.data()); +# endif buffer1 = workSpace2_.data(); buffer2 = workSpace1_.data(); } @@ -3599,7 +3606,7 @@ class PMEInstance { int kx = startX + xz / nz; // We only loop over the first nx/2+1 x values; this // accounts for the "missing" complex conjugate values. - Real permPrefac = kx != 0 && kx != halfNx - 1 ? 2 : 1; + Real permPrefac = kx != 0 && kx != (int)halfNx - 1 ? 2 : 1; Real structFactorNorm = std::norm(transformedGrid[yxz]); energy += permPrefac * structFactorNorm * influenceFunction[yxz]; transformedGrid[yxz] *= influenceFunction[yxz]; From cf7896b787094cbb0eb557b0b33c98f8497f5be4 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 15 Feb 2019 12:38:01 -0500 Subject: [PATCH 36/36] DRR - Cpptraj: No need to pass the adjustment energy out - just report it when debug info is on. --- src/Ewald.cpp | 29 +++++++++++++++++++---------- src/Ewald.h | 6 +++--- src/Ewald_ParticleMesh.cpp | 11 +++++------ src/Ewald_Regular.cpp | 11 +++++------ 4 files changed, 32 insertions(+), 25 deletions(-) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index 3d5f170ea5..f1fe2aaf14 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -377,7 +377,10 @@ static inline double switch_fn(double rij2, double cut2_0, double cut2_1) } } -double Ewald::Direct_VDW_LongRangeCorrection(PairList const& PL, double& e_adjust_out, double& evdw_out) +/** Nonbond direct-space calculation for Coulomb electrostatics and Lennard-Jones, + * intended for use with long-range LJ correction. + */ +double Ewald::Direct_VDW_LongRangeCorrection(PairList const& PL, double& evdw_out) { t_direct_.Start(); double Eelec = 0.0; @@ -394,15 +397,19 @@ double Ewald::Direct_VDW_LongRangeCorrection(PairList const& PL, double& e_adjus } // END pragma omp parallel # endif t_direct_.Stop(); - e_adjust_out = e_adjust; # ifdef DEBUG_PAIRLIST + mprintf("DEBUG: Elec = %16.8f\n", Eelec); + mprintf("DEBUG: Eadjust = %16.8f\n", e_adjust); mprintf("DEBUG: LJ vdw = %16.8f\n", Evdw); # endif evdw_out = Evdw; - return Eelec; + return Eelec + e_adjust; } -double Ewald::Direct_VDW_LJPME(PairList const& PL, double& e_adjust_out, double& evdw_out) +/** Nonbond direct-space calculation for Coulomb electrostatics and Lennard-Jones + * calculated via PME. + */ +double Ewald::Direct_VDW_LJPME(PairList const& PL, double& evdw_out) { t_direct_.Start(); double Eelec = 0.0; @@ -423,14 +430,15 @@ double Ewald::Direct_VDW_LJPME(PairList const& PL, double& e_adjust_out, double& # endif # undef CPPTRAJ_EKERNEL_LJPME t_direct_.Stop(); - e_adjust_out = e_adjust; # ifdef DEBUG_PAIRLIST + mprintf("DEBUG: Elec = %16.8f\n", Eelec); + mprintf("DEBUG: Eadjust = %16.8f\n", e_adjust); mprintf("DEBUG: LJ vdw = %16.8f\n", Evdw); mprintf("DEBUG: LJ vdw PME correction = %16.8f\n", Eljpme_correction); mprintf("DEBUG: LJ vdw PME correction (excluded) = %16.8f\n", Eljpme_correction_excl); # endif evdw_out = Evdw + Eljpme_correction + Eljpme_correction_excl; - return Eelec; + return Eelec + e_adjust; } @@ -440,14 +448,15 @@ double Ewald::Direct_VDW_LJPME(PairList const& PL, double& e_adjust_out, double& * atoms. * \param PL The pairlist used to calculate energy. * \param e_adjust_out The electrostatic adjust energy for excluded atoms. - * \param evdw_out The direct space van der Waals term + * \param evdw_out The direct space van der Waals term (corrected for exclusion if LJ PME). + * \return The electrostatics term plus exclusion adjustment. */ -double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) +double Ewald::Direct(PairList const& PL, double& evdw_out) { if (lw_coeff_ > 0.0) - return Direct_VDW_LJPME(PL, e_adjust_out, evdw_out); + return Direct_VDW_LJPME(PL, evdw_out); else - return Direct_VDW_LongRangeCorrection(PL, e_adjust_out, evdw_out); + return Direct_VDW_LongRangeCorrection(PL, evdw_out); } /** Determine VDW long range correction prefactor. */ diff --git a/src/Ewald.h b/src/Ewald.h index 825f48f770..2355dae973 100644 --- a/src/Ewald.h +++ b/src/Ewald.h @@ -59,7 +59,7 @@ class Ewald { double Direct(Matrix_3x3 const&, Topology const&, AtomMask const&); # endif /// Fast version of direct space energy using a pairlist - double Direct(PairList const&, double&, double&); + double Direct(PairList const&, double&); /// \return adjusted energy for excluded atom pair # ifdef _OPENMP inline double Adjust(double,double,double) const; @@ -101,7 +101,7 @@ class Ewald { Timer t_erfc_; Timer t_adjust_; private: - double Direct_VDW_LongRangeCorrection(PairList const&, double&, double&); - double Direct_VDW_LJPME(PairList const&, double&, double&); + double Direct_VDW_LongRangeCorrection(PairList const&, double&); + double Direct_VDW_LJPME(PairList const&, double&); }; #endif diff --git a/src/Ewald_ParticleMesh.cpp b/src/Ewald_ParticleMesh.cpp index 3c6f5a3195..94777711c3 100644 --- a/src/Ewald_ParticleMesh.cpp +++ b/src/Ewald_ParticleMesh.cpp @@ -259,15 +259,14 @@ double Ewald_ParticleMesh::CalcEnergy(Frame const& frameIn, AtomMask const& mask e_vdw_lr_correction = Vdw_Correction( volume ); } - double e_adjust = 0.0; - e_vdw = 0.0; - double e_direct = Direct( pairList_, e_adjust, e_vdw ); + e_vdw = 0.0; + double e_direct = Direct( pairList_, e_vdw ); if (debug_ > 0) - mprintf("DEBUG: Eself= %20.10f Erecip= %20.10f Edirect= %20.10f Eadjust= %20.10f Evdw= %20.10f\n", - e_self, e_recip, e_direct, e_adjust, e_vdw); + mprintf("DEBUG: Eself= %20.10f Erecip= %20.10f Edirect= %20.10f Evdw= %20.10f\n", + e_self, e_recip, e_direct, e_vdw); e_vdw += (e_vdw_lr_correction + e_vdw6self + e_vdw6recip); t_total_.Stop(); - return e_self + e_recip + e_direct + e_adjust; + return e_self + e_recip + e_direct; } #endif /* LIBPME */ diff --git a/src/Ewald_Regular.cpp b/src/Ewald_Regular.cpp index 0f99aee0a7..f37d34f8c8 100644 --- a/src/Ewald_Regular.cpp +++ b/src/Ewald_Regular.cpp @@ -390,13 +390,12 @@ double Ewald_Regular::CalcEnergy(Frame const& frameIn, AtomMask const& maskIn, d // MapCoords(frameIn, ucell, recip, maskIn); double e_recip = Recip_Regular( recip, volume ); - double e_adjust = 0.0; - e_vdw = 0.0; - double e_direct = Direct( pairList_, e_adjust, e_vdw ); + e_vdw = 0.0; + double e_direct = Direct( pairList_, e_vdw ); if (debug_ > 0) - mprintf("DEBUG: Eself= %20.10f Erecip= %20.10f Edirect= %20.10f Eadjust= %20.10f Evdw= %20.10f\n", - e_self, e_recip, e_direct, e_adjust, e_vdw); + mprintf("DEBUG: Eself= %20.10f Erecip= %20.10f Edirect= %20.10f Evdw= %20.10f\n", + e_self, e_recip, e_direct, e_vdw); e_vdw += e_vdwr; t_total_.Stop(); - return e_self + e_recip + e_direct + e_adjust; + return e_self + e_recip + e_direct; }