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/src/Action_Energy.cpp b/src/Action_Energy.cpp index 4677a01e69..3f10507117 100644 --- a/src/Action_Energy.cpp +++ b/src/Action_Energy.cpp @@ -23,8 +23,8 @@ 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 [ewcoeff ] [skinnb ]\n" + "\t pme [cut ] [dsumtol ] [order ] [ljswidth ]\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_); + ljswidth_ = actionArgs.getKeyDouble("ljswidth", 0.0); 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,18 @@ 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 (ljswidth_ > 0.0) + mprintf("\tWidth of LJ switch region: %.4f Ang.\n", ljswidth_); + } if (KEtype_ != KE_NONE) { if (KEtype_ == KE_AUTO) mprintf("\tIf forces and velocities present KE will be calculated assuming\n" @@ -335,7 +352,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_, ljswidth_, skinnb_, erfcDx_, npoints_, + debug_, mlimits_)) return Action::ERR; EW_->Setup( setup.Top(), Imask_ ); } @@ -356,6 +374,7 @@ Action::RetType Action_Energy::Setup(ActionSetup& setup) { "Warning: 'ketype vv' to estimate kinetic energy.\n"); } } + currentParm_ = setup.TopAddress(); return Action::OK; } diff --git a/src/Action_Energy.h b/src/Action_Energy.h index 668a7acdd3..979e637cba 100644 --- a/src/Action_Energy.h +++ b/src/Action_Energy.h @@ -45,6 +45,8 @@ 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 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/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/EnergyKernel_Adjust.h b/src/EnergyKernel_Adjust.h new file mode 100644 index 0000000000..8ea81e953a --- /dev/null +++ b/src/EnergyKernel_Adjust.h @@ -0,0 +1,13 @@ + 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; + 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; +# endif diff --git a/src/EnergyKernel_Nonbond.h b/src/EnergyKernel_Nonbond.h new file mode 100644 index 0000000000..d8c58bd91b --- /dev/null +++ b/src/EnergyKernel_Nonbond.h @@ -0,0 +1,43 @@ + 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) { + double vswitch = switch_fn(rij2, cut2_0_, cut2_); + 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); +# ifdef CPPTRAJ_EKERNEL_LJPME + // 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; +# endif + } diff --git a/src/Ewald.cpp b/src/Ewald.cpp index 50537f115a..f1fe2aaf14 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -13,7 +13,11 @@ Ewald::Ewald() : sumq_(0.0), sumq2_(0.0), ew_coeff_(0.0), + 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), @@ -177,6 +181,21 @@ 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(); + 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. */ void Ewald::SetupExcluded(Topology const& topIn, AtomMask const& maskIn) { @@ -219,12 +238,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) { @@ -243,6 +265,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) @@ -251,7 +278,18 @@ 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_; + + // Calculate some common factors. + cut2_ = cutoff_ * cutoff_; + double cut0 = cutoff_ - switch_width_; + cut2_0_ = cut0 * cut0; return 0; } @@ -280,7 +318,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_; @@ -293,6 +331,18 @@ double Ewald::Self(double volume) { return ene; } +/** Lennard-Jones self energy. */ +double Ewald::Self6() { + t_self_.Start(); // TODO precalc + 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) + 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 { @@ -313,15 +363,26 @@ double Ewald::Adjust(double q0, double q1, double rij) { } # endif -// Ewald::Direct() -/** Calculate direct space energy. This is the faster version that uses - * a pair list. Also calculate the energy adjustment for excluded - * atoms. +/** Switching function for Lennard-Jones. */ +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); + } +} + +/** Nonbond direct-space calculation for Coulomb electrostatics and Lennard-Jones, + * intended for use with long-range LJ correction. */ -double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) +double Ewald::Direct_VDW_LongRangeCorrection(PairList const& PL, double& evdw_out) { t_direct_.Start(); - double cut2 = cutoff_ * cutoff_; double Eelec = 0.0; double e_adjust = 0.0; double Evdw = 0.0; @@ -331,153 +392,71 @@ double Ewald::Direct(PairList const& PL, double& e_adjust_out, double& evdw_out) { # 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 ) { - 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) { - 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; - //mprintf("PVDW %8i%8i%20.6f%20.6f\n", ta0+1, ta1+1, e_vdw, r2); - } - } - } else - e_adjust += Adjust(q0, q1, sqrt(rij2)); - } // 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 ) { - 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) { - 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; - //mprintf("PVDW %8i%8i%20.6f%20.6f\n", ta0+1, ta1+1, e_vdw, r2); - } - } - } else - e_adjust += Adjust(q0, q1, sqrt(rij2)); - } // END loop over neighbor cell atoms - } // END Loop over neighbor cells - } // Loop over thisCell atoms - } // END if thisCell is not empty - } // Loop over cells +# include "PairListLoop.h" # ifdef _OPENMP } // 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; +} + +/** 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; + double e_adjust = 0.0; + double Evdw = 0.0; + double Eljpme_correction = 0.0; + double Eljpme_correction_excl = 0.0; + int cidx; +# 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(); +# 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 + e_adjust; +} + + +// Ewald::Direct() +/** 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 (corrected for exclusion if LJ PME). + * \return The electrostatics term plus exclusion adjustment. + */ +double Ewald::Direct(PairList const& PL, double& evdw_out) +{ + if (lw_coeff_ > 0.0) + return Direct_VDW_LJPME(PL, evdw_out); + else + return Direct_VDW_LongRangeCorrection(PL, evdw_out); } /** Determine VDW long range correction prefactor. */ diff --git a/src/Ewald.h b/src/Ewald.h index f4458c7974..2355dae973 100644 --- a/src/Ewald.h +++ b/src/Ewald.h @@ -37,14 +37,18 @@ 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 - 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 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 @@ -55,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; @@ -68,6 +72,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. @@ -77,19 +82,26 @@ 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 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 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_. 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 Direct_VDW_LJPME(PairList const&, double&); }; #endif diff --git a/src/Ewald_ParticleMesh.cpp b/src/Ewald_ParticleMesh.cpp index e5ce9eca5f..94777711c3 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; @@ -88,10 +87,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 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, 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_); @@ -105,6 +106,10 @@ 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_); + 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 "); @@ -127,6 +132,9 @@ 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); + // 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); SetupExcluded(topIn, maskIn); @@ -167,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. @@ -177,14 +185,40 @@ 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; } +/** 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, 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); + t_recip_.Stop(); + return evdwrecip; +} + /** Calculate full nonbonded energy with PME */ double Ewald_ParticleMesh::CalcEnergy(Frame const& frameIn, AtomMask const& maskIn, double& e_vdw) { @@ -192,7 +226,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 ); + double e_vdw_lr_correction; pairList_.CreatePairList(frameIn, ucell, recip, maskIn); @@ -208,15 +242,31 @@ double Ewald_ParticleMesh::CalcEnergy(Frame const& frameIn, AtomMask const& mask // MapCoords(frameIn, ucell, recip, maskIn); double e_recip = Recip_ParticleMesh( frameIn.BoxCrd() ); - double e_adjust = 0.0; - e_vdw = 0.0; - double e_direct = Direct( pairList_, e_adjust, e_vdw ); + + // TODO branch + double e_vdw6self, e_vdw6recip; + if (lw_coeff_ > 0.0) { + e_vdw6self = Self6(); + 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; + e_vdw6recip = 0.0; + e_vdw_lr_correction = Vdw_Correction( volume ); + } + + 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); - e_vdw += e_vdwr; + 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_ParticleMesh.h b/src/Ewald_ParticleMesh.h index 73c625bbf7..1e3b4acea8 100644 --- a/src/Ewald_ParticleMesh.h +++ b/src/Ewald_ParticleMesh.h @@ -2,12 +2,14 @@ #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: 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, 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? @@ -19,11 +21,16 @@ 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 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/Ewald_Regular.cpp b/src/Ewald_Regular.cpp index 0a713040a6..f37d34f8c8 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; @@ -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. @@ -388,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; } diff --git a/src/Exec_CatCrd.cpp b/src/Exec_CatCrd.cpp new file mode 100644 index 0000000000..623fecd280 --- /dev/null +++ b/src/Exec_CatCrd.cpp @@ -0,0 +1,93 @@ +#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 ); + } + } + dsarg = argIn.GetStringNext(); + } + 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 diff --git a/src/PairListLoop.h b/src/PairListLoop.h new file mode 100644 index 0000000000..b799cc0e5c --- /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/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; 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 diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 178553e614..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 @@ -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 @@ -230,11 +230,12 @@ 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 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 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 \ diff --git a/src/helpme_standalone.h b/src/helpme_standalone.h index 3a3ea56520..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"); } @@ -2881,7 +2882,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 +2979,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; @@ -3280,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(); } @@ -3587,13 +3598,15 @@ 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; // 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]; @@ -3765,7 +3778,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; diff --git a/test/Makefile b/test/Makefile index f2e7ef3939..8406e62b0c 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) @@ -443,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 \ @@ -585,7 +589,8 @@ COMPLETETESTS=test.general \ test.cphstats \ test.cmdline \ test.xyzfmt \ - test.hausdorff + test.hausdorff \ + test.catcrd test.all: $(MAKE) test.complete summary 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 < lj.in < lj.in < lj.in <