sfree(td->f);
}
-static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr)
+static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr,
+ gmx_bool b14only)
{
/* Fill the table according to the formulas in the manual.
* In principle, we only need the potential and the second
int i;
double reppow, p;
double r1, rc, r12, r13;
- double r, r2, r6, rc6;
+ double r, r2, r6, rc2, rc6, rc12;
double expr, Vtab, Ftab;
/* Parameters for David's function */
double A = 0, B = 0, C = 0, A_3 = 0, B_4 = 0;
/* Parameters for the switching function */
double ksw, swi, swi1;
/* Temporary parameters */
- gmx_bool bSwitch, bShift;
+ gmx_bool bPotentialSwitch, bForceSwitch, bPotentialShift;
double ewc = fr->ewaldcoeff_q;
double ewclj = fr->ewaldcoeff_lj;
+ double Vcut = 0;
- bSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
- (tp == etabCOULSwitch) ||
- (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch));
-
- bShift = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
- (tp == etabShift));
+ if (b14only)
+ {
+ bPotentialSwitch = FALSE;
+ bForceSwitch = FALSE;
+ bPotentialShift = FALSE;
+ }
+ else
+ {
+ bPotentialSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
+ (tp == etabCOULSwitch) ||
+ (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch) ||
+ (tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSWITCH)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSWITCH)));
+ bForceSwitch = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
+ (tp == etabShift) ||
+ (tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodFORCESWITCH)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodFORCESWITCH)));
+ bPotentialShift = ((tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSHIFT)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSHIFT)));
+ }
reppow = fr->reppow;
r1 = fr->rvdw_switch;
rc = fr->rvdw;
}
- if (bSwitch)
+ if (bPotentialSwitch)
{
ksw = 1.0/(pow5(rc-r1));
}
{
ksw = 0.0;
}
- if (bShift)
+ if (bForceSwitch)
{
if (tp == etabShift)
{
fp = xvgropen("switch.xvg", "switch", "r", "s");
#endif
+ if (bPotentialShift)
+ {
+ rc2 = rc*rc;
+ rc6 = 1.0/(rc2*rc2*rc2);
+ if (gmx_within_tol(reppow, 12.0, 10*GMX_DOUBLE_EPS))
+ {
+ rc12 = rc6*rc6;
+ }
+ else
+ {
+ rc12 = pow(rc, -reppow);
+ }
+
+ switch (tp)
+ {
+ case etabLJ6:
+ /* Dispersion */
+ Vcut = -rc6;
+ break;
+ case etabLJ6Ewald:
+ Vcut = -rc6*exp(-ewclj*ewclj*rc2)*(1 + ewclj*ewclj*rc2 + pow4(ewclj)*rc2*rc2/2);
+ break;
+ case etabLJ12:
+ /* Repulsion */
+ Vcut = rc12;
+ break;
+ case etabCOUL:
+ Vcut = 1.0/rc;
+ break;
+ case etabEwald:
+ case etabEwaldSwitch:
+ Vtab = gmx_erfc(ewc*rc)/rc;
+ break;
+ case etabEwaldUser:
+ /* Only calculate minus the reciprocal space contribution */
+ Vtab = -gmx_erf(ewc*rc)/rc;
+ break;
+ case etabRF:
+ case etabRF_ZERO:
+ /* No need for preventing the usage of modifiers with RF */
+ Vcut = 0.0;
+ break;
+ case etabEXPMIN:
+ Vcut = exp(-rc);
+ break;
+ default:
+ gmx_fatal(FARGS, "Cannot apply new potential-shift modifier to interaction type '%s' yet. (%s,%d)",
+ tprops[tp].name, __FILE__, __LINE__);
+ }
+ }
+
for (i = td->nx0; (i < td->nx); i++)
{
r = td->x[i];
}
Vtab = 0.0;
Ftab = 0.0;
- if (bSwitch)
+ if (bPotentialSwitch)
{
/* swi is function, swi1 1st derivative and swi2 2nd derivative */
/* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
gmx_fatal(FARGS, "Table type %d not implemented yet. (%s,%d)",
tp, __FILE__, __LINE__);
}
- if (bShift)
+ if (bForceSwitch)
{
/* Normal coulomb with cut-off correction for potential */
if (r < rc)
Ftab += A*r12 + B*r13;
}
}
+ else
+ {
+ /* Make sure interactions are zero outside cutoff with modifiers */
+ Vtab = 0;
+ Ftab = 0;
+ }
+ }
+ if (bPotentialShift)
+ {
+ if (r < rc)
+ {
+ Vtab -= Vcut;
+ }
+ else
+ {
+ /* Make sure interactions are zero outside cutoff with modifiers */
+ Vtab = 0;
+ Ftab = 0;
+ }
}
if (ETAB_USER(tp))
Ftab += td->f[i];
}
- if ((r > r1) && bSwitch)
+ if (bPotentialSwitch)
{
- Ftab = Ftab*swi - Vtab*swi1;
- Vtab = Vtab*swi;
+ if (r >= rc)
+ {
+ /* Make sure interactions are zero outside cutoff with modifiers */
+ Vtab = 0;
+ Ftab = 0;
+ }
+ else if (r > r1)
+ {
+ Ftab = Ftab*swi - Vtab*swi1;
+ Vtab = Vtab*swi;
+ }
}
-
/* Convert to single precision when we store to mem */
td->v[i] = Vtab;
td->f[i] = Ftab;
gmx_incons("Potential modifiers other than potential-shift are only implemented for LJ cut-off");
}
- switch (fr->vdw_modifier)
+ /* LJ-PME and other (shift-only) modifiers are handled by applying the modifiers
+ * to the original interaction forms when we fill the table, so we only check cutoffs here.
+ */
+ if (fr->vdwtype == evdwCUT)
{
- case eintmodNONE:
- case eintmodPOTSHIFT:
- case eintmodEXACTCUTOFF:
- /* No modification */
- break;
- case eintmodPOTSWITCH:
- tabsel[etiLJ6] = etabLJ6Switch;
- tabsel[etiLJ12] = etabLJ12Switch;
- break;
- case eintmodFORCESWITCH:
- tabsel[etiLJ6] = etabLJ6Shift;
- tabsel[etiLJ12] = etabLJ12Shift;
- break;
- default:
- gmx_incons("Unsupported vdw_modifier");
+ switch (fr->vdw_modifier)
+ {
+ case eintmodNONE:
+ case eintmodPOTSHIFT:
+ case eintmodEXACTCUTOFF:
+ /* No modification */
+ break;
+ case eintmodPOTSWITCH:
+ tabsel[etiLJ6] = etabLJ6Switch;
+ tabsel[etiLJ12] = etabLJ12Switch;
+ break;
+ case eintmodFORCESWITCH:
+ tabsel[etiLJ6] = etabLJ6Shift;
+ tabsel[etiLJ12] = etabLJ12Shift;
+ break;
+ default:
+ gmx_incons("Unsupported vdw_modifier");
+ }
}
}
}
init_table(nx, nx0,
(tabsel[k] == etabEXPMIN) ? table.scale_exp : table.scale,
&(td[k]), !bReadTab);
- fill_table(&(td[k]), tabsel[k], fr);
+ fill_table(&(td[k]), tabsel[k], fr, b14only);
if (out)
{
fprintf(out, "%s table with %d data points for %s%s.\n"
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff