Skip to content

Commit 18fb2a8

Browse files
drroeDaniel R. Roe
andauthored
Revert #1159, LJ 12-6-4 term needs more testing (#1160)
* Revert #1159. Needs more testing * Revert #1159 * Revert #1159 * Revert #1159 * Revert #1159 * Revert #1159 * Revert #1159 * Revert #1159 * Revert #1159 * Revert #1159 --------- Co-authored-by: Daniel R. Roe <daniel.roe@nih.gov>
1 parent f109533 commit 18fb2a8

10 files changed

Lines changed: 32 additions & 62 deletions

src/Energy.cpp

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -175,7 +175,7 @@ double Energy_Amber::Calc_14_Energy(Frame const& fIn, DihedralArray const& Dihed
175175
double rij = sqrt( rij2 );
176176
// VDW
177177
NonbondType const& LJ = tIn.GetLJparam(d->A1(), d->A4());
178-
double e_vdw = Cpptraj::Energy::Ene_LJ_12_6_4( rij2, LJ.A(), LJ.B(), LJ.C() );
178+
double e_vdw = Cpptraj::Energy::Ene_LJ_6_12( rij2, LJ.A(), LJ.B() );
179179
/* double r2 = 1.0 / rij2;
180180
double r6 = r2 * r2 * r2;
181181
double r12 = r6 * r6;
@@ -235,7 +235,7 @@ double Energy_Amber::E_Nonbond(Frame const& fIn, Topology const& tIn, AtomMask c
235235
double rij = sqrt( rij2 );
236236
// VDW
237237
NonbondType const& LJ = tIn.GetLJparam(atom1, atom2);
238-
double e_vdw = Cpptraj::Energy::Ene_LJ_12_6_4( rij2, LJ.A(), LJ.B(), LJ.C() );
238+
double e_vdw = Cpptraj::Energy::Ene_LJ_6_12( rij2, LJ.A(), LJ.B() );
239239
/*double r2 = 1.0 / rij2;
240240
double r6 = r2 * r2 * r2;
241241
double r12 = r6 * r6;
@@ -295,7 +295,7 @@ double Energy_Amber::E_VDW(Frame const& fIn, Topology const& tIn, AtomMask const
295295
double rij2 = DIST2_NoImage( crd1, fIn.XYZ( atom2 ) );
296296
// VDW
297297
NonbondType const& LJ = tIn.GetLJparam(atom1, atom2);
298-
double e_vdw = Cpptraj::Energy::Ene_LJ_12_6_4( rij2, LJ.A(), LJ.B(), LJ.C() );
298+
double e_vdw = Cpptraj::Energy::Ene_LJ_6_12( rij2, LJ.A(), LJ.B() );
299299
/*double r2 = 1.0 / rij2;
300300
double r6 = r2 * r2 * r2;
301301
double r12 = r6 * r6;

src/Energy/Ene_Decomp_Nonbond.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,7 @@ void Ene_Decomp_Nonbond(Frame const& fIn, Topology const& tIn, CharMask const& s
3838
T rij2 = (dx*dx) + (dy*dy) + (dz*dz);
3939
// VDW
4040
NonbondType const& LJ = tIn.GetLJparam(atom1, atom2);
41-
T e_vdw = Ene_LJ_12_6_4<T>( rij2, LJ.A(), LJ.B(), LJ.C() );
41+
T e_vdw = Ene_LJ_6_12<T>( rij2, LJ.A(), LJ.B() );
4242
# ifdef CPPTRAJ_DEBUG_ENEDECOMP
4343
mprintf("DEBUG: VDW %f\n", e_vdw);
4444
# endif

src/Energy/Ene_LJ_6_12.h

Lines changed: 0 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -14,21 +14,6 @@ T Ene_LJ_6_12(T const& rij2, T const& LJA, T const& LJB)
1414
T e_vdw = f12 - f6; // (A/r^12)-(B/r^6)
1515
return e_vdw;
1616
}
17-
18-
/// \return LJ 12-6-4 energy
19-
template <typename T>
20-
T Ene_LJ_12_6_4(T const& rij2, T const& LJA, T const& LJB, T const& LJC)
21-
{
22-
T r2 = 1.0 / rij2;
23-
T r4 = r2 * r2;
24-
T r6 = r4 * r2;
25-
T r12 = r6 * r6;
26-
T f12 = LJA * r12; // A/r^12
27-
T f6 = LJB * r6; // B/r^6
28-
T f4 = LJC * r4; // C/r^4
29-
T e_vdw = f12 - f6 - f4; // (A/r^12)-(B/r^6)-(C/r^4)
30-
return e_vdw;
31-
}
3217
}
3318
}
3419
#endif

src/Energy/EnergyDecomposer.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -370,7 +370,7 @@ void EnergyDecomposer::calcDihedrals( Frame const& frameIn ) {
370370
// 1-4 vdw energy
371371
double rij2 = DIST2_NoImage( frameIn.XYZ(dih->A1()), frameIn.XYZ(dih->A4()) );
372372
NonbondType const& LJ = currentTop_->GetLJparam(dih->A1(), dih->A4());
373-
double e_vdw = Ene_LJ_12_6_4( rij2, LJ.A(), LJ.B(), LJ.C() );
373+
double e_vdw = Ene_LJ_6_12( rij2, LJ.A(), LJ.B() );
374374
e_vdw /= DP.SCNB();
375375
# ifdef CPPTRAJ_DEBUG_ENEDECOMP
376376
mprintf("DEBUG: V14 %f\n", e_vdw);

src/EnergyKernel_NonBond_Simple.h

Lines changed: 5 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -12,14 +12,14 @@ template <class REAL> class EnergyKernel_NonBond_Simple {
1212
/// CONSTRUCTOR - specified cutoff^2
1313
EnergyKernel_NonBond_Simple(REAL cutIn) : cutoff2_(cutIn) {}
1414
/// Energy/forces
15-
void Calc_F_E(Frame&, int, int, double, double, double, double, double, double, double, double, CharMask const&, double&, double&);
15+
void Calc_F_E(Frame&, int, int, double, double, double, double, double, double, double, CharMask const&, double&, double&);
1616
private:
1717
REAL cutoff2_; ///< Only calculate interactions within this cutoff (squared)
1818
};
1919

2020
template<class REAL>
2121
void EnergyKernel_NonBond_Simple<REAL>::Calc_F_E(Frame& frameIn, int idx, int jdx,
22-
double LJA, double LJB, double LJC,
22+
double LJA, double LJB,
2323
double QFAC, double qi, double qj,
2424
double enbfac, double eelfac,
2525
CharMask const& maskIn,
@@ -35,13 +35,11 @@ void EnergyKernel_NonBond_Simple<REAL>::Calc_F_E(Frame& frameIn, int idx, int jd
3535
REAL rij = sqrt( rij2 );
3636
// VDW
3737
REAL r2 = 1.0 / rij2;
38-
REAL r4 = r2 * r2;
39-
REAL r6 = r4 * r2;
38+
REAL r6 = r2 * r2 * r2;
4039
REAL r12 = r6 * r6;
4140
REAL f12 = LJA * r12; // A/r^12
4241
REAL f6 = LJB * r6; // B/r^6
43-
REAL f4 = LJC * r4; // C/r^4
44-
REAL e_vdw = f12 - f6 - f4; // (A/r^12)-(B/r^6)-(C/r^4)
42+
REAL e_vdw = f12 - f6; // (A/r^12)-(B/r^6)
4543
//mprintf("DBG:\t\t%8i %8i %12.4f\n", e_vdw);
4644
E_vdw += (e_vdw * enbfac);
4745
// VDW force
@@ -58,11 +56,7 @@ void EnergyKernel_NonBond_Simple<REAL>::Calc_F_E(Frame& frameIn, int idx, int jd
5856
// dfy += ry * felec;
5957
// dfz += rz * felec;
6058
// COMBINED Vdw and Coulomb force
61-
// Force = -dU/dr
62-
// U = A/r^12 - B/r^6 - C/r^4
63-
// F = 12A/r^13 - 6B/r^7 - 4C/r^5
64-
// F_term = F * r = 12A/r^12 - 6B/r^6 - 4C/r^4 = 12f12 - 6f6 - 4f4
65-
REAL dfn = ( (((12*f12) - (6*f6) - (4*f4)) * enbfac) + e_elec ) * r2;
59+
REAL dfn = ( (((12*f12) - (6*f6)) * enbfac) + e_elec ) * r2;
6660
REAL dfx = rx * dfn;
6761
REAL dfy = ry * dfn;
6862
REAL dfz = rz * dfn;

src/PairListEngine_Ewald_Decomp_LJLR.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,7 @@ class PairListEngine_Ewald_Decomp_LJLR {
4545
if (nbindex > -1) {
4646
double vswitch = EW_.Switch_Fn(rij2);
4747
NonbondType const& LJ = EW_.GetLJ( nbindex );
48-
T e_vdw = Cpptraj::Energy::Ene_LJ_12_6_4<T>(rij2, LJ.A(), LJ.B(), LJ.C());
48+
T e_vdw = Cpptraj::Energy::Ene_LJ_6_12<T>(rij2, LJ.A(), LJ.B());
4949
e_vdw *= vswitch;
5050
Evdw_ += e_vdw;
5151
e_half = e_vdw * 0.5;

src/ParameterTypes.h

Lines changed: 13 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -279,15 +279,12 @@ class NonbondType {
279279
# define tol_ 0.00000001
280280
//static const double tol_ = 0.00000001;
281281
public:
282-
NonbondType() : A_(0), B_(0), C_(0) {}
283-
NonbondType(double a, double b, double c=0.0) : A_(a), B_(b), C_(c) {}
282+
NonbondType() : A_(0), B_(0) {}
283+
NonbondType(double a, double b) : A_(a), B_(b) {}
284284
inline double A() const { return A_; }
285285
inline double B() const { return B_; }
286-
inline double C() const { return C_; }
287286
void SetA(double a) { A_ = a; }
288287
void SetB(double b) { B_ = b; }
289-
void SetC(double c) { C_ = c; }
290-
291288
double Radius() const {
292289
if (B_ > 0.0)
293290
return (0.5 * pow(2.0 * A_ / B_, (1.0/6.0)));
@@ -303,31 +300,26 @@ class NonbondType {
303300
/// \return True if A and B match
304301
bool operator==(NonbondType const& rhs) const {
305302
return ( (fabs(A_ - rhs.A_) < tol_) &&
306-
(fabs(B_ - rhs.B_) < tol_) &&
307-
(fabs(C_ - rhs.C_) < tol_) );
303+
(fabs(B_ - rhs.B_) < tol_) );
308304
}
309305
/// \return True if A and B do not match
310306
bool operator!=(NonbondType const& rhs) const {
311307
return ( (fabs(A_ - rhs.A_) > tol_) ||
312-
(fabs(B_ - rhs.B_) > tol_) ||
313-
(fabs(C_ - rhs.C_) > tol_) );
308+
(fabs(B_ - rhs.B_) > tol_) );
314309
}
315310
/// \return True if A less than zero, or B if A is equal.
316311
bool operator<(NonbondType const& rhs) const {
317312
if (*this != rhs) {
318-
if ( (fabs(A_ - rhs.A_) < tol_) ) {
319-
if (fabs(B_ - rhs.B_) < tol_)
320-
return (C_ < rhs.C_);
313+
if ( (fabs(A_ - rhs.A_) < tol_) )
321314
return (B_ < rhs.B_);
322-
} else
315+
else
323316
return (A_ < rhs.A_);
324317
} else
325318
return false;
326319
}
327320
private:
328321
double A_;
329322
double B_;
330-
double C_;
331323
# undef tol_
332324
};
333325
typedef std::vector<NonbondType> NonbondArray;
@@ -379,20 +371,20 @@ typedef std::vector<LJparmType> LJparmArray;
379371
*/
380372
class NonbondParmType {
381373
public:
382-
NonbondParmType() : ntypes_(0), has_c_param_(false) {}
374+
NonbondParmType() : ntypes_(0) {}
383375
// NonbondParmType(int n, std::vector<int> const& nbi, NonbondArray const& nba,
384376
// HB_ParmArray const& hba) :
385377
// ntypes_(n), nbindex_(nbi), nbarray_(nba), hbarray_(hba) {}
386378
inline bool HasNonbond() const { return ntypes_ > 0; }
387379
inline int Ntypes() const { return ntypes_; }
388-
bool Has_C_Coeff() const { return has_c_param_; }
380+
bool Has_C_Coeff() const { return !ccoef_.empty(); }
389381
std::vector<int> const& NBindex() const { return nbindex_; }
390382
/// \return Array of LJ 12-6 A and B parameters
391383
NonbondArray const& NBarray() const { return nbarray_; }
392384
/// \return Array of LJ 10-12 (hbond) parameters
393385
HB_ParmArray const& HBarray() const { return hbarray_; }
394386
/// \return Array of LJ 12-6-4 C parameters
395-
//std::vector<double> const& LJC_Array() const { return ccoef_; }
387+
std::vector<double> const& LJC_Array() const { return ccoef_; }
396388
/// \return LJ 6-12 A and B parameter at specified index
397389
NonbondType const& NBarray(int i) const { return nbarray_[i]; }
398390
/// \return LJ 10-12 (hbond) parameter at specified index
@@ -414,8 +406,8 @@ class NonbondParmType {
414406
void SetNHBterms(int n) { hbarray_.assign( n, HB_ParmType() ); }
415407
/// Set specified HB term
416408
HB_ParmType& SetHB(int i) { return hbarray_[i]; }
417-
/// Set a LJ C parameter
418-
void SetLJC(int i, double c) { nbarray_[i].SetC(c); has_c_param_ = true; }
409+
/// Add a LJ C parameter
410+
void AddLJC(double c) { ccoef_.push_back( c ); }
419411
/// Set specified nbindex location to given value.
420412
void SetNbIdx(int idx, int nbidx) { nbindex_[idx] = nbidx; }
421413
/// Add given LJ term to nonbond array and update nonbond index array.
@@ -447,13 +439,13 @@ class NonbondParmType {
447439
nbindex_[ntypes_ * type2 + type1] = ndx;
448440
hbarray_.push_back( HB );
449441
}
450-
void Clear() { ntypes_ = 0; nbindex_.clear(); nbarray_.clear(); hbarray_.clear(); has_c_param_ = false; }
442+
void Clear() { ntypes_ = 0; nbindex_.clear(); nbarray_.clear(); hbarray_.clear(); }
451443
private:
452444
int ntypes_; ///< Number of unique atom types
453445
std::vector<int> nbindex_; ///< Hold indices into arrays nbarray/hbarray for atom type pairs
454446
NonbondArray nbarray_; ///< Hold Lennard-Jones 6-12 A and B parameters for all pairs.
455447
HB_ParmArray hbarray_; ///< Hold 10-12 Amber HBond params for all pairs.
456-
bool has_c_param_; ///< True if Lennard-Jones C parameters for 12-6-4 LJ potential present.
448+
std::vector<double> ccoef_; ///< Hold Lennard-Jones C parameters for 12-6-4 LJ potential.
457449
};
458450
// ----- LES PARAMETERS --------------------------------------------------------
459451
/// Hold LES atom parameters

src/Parm_Amber.cpp

Lines changed: 5 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -890,7 +890,7 @@ int Parm_Amber::ReadLJC(Topology& TopIn, FortranData const& FMT) {
890890
if (SetupBuffer(F_LJ_C, numLJparm_, FMT)) return 1;
891891
for (int idx = 0; idx != numLJparm_; idx++)
892892
{
893-
TopIn.SetNonbond().SetLJC( idx, FileBufferToDouble(F_LJ_C, idx, numLJparm_) );
893+
TopIn.SetNonbond().AddLJC( FileBufferToDouble(F_LJ_C, idx, numLJparm_) );
894894
if (atProblemFlag_) break;
895895
}
896896
return 0;
@@ -2094,7 +2094,6 @@ int Parm_Amber::WriteParm(FileName const& fname, Topology const& TopOut) {
20942094

20952095
// LJ A and B terms
20962096
if (WriteLJ(F_LJ_A, F_LJ_B, TopOut.Nonbond().NBarray())) return 1;
2097-
20982097

20992098
// CHAMBER only - LJ 1-4 terms
21002099
if (ptype_ == CHAMBER) {
@@ -2481,10 +2480,10 @@ int Parm_Amber::WriteParm(FileName const& fname, Topology const& TopOut) {
24812480
// LJ C terms.
24822481
// NOTE: I put this at the very end to match behavior of parmed add_12_6_4
24832482
if (TopOut.Nonbond().Has_C_Coeff()) {
2484-
if (BufferAlloc(F_LJ_C, TopOut.Nonbond().NBarray().size())) return 1;
2485-
for (NonbondArray::const_iterator it = TopOut.Nonbond().NBarray().begin();
2486-
it != TopOut.Nonbond().NBarray().end(); ++it)
2487-
file_.DblToBuffer( it->C() );
2483+
std::vector<double> const& ccoef = TopOut.Nonbond().LJC_Array();
2484+
if (BufferAlloc(F_LJ_C, ccoef.size())) return 1;
2485+
for (std::vector<double>::const_iterator it = ccoef.begin(); it != ccoef.end(); ++it)
2486+
file_.DblToBuffer( *it );
24882487
file_.FlushBuffer();
24892488
}
24902489

src/PotentialTerm_Dihedral.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,7 @@ void PotentialTerm_Dihedral::CalcForce(Frame& frameIn, CharMask const& maskIn) c
7777
Atom const& A4 = (*atoms_)[dih->A4()];
7878
int nbindex = nonbond_->GetLJindex( A1.TypeIndex(), A4.TypeIndex() );
7979
NonbondType const& LJ = nonbond_->NBarray( nbindex );
80-
NB14.Calc_F_E( frameIn, dih->A1(), dih->A4(), LJ.A(), LJ.B(), LJ.C(),
80+
NB14.Calc_F_E( frameIn, dih->A1(), dih->A4(), LJ.A(), LJ.B(),
8181
Constants::COULOMBFACTOR, A1.Charge(), A4.Charge(),
8282
1.0/DP.SCNB(), 1.0/DP.SCEE(), maskIn,
8383
*Enb14_, *Eq14_);

src/PotentialTerm_LJ_Coulomb.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -96,7 +96,7 @@ void PotentialTerm_LJ_Coulomb::CalcForce(Frame& frameIn, CharMask const& maskIn)
9696
{
9797
NonbondType LJ = GetLJparam(*atoms_, *nonbond_, *idx, *jdx);
9898
//NonBondKernel(frameIn, *idx, *jdx, LJ.A(), LJ.B(), maskIn, *E_vdw_, *E_elec_);
99-
nonbond.Calc_F_E(frameIn, *idx, *jdx, LJ.A(), LJ.B(), LJ.C(), QFAC_, (*atoms_)[*idx].Charge(), (*atoms_)[*jdx].Charge(), 1, 1, maskIn, *E_vdw_, *E_elec_);
99+
nonbond.Calc_F_E(frameIn, *idx, *jdx, LJ.A(), LJ.B(), QFAC_, (*atoms_)[*idx].Charge(), (*atoms_)[*jdx].Charge(), 1, 1, maskIn, *E_vdw_, *E_elec_);
100100
} // END idx not bonded to jdx
101101
} // END inner loop over jdx
102102
} // END outer loop over idx
@@ -112,7 +112,7 @@ void PotentialTerm_LJ_Coulomb::CalcForce(Frame& frameIn, CharMask const& maskIn)
112112
if (excluded.find( *idx1 ) == excluded.end())
113113
{
114114
NonbondType LJ = GetLJparam(*atoms_, *nonbond_, *idx0, *idx1);
115-
nonbond.Calc_F_E(frameIn, *idx0, *idx1, LJ.A(), LJ.B(), LJ.C(), QFAC_, (*atoms_)[*idx0].Charge(), (*atoms_)[*idx1].Charge(), 1, 1, maskIn, *E_vdw_, *E_elec_);
115+
nonbond.Calc_F_E(frameIn, *idx0, *idx1, LJ.A(), LJ.B(), QFAC_, (*atoms_)[*idx0].Charge(), (*atoms_)[*idx1].Charge(), 1, 1, maskIn, *E_vdw_, *E_elec_);
116116
//NonBondKernel(frameIn, *idx0, *idx1, LJ.A(), LJ.B(), maskIn, *E_vdw_, *E_elec_);
117117
}
118118
} // END inner loop over idx1

0 commit comments

Comments
 (0)