Add support for Lennard-Jones PME to 'energy' command.

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) {

src/Action_Energy.cpp

@@ -23,8 +23,8 @@ void Action_Energy::Help() const {
           "\t          ewald [cut <cutoff>] [dsumtol <dtol>] [rsumtol <rtol>]\n"
           "\t                [ewcoeff <coeff>] [maxexp <max>] [skinnb <skinnb>]\n"
           "\t                [mlimits <X>,<Y>,<Z>] [erfcdx <dx>]\n"
-          "\t          pme [cut <cutoff>] [dsumtol <dtol>] [order <order>]\n"
-          "\t              [ewcoeff <coeff>] [skinnb <skinnb>]\n"
+          "\t          pme [cut <cutoff>] [dsumtol <dtol>] [order <order>] [ljswidth <width>]\n"
+          "\t              [ewcoeff <coeff>] [ljpme] [ewcoefflj] [skinnb <skinnb>]\n"
           "\t              [nfft <nfft1>,<nfft2>,<nfft3>] [erfcdx <dx>]\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;
 }

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)

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" );

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

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
+                }