diff --git a/src/ASPHERE/pair_gayberne.cpp b/src/ASPHERE/pair_gayberne.cpp index 25bdae14f1..9ff87326ed 100644 --- a/src/ASPHERE/pair_gayberne.cpp +++ b/src/ASPHERE/pair_gayberne.cpp @@ -391,7 +391,7 @@ double PairGayBerne::init_one(int i, int j) lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double ratio = sigma[i][j] / cut[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/ASPHERE/pair_resquared.cpp b/src/ASPHERE/pair_resquared.cpp index ed9d9b36c4..caa031a1e8 100644 --- a/src/ASPHERE/pair_resquared.cpp +++ b/src/ASPHERE/pair_resquared.cpp @@ -391,7 +391,7 @@ double PairRESquared::init_one(int i, int j) lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double ratio = sigma[i][j] / cut[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/CLASS2/pair_lj_class2.cpp b/src/CLASS2/pair_lj_class2.cpp index e79dc0c6de..0b90b2717e 100644 --- a/src/CLASS2/pair_lj_class2.cpp +++ b/src/CLASS2/pair_lj_class2.cpp @@ -235,7 +235,7 @@ double PairLJClass2::init_one(int i, int j) lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0); lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double ratio = sigma[i][j] / cut[i][j]; offset[i][j] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.0*pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/CLASS2/pair_lj_class2_coul_cut.cpp b/src/CLASS2/pair_lj_class2_coul_cut.cpp index bec7f1da15..395953e0a9 100644 --- a/src/CLASS2/pair_lj_class2_coul_cut.cpp +++ b/src/CLASS2/pair_lj_class2_coul_cut.cpp @@ -286,7 +286,7 @@ double PairLJClass2CoulCut::init_one(int i, int j) lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0); lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; offset[i][j] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.0*pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/CLASS2/pair_lj_class2_coul_long.cpp b/src/CLASS2/pair_lj_class2_coul_long.cpp index 5f7d738e92..5d2ae891d0 100644 --- a/src/CLASS2/pair_lj_class2_coul_long.cpp +++ b/src/CLASS2/pair_lj_class2_coul_long.cpp @@ -329,7 +329,7 @@ double PairLJClass2CoulLong::init_one(int i, int j) lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0); lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; offset[i][j] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.0*pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/COLLOID/pair_colloid.cpp b/src/COLLOID/pair_colloid.cpp index 68150f6eff..983b973e0e 100644 --- a/src/COLLOID/pair_colloid.cpp +++ b/src/COLLOID/pair_colloid.cpp @@ -354,7 +354,7 @@ double PairColloid::init_one(int i, int j) lj4[j][i] = lj4[i][j] = 4.0 * epsilon * sigma6[i][j]; offset[j][i] = offset[i][j] = 0.0; - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double tmp; offset[j][i] = offset[i][j] = single(0,0,i,j,cut[i][j]*cut[i][j],0.0,1.0,tmp); diff --git a/src/COLLOID/pair_yukawa_colloid.cpp b/src/COLLOID/pair_yukawa_colloid.cpp index 9b8d0ecaad..87fa7f5422 100644 --- a/src/COLLOID/pair_yukawa_colloid.cpp +++ b/src/COLLOID/pair_yukawa_colloid.cpp @@ -147,7 +147,7 @@ double PairYukawaColloid::init_one(int i, int j) cut[i][j] = mix_distance(cut[i][i],cut[j][j]); } - if (offset_flag) { + if (offset_flag && (kappa != 0.0)) { double screening = exp(-kappa * (cut[i][j] - (rad[i]+rad[j]))); offset[i][j] = a[i][j]/kappa * screening; } else offset[i][j] = 0.0; diff --git a/src/DIPOLE/pair_lj_cut_dipole_cut.cpp b/src/DIPOLE/pair_lj_cut_dipole_cut.cpp index addd02e505..af987bf258 100644 --- a/src/DIPOLE/pair_lj_cut_dipole_cut.cpp +++ b/src/DIPOLE/pair_lj_cut_dipole_cut.cpp @@ -387,7 +387,7 @@ double PairLJCutDipoleCut::init_one(int i, int j) lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/DIPOLE/pair_lj_cut_dipole_long.cpp b/src/DIPOLE/pair_lj_cut_dipole_long.cpp index 78922e356f..63f48bea81 100644 --- a/src/DIPOLE/pair_lj_cut_dipole_long.cpp +++ b/src/DIPOLE/pair_lj_cut_dipole_long.cpp @@ -420,7 +420,7 @@ double PairLJCutDipoleLong::init_one(int i, int j) lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/DIPOLE/pair_lj_long_dipole_long.cpp b/src/DIPOLE/pair_lj_long_dipole_long.cpp index 15ac2e788c..b833b250d4 100644 --- a/src/DIPOLE/pair_lj_long_dipole_long.cpp +++ b/src/DIPOLE/pair_lj_long_dipole_long.cpp @@ -314,7 +314,7 @@ double PairLJLongDipoleLong::init_one(int i, int j) //if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) //error->all(FLERR,"Pair cutoff < Respa interior cutoff"); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/KSPACE/pair_born_coul_long.cpp b/src/KSPACE/pair_born_coul_long.cpp index e588a30b55..479128ef2b 100644 --- a/src/KSPACE/pair_born_coul_long.cpp +++ b/src/KSPACE/pair_born_coul_long.cpp @@ -311,7 +311,7 @@ double PairBornCoulLong::init_one(int i, int j) born2[i][j] = 6.0*c[i][j]; born3[i][j] = 8.0*d[i][j]; - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double rexp = exp((sigma[i][j]-cut_lj[i][j])*rhoinv[i][j]); offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut_lj[i][j],6.0) + d[i][j]/pow(cut_lj[i][j],8.0); diff --git a/src/KSPACE/pair_buck_coul_long.cpp b/src/KSPACE/pair_buck_coul_long.cpp index 476e3c716a..95496409b9 100644 --- a/src/KSPACE/pair_buck_coul_long.cpp +++ b/src/KSPACE/pair_buck_coul_long.cpp @@ -297,7 +297,7 @@ double PairBuckCoulLong::init_one(int i, int j) buck1[i][j] = a[i][j]/rho[i][j]; buck2[i][j] = 6.0*c[i][j]; - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double rexp = exp(-cut_lj[i][j]/rho[i][j]); offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut_lj[i][j],6.0); } else offset[i][j] = 0.0; diff --git a/src/KSPACE/pair_buck_long_coul_long.cpp b/src/KSPACE/pair_buck_long_coul_long.cpp index 8aa4d72083..4cfb9b7267 100644 --- a/src/KSPACE/pair_buck_long_coul_long.cpp +++ b/src/KSPACE/pair_buck_long_coul_long.cpp @@ -330,7 +330,7 @@ double PairBuckLongCoulLong::init_one(int i, int j) if (cut_respa && MIN(cut_buck[i][j],cut_coul) < cut_respa[3]) error->all(FLERR,"Pair cutoff < Respa interior cutoff"); - if (offset_flag) { + if (offset_flag && (cut_buck[i][j] > 0.0)) { double rexp = exp(-cut_buck[i][j]/buck_rho[i][j]); offset[i][j] = buck_a[i][j]*rexp - buck_c[i][j]/pow(cut_buck[i][j],6.0); } else offset[i][j] = 0.0; diff --git a/src/KSPACE/pair_lj_cut_coul_long.cpp b/src/KSPACE/pair_lj_cut_coul_long.cpp index e9799843fc..f8be9fdb79 100644 --- a/src/KSPACE/pair_lj_cut_coul_long.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long.cpp @@ -743,7 +743,7 @@ double PairLJCutCoulLong::init_one(int i, int j) lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/KSPACE/pair_lj_long_coul_long.cpp b/src/KSPACE/pair_lj_long_coul_long.cpp index 44256a9fbb..7c6adfcb41 100644 --- a/src/KSPACE/pair_lj_long_coul_long.cpp +++ b/src/KSPACE/pair_lj_long_coul_long.cpp @@ -333,7 +333,7 @@ double PairLJLongCoulLong::init_one(int i, int j) if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) error->all(FLERR,"Pair cutoff < Respa interior cutoff"); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/MISC/pair_nm_cut.cpp b/src/MISC/pair_nm_cut.cpp index 0163cdcf58..b9bf6ac477 100644 --- a/src/MISC/pair_nm_cut.cpp +++ b/src/MISC/pair_nm_cut.cpp @@ -243,7 +243,7 @@ double PairNMCut::init_one(int i, int j) r0n[i][j] = pow(r0[i][j],nn[i][j]); r0m[i][j] = pow(r0[i][j],mm[i][j]); - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { offset[i][j] = e0nm[i][j] * ((mm[i][j]*r0n[i][j] / pow(cut[i][j],nn[i][j])) - (nn[i][j]*r0m[i][j] / pow(cut[i][j],mm[i][j]))); diff --git a/src/MISC/pair_nm_cut_coul_cut.cpp b/src/MISC/pair_nm_cut_coul_cut.cpp index 5cb2452906..78c77a648f 100644 --- a/src/MISC/pair_nm_cut_coul_cut.cpp +++ b/src/MISC/pair_nm_cut_coul_cut.cpp @@ -291,7 +291,7 @@ double PairNMCutCoulCut::init_one(int i, int j) r0n[i][j] = pow(r0[i][j],nn[i][j]); r0m[i][j] = pow(r0[i][j],mm[i][j]); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { offset[i][j] = e0nm[i][j] * ((mm[i][j]*r0n[i][j] / pow(cut_lj[i][j],nn[i][j])) - (nn[i][j]*r0m[i][j] / pow(cut_lj[i][j],mm[i][j]))); diff --git a/src/MISC/pair_nm_cut_coul_long.cpp b/src/MISC/pair_nm_cut_coul_long.cpp index 15d5d03757..8e0da40eac 100644 --- a/src/MISC/pair_nm_cut_coul_long.cpp +++ b/src/MISC/pair_nm_cut_coul_long.cpp @@ -339,7 +339,7 @@ double PairNMCutCoulLong::init_one(int i, int j) r0n[i][j] = pow(r0[i][j],nn[i][j]); r0m[i][j] = pow(r0[i][j],mm[i][j]); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { offset[i][j] = e0nm[i][j] * ((mm[i][j]*r0n[i][j] / pow(cut_lj[i][j],nn[i][j])) - (nn[i][j]*r0m[i][j] / pow(cut_lj[i][j],mm[i][j]))); diff --git a/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp b/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp index e3093e4d10..c616a9fa83 100644 --- a/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp +++ b/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp @@ -531,7 +531,7 @@ double PairLJCutTIP4PCut::init_one(int i, int j) lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/USER-CGSDK/pair_lj_sdk.cpp b/src/USER-CGSDK/pair_lj_sdk.cpp index 23b0f47a6d..56e56c9605 100644 --- a/src/USER-CGSDK/pair_lj_sdk.cpp +++ b/src/USER-CGSDK/pair_lj_sdk.cpp @@ -311,7 +311,7 @@ double PairLJSDK::init_one(int i, int j) lj3[i][j] = lj_prefact[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow1[ljt]); lj4[i][j] = lj_prefact[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow2[ljt]); - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double ratio = sigma[i][j] / cut[i][j]; offset[i][j] = lj_prefact[ljt] * epsilon[i][j] * (pow(ratio,lj_pow1[ljt]) - pow(ratio,lj_pow2[ljt])); } else offset[i][j] = 0.0; diff --git a/src/USER-CGSDK/pair_lj_sdk_coul_long.cpp b/src/USER-CGSDK/pair_lj_sdk_coul_long.cpp index 845c5822a7..4c9e42f778 100644 --- a/src/USER-CGSDK/pair_lj_sdk_coul_long.cpp +++ b/src/USER-CGSDK/pair_lj_sdk_coul_long.cpp @@ -401,7 +401,7 @@ double PairLJSDKCoulLong::init_one(int i, int j) lj3[i][j] = lj_prefact[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow1[ljt]); lj4[i][j] = lj_prefact[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow2[ljt]); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; offset[i][j] = lj_prefact[ljt] * epsilon[i][j] * (pow(ratio,lj_pow1[ljt]) - pow(ratio,lj_pow2[ljt])); diff --git a/src/USER-FEP/pair_lj_charmm_coul_long_soft.cpp b/src/USER-FEP/pair_lj_charmm_coul_long_soft.cpp index 81b82e9774..d8bfd698be 100644 --- a/src/USER-FEP/pair_lj_charmm_coul_long_soft.cpp +++ b/src/USER-FEP/pair_lj_charmm_coul_long_soft.cpp @@ -994,7 +994,7 @@ double PairLJCharmmCoulLongSoft::single(int i, int j, int itype, int jtype, if (rsq < cut_ljsq) { philj = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] * - (1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype]; + (1.0/(denlj*denlj) - 1.0/denlj); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; diff --git a/src/USER-FEP/pair_lj_cut_coul_cut_soft.cpp b/src/USER-FEP/pair_lj_cut_coul_cut_soft.cpp index b2e781c57b..9e52a16a02 100644 --- a/src/USER-FEP/pair_lj_cut_coul_cut_soft.cpp +++ b/src/USER-FEP/pair_lj_cut_coul_cut_soft.cpp @@ -236,6 +236,8 @@ void PairLJCutCoulCutSoft::coeff(int narg, char **arg) double epsilon_one = force->numeric(FLERR,arg[2]); double sigma_one = force->numeric(FLERR,arg[3]); double lambda_one = force->numeric(FLERR,arg[4]); + if (sigma_one <= 0.0) + error->all(FLERR,"Incorrect args for pair coefficients"); double cut_lj_one = cut_lj_global; double cut_coul_one = cut_coul_global; @@ -296,7 +298,7 @@ double PairLJCutCoulCutSoft::init_one(int i, int j) lj3[i][j] = alphalj * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]); lj4[i][j] = alphac * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double denlj = lj3[i][j] + pow(cut_lj[i][j] / sigma[i][j], 6.0); offset[i][j] = lj1[i][j] * 4.0 * epsilon[i][j] * (1.0/(denlj*denlj) - 1.0/denlj); } else offset[i][j] = 0.0; diff --git a/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp b/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp index 3b80729b0b..f7c4084fe2 100644 --- a/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp +++ b/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp @@ -604,6 +604,8 @@ void PairLJCutCoulLongSoft::coeff(int narg, char **arg) double epsilon_one = force->numeric(FLERR,arg[2]); double sigma_one = force->numeric(FLERR,arg[3]); double lambda_one = force->numeric(FLERR,arg[4]); + if (sigma_one <= 0.0) + error->all(FLERR,"Incorrect args for pair coefficients"); double cut_lj_one = cut_lj_global; if (narg == 6) cut_lj_one = force->numeric(FLERR,arg[5]); @@ -723,7 +725,7 @@ double PairLJCutCoulLongSoft::init_one(int i, int j) lj3[i][j] = alphalj * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]); lj4[i][j] = alphac * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double denlj = lj3[i][j] + pow(cut_lj[i][j] / sigma[i][j], 6.0); offset[i][j] = lj1[i][j] * 4.0 * epsilon[i][j] * (1.0/(denlj*denlj) - 1.0/denlj); } else offset[i][j] = 0.0; diff --git a/src/USER-FEP/pair_lj_cut_soft.cpp b/src/USER-FEP/pair_lj_cut_soft.cpp index 800fdfcde8..8b6280a61a 100644 --- a/src/USER-FEP/pair_lj_cut_soft.cpp +++ b/src/USER-FEP/pair_lj_cut_soft.cpp @@ -484,6 +484,8 @@ void PairLJCutSoft::coeff(int narg, char **arg) double epsilon_one = force->numeric(FLERR,arg[2]); double sigma_one = force->numeric(FLERR,arg[3]); double lambda_one = force->numeric(FLERR,arg[4]); + if (sigma_one <= 0.0) + error->all(FLERR,"Incorrect args for pair coefficients"); double cut_one = cut_global; if (narg == 6) cut_one = force->numeric(FLERR,arg[5]); @@ -584,7 +586,7 @@ double PairLJCutSoft::init_one(int i, int j) lj2[i][j] = pow(sigma[i][j], 6.0); lj3[i][j] = alphalj * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]); - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double denlj = lj3[i][j] + pow(cut[i][j] / sigma[i][j], 6.0); offset[i][j] = lj1[i][j] * 4.0 * epsilon[i][j] * (1.0/(denlj*denlj) - 1.0/denlj); } else offset[i][j] = 0.0; diff --git a/src/USER-MISC/pair_buck_mdf.cpp b/src/USER-MISC/pair_buck_mdf.cpp index 6c3dcbd7ee..372fbaf8c0 100644 --- a/src/USER-MISC/pair_buck_mdf.cpp +++ b/src/USER-MISC/pair_buck_mdf.cpp @@ -259,7 +259,7 @@ double PairBuckMDF::init_one(int i, int j) buck1[i][j] = a[i][j]/rho[i][j]; buck2[i][j] = 6.0*c[i][j]; - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double rexp = exp(-cut[i][j]/rho[i][j]); offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut[i][j],6.0); } else offset[i][j] = 0.0; diff --git a/src/USER-MISC/pair_coul_diel.cpp b/src/USER-MISC/pair_coul_diel.cpp index a62362aa6f..c653e9abb2 100644 --- a/src/USER-MISC/pair_coul_diel.cpp +++ b/src/USER-MISC/pair_coul_diel.cpp @@ -235,7 +235,7 @@ double PairCoulDiel::init_one(int i, int j) double *q = atom->q; double qqrd2e = force->qqrd2e; - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double rarg = (cut[i][j]-rme[i][j])/sigmae[i][j]; double epsr=a_eps+b_eps*tanh(rarg); offset[i][j] = qqrd2e*q[i]*q[j]*((eps_s/epsr) -1.)/cut[i][j]; diff --git a/src/USER-MISC/pair_gauss_cut.cpp b/src/USER-MISC/pair_gauss_cut.cpp index 3836187a64..4408545c34 100644 --- a/src/USER-MISC/pair_gauss_cut.cpp +++ b/src/USER-MISC/pair_gauss_cut.cpp @@ -196,6 +196,9 @@ void PairGaussCut::coeff(int narg, char **arg) double hgauss_one = force->numeric(FLERR,arg[2]); double rmh_one = force->numeric(FLERR,arg[3]); double sigmah_one = force->numeric(FLERR,arg[4]); + if (sigmah_one <= 0.0) + error->all(FLERR,"Incorrect args for pair coefficients"); + double cut_one = cut_global; if (narg == 6) cut_one = force->numeric(FLERR,arg[5]); diff --git a/src/USER-MISC/pair_kolmogorov_crespi_z.cpp b/src/USER-MISC/pair_kolmogorov_crespi_z.cpp index 15a325e106..c03aee6362 100644 --- a/src/USER-MISC/pair_kolmogorov_crespi_z.cpp +++ b/src/USER-MISC/pair_kolmogorov_crespi_z.cpp @@ -53,7 +53,7 @@ PairKolmogorovCrespiZ::PairKolmogorovCrespiZ(LAMMPS *lmp) : Pair(lmp) map = NULL; // always compute energy offset - offset_flag = true; + offset_flag = 1; } /* ---------------------------------------------------------------------- */ @@ -285,7 +285,7 @@ double PairKolmogorovCrespiZ::init_one(int i, int j) { if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { int iparam_ij = elem2param[map[i]][map[j]]; Param& p = params[iparam_ij]; offset[i][j] = -p.A*pow(p.z0/cut[i][j],6); diff --git a/src/USER-MISC/pair_list.cpp b/src/USER-MISC/pair_list.cpp index 3adbe8b69d..458c228d58 100644 --- a/src/USER-MISC/pair_list.cpp +++ b/src/USER-MISC/pair_list.cpp @@ -81,7 +81,7 @@ void PairList::compute(int eflag, int vflag) { if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = eflag_global = - vflag_global = eflag_atom = vflag_atom = 0; + vflag_global = eflag_atom = vflag_atom = 0; const int nlocal = atom->nlocal; const int newton_pair = force->newton_pair; @@ -126,46 +126,46 @@ void PairList::compute(int eflag, int vflag) const double r2inv = 1.0/rsq; if (style[n] == HARM) { - const double r = sqrt(rsq); - const double dr = par.parm.harm.r0 - r; - fpair = 2.0*par.parm.harm.k*dr/r; + const double r = sqrt(rsq); + const double dr = par.parm.harm.r0 - r; + fpair = 2.0*par.parm.harm.k*dr/r; - if (eflag_either) - epair = par.parm.harm.k*dr*dr - par.offset; + if (eflag_either) + epair = par.parm.harm.k*dr*dr - par.offset; } else if (style[n] == MORSE) { - const double r = sqrt(rsq); - const double dr = par.parm.morse.r0 - r; - const double dexp = exp(par.parm.morse.alpha * dr); - fpair = 2.0*par.parm.morse.d0*par.parm.morse.alpha - * (dexp*dexp - dexp) / r; + const double r = sqrt(rsq); + const double dr = par.parm.morse.r0 - r; + const double dexp = exp(par.parm.morse.alpha * dr); + fpair = 2.0*par.parm.morse.d0*par.parm.morse.alpha + * (dexp*dexp - dexp) / r; - if (eflag_either) - epair = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp) - par.offset; + if (eflag_either) + epair = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp) - par.offset; } else if (style[n] == LJ126) { - const double r6inv = r2inv*r2inv*r2inv; - const double sig6 = mypow(par.parm.lj126.sigma,6); - fpair = 24.0*par.parm.lj126.epsilon*r6inv - * (2.0*sig6*sig6*r6inv - sig6) * r2inv; + const double r6inv = r2inv*r2inv*r2inv; + const double sig6 = mypow(par.parm.lj126.sigma,6); + fpair = 24.0*par.parm.lj126.epsilon*r6inv + * (2.0*sig6*sig6*r6inv - sig6) * r2inv; - if (eflag_either) - epair = 4.0*par.parm.lj126.epsilon*r6inv - * (sig6*sig6*r6inv - sig6) - par.offset; + if (eflag_either) + epair = 4.0*par.parm.lj126.epsilon*r6inv + * (sig6*sig6*r6inv - sig6) - par.offset; } if (newton_pair || i < nlocal) { - f[i].x += dx*fpair; - f[i].y += dy*fpair; - f[i].z += dz*fpair; + f[i].x += dx*fpair; + f[i].y += dy*fpair; + f[i].z += dz*fpair; } if (newton_pair || j < nlocal) { - f[j].x -= dx*fpair; - f[j].y -= dy*fpair; - f[j].z -= dz*fpair; + f[j].x -= dx*fpair; + f[j].y -= dy*fpair; + f[j].z -= dz*fpair; } if (evflag) ev_tally(i,j,nlocal,newton_pair,epair,0.0,fpair,dx,dy,dz); @@ -174,10 +174,10 @@ void PairList::compute(int eflag, int vflag) if (vflag_fdotr) virial_fdotr_compute(); if (check_flag) { - int tmp; - MPI_Allreduce(&pc,&tmp,1,MPI_INT,MPI_SUM,world); - if (tmp != 2*npairs) - error->all(FLERR,"Not all pairs processed in pair_style list"); + int tmp; + MPI_Allreduce(&pc,&tmp,1,MPI_INT,MPI_SUM,world); + if (tmp != 2*npairs) + error->all(FLERR,"Not all pairs processed in pair_style list"); } } @@ -263,12 +263,12 @@ void PairList::settings(int narg, char **arg) ptr = strtok(NULL," \t\n\r\f"); if ((ptr == NULL) || (*ptr == '#')) - error->all(FLERR,"Incorrectly formatted harmonic pair parameters"); + error->all(FLERR,"Incorrectly formatted harmonic pair parameters"); par.parm.harm.k = force->numeric(FLERR,ptr); ptr = strtok(NULL," \t\n\r\f"); if ((ptr == NULL) || (*ptr == '#')) - error->all(FLERR,"Incorrectly formatted harmonic pair parameters"); + error->all(FLERR,"Incorrectly formatted harmonic pair parameters"); par.parm.harm.r0 = force->numeric(FLERR,ptr); ++nharm; @@ -279,17 +279,17 @@ void PairList::settings(int narg, char **arg) ptr = strtok(NULL," \t\n\r\f"); if (!ptr) - error->all(FLERR,"Incorrectly formatted morse pair parameters"); + error->all(FLERR,"Incorrectly formatted morse pair parameters"); par.parm.morse.d0 = force->numeric(FLERR,ptr); ptr = strtok(NULL," \t\n\r\f"); if (!ptr) - error->all(FLERR,"Incorrectly formatted morse pair parameters"); + error->all(FLERR,"Incorrectly formatted morse pair parameters"); par.parm.morse.alpha = force->numeric(FLERR,ptr); ptr = strtok(NULL," \t\n\r\f"); if (!ptr) - error->all(FLERR,"Incorrectly formatted morse pair parameters"); + error->all(FLERR,"Incorrectly formatted morse pair parameters"); par.parm.morse.r0 = force->numeric(FLERR,ptr); ++nmorse; @@ -300,12 +300,12 @@ void PairList::settings(int narg, char **arg) ptr = strtok(NULL," \t\n\r\f"); if (!ptr) - error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters"); + error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters"); par.parm.lj126.epsilon = force->numeric(FLERR,ptr); ptr = strtok(NULL," \t\n\r\f"); if (!ptr) - error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters"); + error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters"); par.parm.lj126.sigma = force->numeric(FLERR,ptr); ++nlj126; @@ -332,10 +332,10 @@ void PairList::settings(int narg, char **arg) if (comm->me == 0) { if (screen) fprintf(screen,"Read %d (%d/%d/%d) interacting pairs from %s\n", - npairs, nharm, nmorse, nlj126, arg[0]); + npairs, nharm, nmorse, nlj126, arg[0]); if (logfile) fprintf(logfile,"Read %d (%d/%d/%d) interacting pairs from %s\n", - npairs, nharm, nmorse, nlj126, arg[0]); + npairs, nharm, nmorse, nlj126, arg[0]); } } @@ -380,18 +380,18 @@ void PairList::init_style() list_parm_t &par = params[n]; if (style[n] == HARM) { - const double dr = sqrt(par.cutsq) - par.parm.harm.r0; - par.offset = par.parm.harm.k*dr*dr; + const double dr = sqrt(par.cutsq) - par.parm.harm.r0; + par.offset = par.parm.harm.k*dr*dr; } else if (style[n] == MORSE) { - const double dr = par.parm.morse.r0 - sqrt(par.cutsq); - const double dexp = exp(par.parm.morse.alpha * dr); - par.offset = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp); + const double dr = par.parm.morse.r0 - sqrt(par.cutsq); + const double dexp = exp(par.parm.morse.alpha * dr); + par.offset = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp); } else if (style[n] == LJ126) { - const double r6inv = par.cutsq*par.cutsq*par.cutsq; - const double sig6 = mypow(par.parm.lj126.sigma,6); - par.offset = 4.0*par.parm.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6); + const double r6inv = par.cutsq*par.cutsq*par.cutsq; + const double sig6 = mypow(par.parm.lj126.sigma,6); + par.offset = 4.0*par.parm.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6); } } } diff --git a/src/pair_born.cpp b/src/pair_born.cpp index 6d420fb36b..979499488e 100644 --- a/src/pair_born.cpp +++ b/src/pair_born.cpp @@ -243,7 +243,7 @@ double PairBorn::init_one(int i, int j) born2[i][j] = 6.0*c[i][j]; born3[i][j] = 8.0*d[i][j]; - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double rexp = exp((sigma[i][j]-cut[i][j])*rhoinv[i][j]); offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut[i][j],6.0) + d[i][j]/pow(cut[i][j],8.0); diff --git a/src/pair_born_coul_dsf.cpp b/src/pair_born_coul_dsf.cpp index caec95759a..d53e9cf00e 100644 --- a/src/pair_born_coul_dsf.cpp +++ b/src/pair_born_coul_dsf.cpp @@ -306,7 +306,7 @@ double PairBornCoulDSF::init_one(int i, int j) born2[i][j] = 6.0*c[i][j]; born3[i][j] = 8.0*d[i][j]; - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double rexp = exp((sigma[i][j]-cut_lj[i][j])*rhoinv[i][j]); offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut_lj[i][j],6.0) + d[i][j]/pow(cut_lj[i][j],8.0); diff --git a/src/pair_born_coul_wolf.cpp b/src/pair_born_coul_wolf.cpp index bad0c5ed3e..fbec4f3da7 100644 --- a/src/pair_born_coul_wolf.cpp +++ b/src/pair_born_coul_wolf.cpp @@ -305,7 +305,7 @@ double PairBornCoulWolf::init_one(int i, int j) born2[i][j] = 6.0*c[i][j]; born3[i][j] = 8.0*d[i][j]; - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double rexp = exp((sigma[i][j]-cut_lj[i][j])*rhoinv[i][j]); offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut_lj[i][j],6.0) + d[i][j]/pow(cut_lj[i][j],8.0); diff --git a/src/pair_buck.cpp b/src/pair_buck.cpp index e4da772e0a..ac7ac0aeb5 100644 --- a/src/pair_buck.cpp +++ b/src/pair_buck.cpp @@ -230,7 +230,7 @@ double PairBuck::init_one(int i, int j) buck1[i][j] = a[i][j]/rho[i][j]; buck2[i][j] = 6.0*c[i][j]; - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double rexp = exp(-cut[i][j]/rho[i][j]); offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut[i][j],6.0); } else offset[i][j] = 0.0; diff --git a/src/pair_buck_coul_cut.cpp b/src/pair_buck_coul_cut.cpp index c052c3100a..2764a8e4e4 100644 --- a/src/pair_buck_coul_cut.cpp +++ b/src/pair_buck_coul_cut.cpp @@ -281,7 +281,7 @@ double PairBuckCoulCut::init_one(int i, int j) buck1[i][j] = a[i][j]/rho[i][j]; buck2[i][j] = 6.0*c[i][j]; - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double rexp = exp(-cut_lj[i][j]/rho[i][j]); offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut_lj[i][j],6.0); } else offset[i][j] = 0.0; diff --git a/src/pair_lj96_cut.cpp b/src/pair_lj96_cut.cpp index f4b2747d40..83fc5bcdda 100644 --- a/src/pair_lj96_cut.cpp +++ b/src/pair_lj96_cut.cpp @@ -557,7 +557,7 @@ double PairLJ96Cut::init_one(int i, int j) lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],9.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double ratio = sigma[i][j] / cut[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,9.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/pair_lj_cut.cpp b/src/pair_lj_cut.cpp index a3ebf414c9..7f838061f1 100644 --- a/src/pair_lj_cut.cpp +++ b/src/pair_lj_cut.cpp @@ -551,7 +551,7 @@ double PairLJCut::init_one(int i, int j) lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double ratio = sigma[i][j] / cut[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/pair_lj_cut_coul_cut.cpp b/src/pair_lj_cut_coul_cut.cpp index 0d62c43dc3..85cf9dc97d 100644 --- a/src/pair_lj_cut_coul_cut.cpp +++ b/src/pair_lj_cut_coul_cut.cpp @@ -278,7 +278,7 @@ double PairLJCutCoulCut::init_one(int i, int j) lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/pair_lj_cut_coul_dsf.cpp b/src/pair_lj_cut_coul_dsf.cpp index 09293a6f4c..5d95d06ed7 100644 --- a/src/pair_lj_cut_coul_dsf.cpp +++ b/src/pair_lj_cut_coul_dsf.cpp @@ -303,7 +303,7 @@ double PairLJCutCoulDSF::init_one(int i, int j) lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/pair_lj_expand.cpp b/src/pair_lj_expand.cpp index 2fd780472a..785733dccc 100644 --- a/src/pair_lj_expand.cpp +++ b/src/pair_lj_expand.cpp @@ -240,7 +240,7 @@ double PairLJExpand::init_one(int i, int j) lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double ratio = sigma[i][j] / cut[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; diff --git a/src/pair_mie_cut.cpp b/src/pair_mie_cut.cpp index 312fb7bc70..320f21248d 100644 --- a/src/pair_mie_cut.cpp +++ b/src/pair_mie_cut.cpp @@ -575,7 +575,7 @@ double PairMIECut::init_one(int i, int j) mie3[i][j] = Cmie[i][j] * epsilon[i][j] * pow(sigma[i][j],gamR[i][j]); mie4[i][j] = Cmie[i][j] * epsilon[i][j] * pow(sigma[i][j],gamA[i][j]); - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double ratio = sigma[i][j] / cut[i][j]; offset[i][j] = Cmie[i][j] * epsilon[i][j] * (pow(ratio,gamR[i][j]) - pow(ratio,gamA[i][j])); diff --git a/src/pair_yukawa.cpp b/src/pair_yukawa.cpp index 0e5fd36cd6..2ba6633d9e 100644 --- a/src/pair_yukawa.cpp +++ b/src/pair_yukawa.cpp @@ -210,7 +210,7 @@ double PairYukawa::init_one(int i, int j) cut[i][j] = mix_distance(cut[i][i],cut[j][j]); } - if (offset_flag) { + if (offset_flag && (cut[i][j] > 0.0)) { double screening = exp(-kappa * cut[i][j]); offset[i][j] = a[i][j] * screening / cut[i][j]; } else offset[i][j] = 0.0;