non-bonded forces must be calculated. The pair list contains particles
$i$, a displacement vector for particle $i$, and all particles $j$ that
are within \verb'rlist' of this particular image of particle $i$. The
-list is updated every \verb'nstlist' steps, where \verb'nstlist' is
-typically 10. There is an option to calculate the total non-bonded
-force on each particle due to all particle in a shell around the
-list cut-off, {\ie} at a distance between \verb'rlist' and
-\verb'rlistlong'. This force is calculated during the pair list update
-and retained during \verb'nstlist' steps.
+list is updated every \verb'nstlist' steps.
To make the \normindex{neighbor list}, all particles that are close
({\ie} within the neighbor list cut-off) to a given particle must be found.
with double precision may be required to get fine details of
thermodynamics correct.
-\subsection{Twin-range cut-offs\index{twin-range!cut-off}}
-To save computation time, slowly varying forces can be calculated
-less often than rapidly varying forces. In {\gromacs}
-such a \normindex{multiple time step} splitting is possible between
-short and long range non-bonded interactions.
-In {\gromacs} versions up to 4.0, an irreversible integration scheme
-was used which is also used by the {\gromos} simulation package:
-every $n$ steps the long range forces are determined and these are
-then also used (without modification) for the next $n-1$ integration steps
-in \eqnref{leapfrogv}. Such an irreversible scheme can result in bad energy
-conservation and, possibly, bad sampling.
-Since version 4.5, a leap-frog version of the reversible Trotter decomposition scheme~\cite{Tuckerman1992a} is used.
-In this integrator the long-range forces are determined every $n$ steps
-and are then integrated into the velocity in \eqnref{leapfrogv} using
-a time step of $\Dt_\mathrm{LR} = n \Dt$:
-\beq
-\ve{v}(t+\hDt) =
-\left\{ \begin{array}{lll} \displaystyle
- \ve{v}(t-\hDt) + \frac{1}{m}\left[\ve{F}_\mathrm{SR}(t) + n \ve{F}_\mathrm{LR}(t)\right] \Dt &,& \mathrm{step} ~\%~ n = 0 \\ \noalign{\medskip} \displaystyle
- \ve{v}(t-\hDt) + \frac{1}{m}\ve{F}_\mathrm{SR}(t)\Dt &,& \mathrm{step} ~\%~ n \neq 0 \\
-\end{array} \right.
-\eeq
-
-The parameter $n$ is equal to the neighbor list update frequency. In
-4.5, the velocity Verlet version of multiple time-stepping is not yet
-fully implemented.
-
+\subsection{Multiple time stepping}
Several other simulation packages uses multiple time stepping for
bonds and/or the PME mesh forces. In {\gromacs} we have not implemented
this (yet), since we use a different philosophy. Bonds can be constrained
which brings the shortest time step up to
PME mesh update frequency of a multiple time stepping scheme.
-As an example we show the energy conservation for integrating
-the equations of motion for SPC/E water at 300 K. To avoid cut-off
-effects, reaction-field electrostatics with $\epsilon_{RF}=\infty$ and
-shifted Lennard-Jones interactions are used, both with a buffer region.
-The long-range interactions were evaluated between 1.0 and 1.4 nm.
-In \figref{leapfrog} one can see that for electrostatics the Trotter scheme
-does an order of magnitude better up to $\Dt_{LR}$ = 16 fs.
-The electrostatics depends strongly on the orientation of the water molecules,
-which changes rapidly.
-For Lennard-Jones interactions, the energy drift is linear in $\Dt_{LR}$
-and roughly two orders of magnitude smaller than for the electrostatics.
-Lennard-Jones forces are smaller than Coulomb forces and
-they are mainly affected by translation of water molecules, not rotation.
-
-\begin{figure}
-\centerline{\includegraphics[width=12cm]{plots/drift-all}}
-\caption{Energy drift per degree of freedom in SPC/E water
-with twin-range cut-offs
-for reaction field (left) and Lennard-Jones interaction (right)
-as a function of the long-range time step length for the irreversible
-``\gromos'' scheme and a reversible Trotter scheme.}
-\label{fig:twinrangeener}
-\end{figure}
-
\subsection{Temperature coupling\index{temperature coupling}}
While direct use of molecular dynamics gives rise to the NVE (constant
number, constant volume, constant energy ensemble), most quantities
parameters, for a total of six non-bonded interaction parameters. See
the User Guide for a complete description of these parameters.
-The neighbor searching (NS) can be performed using a single-range, or a twin-range
-approach. Since the former is merely a special case of the latter, we will
-discuss the more general twin-range. In this case, NS is described by two
-radii: {\tt rlist} and max({\tt rcoulomb},{\tt rvdw}).
-Usually one builds the neighbor list every 10 time steps
-or every 20 fs (parameter {\tt nstlist}). In the neighbor list, all interaction
-pairs that fall within {\tt rlist} are stored. Furthermore, the
-interactions between pairs that do not
-fall within {\tt rlist} but do fall within max({\tt rcoulomb},{\tt rvdw})
-are computed during NS. The
-forces and energy are stored separately and added to short-range forces
-at every time step between successive NS. If {\tt rlist} =
-max({\tt rcoulomb},{\tt rvdw}), no forces
-are evaluated during neighbor list generation.
-The \normindex{virial} is calculated from the sum of the short- and
-long-range forces.
-This means that the virial can be slightly asymmetrical at non-NS steps.
-When mdrun is compiled to use mixed precision, the virial is almost always asymmetrical because the
-off-diagonal elements are about as large as each element in the sum.
-In most cases this is not really a problem, since the fluctuations in the
-virial can be 2 orders of magnitude larger than the average.
-
-Except for the plain cut-off,
-all of the interaction functions in \tabref{funcparm}
-require that neighbor searching be done with a larger radius than the $r_c$
+In the group cut-off scheme, all of the interaction functions in \tabref{funcparm}
+require that neighbor searching be done with a radius at least as large as the $r_c$
specified for the functional form, because of the use of charge groups.
The extra radius is typically of the order of 0.25 nm (roughly the
largest distance between two atoms in a charge group plus the distance a
charge group can diffuse within neighbor list updates).
-%If your charge groups are very large it may be interesting to turn off charge
-%groups, by setting the option
-%{\tt bAtomList = yes} in your {\tt grompp.mdp} file.
-%In this case only a small extra radius to account for diffusion needs to be
-%added (0.1 nm). Do not however use this together with the plain cut-off
-%method, as it will generate large artifacts (\secref{cg}).
-%In summary, there are four parameters that describe NS behavior:
-%{\tt nstlist} (update frequency in number of time steps),
-%{\tt bAtomList} (whether or not charge groups are used to generate neighbor list, the default is to use charge groups, so {\tt bAtomList = no}),
-%{\tt rshort} and {\tt rlong} which are the two radii {\rs} and {\rl}
-%described above.
-
\begin{table}[ht]
\centering
\begin{tabular}{|ll|l|}
User-supplied tabulated interactions yes no
Buckingham VdW interactions yes no
rcoulomb != rvdw yes no
-twin-range yes no
+twin-range no no
==================================== ============ =======
Performance
molecule. When :mdp:`nstlist` is larger than one,
:mdp:`nstlist` insertions are performed in a sphere with radius
:mdp:`rtpi` around a the same random location using the same
- neighborlist (and the same long-range energy when :mdp:`rvdw`
- or :mdp:`rcoulomb` > :mdp:`rlist`, which is only allowed for
- single-atom molecules). Since neighborlist construction is
- expensive, one can perform several extra insertions with the
- same list almost for free. The random seed is set with
+ neighborlist. Since neighborlist construction is expensive,
+ one can perform several extra insertions with the same list
+ almost for free. The random seed is set with
:mdp:`ld-seed`. The temperature for the Boltzmann weighting is
set with :mdp:`ref-t`, this should match the temperature of the
simulation of the original trajectory. Dispersion correction is
(100)
number of steps that elapse between calculating the energies, 0 is
- never. This option is only relevant with dynamics. With a
- twin-range cut-off setup :mdp:`nstcalcenergy` should be equal to
- or a multiple of :mdp:`nstlist`. This option affects the
+ never. This option is only relevant with dynamics. This option affects the
performance in parallel simulations, because calculating energies
requires global communication between all processes which can
become a bottleneck at high parallelization.
.. mdp-value:: >0
- Frequency to update the neighbor list (and the long-range
- forces, when using twin-range cut-offs). When this is 0, the
+ Frequency to update the neighbor list. When this is 0, the
neighbor list is made only once. With energy minimization the
neighborlist will be updated for every energy evaluation when
:mdp:`nstlist` is greater than 0. With :mdp:`Verlet` and
Unused.
-.. mdp:: nstcalclr
-
- (-1) \[steps\]
- Controls the period between calculations of long-range forces when
- using the group cut-off scheme.
-
- .. mdp-value:: 1
-
- Calculate the long-range forces every single step. This is
- useful to have separate neighbor lists with buffers for
- electrostatics and Van der Waals interactions, and in particular
- it makes it possible to have the Van der Waals cutoff longer
- than electrostatics (useful *e.g.* with PME). However, there is
- no point in having identical long-range cutoffs for both
- interaction forms and update them every step - then it will be
- slightly faster to put everything in the short-range list.
-
- .. mdp-value:: >1
-
- Calculate the long-range forces every :mdp:`nstcalclr` steps
- and use a multiple-time-step integrator to combine forces. This
- can now be done more frequently than :mdp:`nstlist` since the
- lists are stored, and it might be a good idea *e.g.* for Van der
- Waals interactions that vary slower than electrostatics.
-
- .. mdp-value:: \-1
-
- Calculate long-range forces on steps where neighbor searching is
- performed. While this is the default value, you might want to
- consider updating the long-range forces more frequently.
-
- Note that twin-range force evaluation might be enabled
- automatically by PP-PME load balancing. This is done in order to
- maintain the chosen Van der Waals interaction radius even if the
- load balancing is changing the electrostatics cutoff. If the
- :ref:`mdp` file already specifies twin-range interactions (*e.g.* to
- evaluate Lennard-Jones interactions with a longer cutoff than
- the PME electrostatics every 2-3 steps), the load balancing will
- have also a small effect on Lennard-Jones, since the short-range
- cutoff (inside which forces are evaluated every step) is
- changed.
-
.. mdp:: ns-type
.. mdp-value:: grid
:mdp:`verlet-buffer-tolerance` option and the value of
:mdp:`rlist` is ignored.
-.. mdp:: rlistlong
-
- (-1) \[nm\]
- Cut-off distance for the long-range neighbor list. This parameter
- is only relevant for a twin-range cut-off setup with switched
- potentials. In that case a buffer region is required to account for
- the size of charge groups. In all other cases this parameter is
- automatically set to the longest cut-off distance.
-
Electrostatics
^^^^^^^^^^^^^^
.. mdp-value:: Cut-off
- Twin range cut-offs with neighborlist cut-off :mdp:`rlist` and
- Coulomb cut-off :mdp:`rcoulomb`, where :mdp:`rcoulomb` >=
- :mdp:`rlist`.
+ Plain cut-off with neighborlist radius :mdp:`rlist` and
+ Coulomb cut-off :mdp:`rcoulomb`, where :mdp:`rlist` >=
+ :mdp:`rcoulomb`.
.. mdp-value:: Ewald
.. mdp-value:: Reaction-Field
Reaction field electrostatics with Coulomb cut-off
- :mdp:`rcoulomb`, where :mdp:`rcoulomb` >= :mdp:`rlist`. The
+ :mdp:`rcoulomb`, where :mdp:`rlist` >= :mdp:`rvdw`. The
dielectric constant beyond the cut-off is
:mdp:`epsilon-rf`. The dielectric constant can be set to
infinity by setting :mdp:`epsilon-rf` =0.
.. mdp-value:: Generalized-Reaction-Field
Generalized reaction field with Coulomb cut-off
- :mdp:`rcoulomb`, where :mdp:`rcoulomb` >= :mdp:`rlist`. The
+ :mdp:`rcoulomb`, where :mdp:`rlist` >= :mdp:`rcoulomb`. The
dielectric constant beyond the cut-off is
:mdp:`epsilon-rf`. The ionic strength is computed from the
number of charged (*i.e.* with non zero charge) charge
dd->bInterCGcons = inter_charge_group_constraints(mtop);
dd->bInterCGsettles = inter_charge_group_settles(mtop);
- if (ir->rlistlong == 0)
+ if (ir->rlist == 0)
{
/* Set the cut-off to some very large value,
* so we don't need if statements everywhere in the code.
}
else
{
- comm->cutoff = ir->rlistlong;
+ comm->cutoff = ir->rlist;
}
comm->cutoff_mbody = 0;
* a cell size, DD cells should be at least the size of the list buffer.
*/
comm->cellsize_limit = std::max(comm->cellsize_limit,
- ir->rlistlong - std::max(ir->rvdw, ir->rcoulomb));
+ ir->rlist - std::max(ir->rvdw, ir->rcoulomb));
if (comm->bInterCGBondeds)
{
copy_ivec(fr->ns->grid->n, ncells_old);
grid_first(fplog, fr->ns->grid, dd, &ddbox,
state_local->box, cell_ns_x0, cell_ns_x1,
- fr->rlistlong, grid_density);
+ fr->rlist, grid_density);
break;
case ecutsVERLET:
nbnxn_get_ncells(fr->nbv->nbs, &ncells_old[XX], &ncells_old[YY]);
}
/* Set the number of atoms required for the force calculation.
- * Forces need to be constrained when using a twin-range setup
- * or with energy minimization. For simple simulations we could
- * avoid some allocation, zeroing and copying, but this is
- * probably not worth the complications ande checking.
+ * Forces need to be constrained when doing energy
+ * minimization. For simple simulations we could avoid some
+ * allocation, zeroing and copying, but this is probably not worth
+ * the complications and checking.
*/
forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
struct pme_setup_t {
real rcut_coulomb; /**< Coulomb cut-off */
real rlist; /**< pair-list cut-off */
- real rlistlong; /**< LR pair-list cut-off */
- int nstcalclr; /**< frequency of evaluating long-range forces for group scheme */
real spacing; /**< (largest) PME grid spacing */
ivec grid; /**< the PME grid dimensions */
real grid_efficiency; /**< ineffiency factor for non-uniform grids <= 1 */
real cut_spacing; /**< the minimum cutoff / PME grid spacing ratio */
real rcut_vdw; /**< Vdw cutoff (does not change) */
real rcut_coulomb_start; /**< Initial electrostatics cutoff */
- int nstcalclr_start; /**< Initial electrostatics cutoff */
real rbuf_coulomb; /**< the pairlist buffer size */
real rbuf_vdw; /**< the pairlist buffer size */
matrix box_start; /**< the initial simulation box */
}
else
{
- if (ic->rcoulomb > ic->rlist)
- {
- pme_lb->rbuf_coulomb = ic->rlistlong - ic->rcoulomb;
- }
- else
- {
- pme_lb->rbuf_coulomb = ic->rlist - ic->rcoulomb;
- }
- if (ic->rvdw > ic->rlist)
- {
- pme_lb->rbuf_vdw = ic->rlistlong - ic->rvdw;
- }
- else
- {
- pme_lb->rbuf_vdw = ic->rlist - ic->rvdw;
- }
+ pme_lb->rbuf_coulomb = ic->rlist - ic->rcoulomb;
+ pme_lb->rbuf_vdw = ic->rlist - ic->rvdw;
}
copy_mat(box, pme_lb->box_start);
pme_lb->rcut_vdw = ic->rvdw;
pme_lb->rcut_coulomb_start = ir->rcoulomb;
- pme_lb->nstcalclr_start = ir->nstcalclr;
pme_lb->cur = 0;
pme_lb->setup[0].rcut_coulomb = ic->rcoulomb;
pme_lb->setup[0].rlist = ic->rlist;
- pme_lb->setup[0].rlistlong = ic->rlistlong;
- pme_lb->setup[0].nstcalclr = ir->nstcalclr;
pme_lb->setup[0].grid[XX] = ir->nkx;
pme_lb->setup[0].grid[YY] = ir->nky;
pme_lb->setup[0].grid[ZZ] = ir->nkz;
if (pme_lb->cutoff_scheme == ecutsVERLET)
{
set->rlist = set->rcut_coulomb + pme_lb->rbuf_coulomb;
- /* We dont use LR lists with Verlet, but this avoids if-statements in further checks */
- set->rlistlong = set->rlist;
}
else
{
tmpr_coulomb = set->rcut_coulomb + pme_lb->rbuf_coulomb;
tmpr_vdw = pme_lb->rcut_vdw + pme_lb->rbuf_vdw;
set->rlist = std::min(tmpr_coulomb, tmpr_vdw);
- set->rlistlong = std::max(tmpr_coulomb, tmpr_vdw);
-
- /* Set the long-range update frequency */
- if (set->rlist == set->rlistlong)
- {
- /* No long-range interactions if the short-/long-range cutoffs are identical */
- set->nstcalclr = 0;
- }
- else if (pme_lb->nstcalclr_start == 0 || pme_lb->nstcalclr_start == 1)
- {
- /* We were not doing long-range before, but now we are since rlist!=rlistlong */
- set->nstcalclr = 1;
- }
- else
- {
- /* We were already doing long-range interactions from the start */
- if (pme_lb->rcut_vdw > pme_lb->rcut_coulomb_start)
- {
- /* We were originally doing long-range VdW-only interactions.
- * If rvdw is still longer than rcoulomb we keep the original nstcalclr,
- * but if the coulomb cutoff has become longer we should update the long-range
- * part every step.
- */
- set->nstcalclr = (tmpr_vdw > tmpr_coulomb) ? pme_lb->nstcalclr_start : 1;
- }
- else
- {
- /* We were not doing any long-range interaction from the start,
- * since it is not possible to do twin-range coulomb for the PME interaction.
- */
- set->nstcalclr = 1;
- }
- }
}
set->spacing = sp;
set = &pme_lb->setup[pme_lb->cur];
set->count++;
- rtab = ir->rlistlong + ir->tabext;
+ rtab = ir->rlist + ir->tabext;
if (set->count % 2 == 1)
{
* better overal performance can be obtained with a slightly
* shorter cut-off and better DD load balancing.
*/
- set_dd_dlb_max_cutoff(cr, pme_lb->setup[pme_lb->fastest].rlistlong);
+ set_dd_dlb_max_cutoff(cr, pme_lb->setup[pme_lb->fastest].rlist);
}
}
cycles_fast = pme_lb->setup[pme_lb->fastest].cycles;
if (OK && ir->ePBC != epbcNONE)
{
- OK = (sqr(pme_lb->setup[pme_lb->cur+1].rlistlong)
+ OK = (sqr(pme_lb->setup[pme_lb->cur+1].rlist)
<= max_cutoff2(ir->ePBC, state->box));
if (!OK)
{
if (DOMAINDECOMP(cr))
{
OK = change_dd_cutoff(cr, state, ir,
- pme_lb->setup[pme_lb->cur].rlistlong);
+ pme_lb->setup[pme_lb->cur].rlist);
if (!OK)
{
/* Failed: do not use this setup */
if (DOMAINDECOMP(cr) && pme_lb->stage > 0)
{
- OK = change_dd_cutoff(cr, state, ir, pme_lb->setup[pme_lb->cur].rlistlong);
+ OK = change_dd_cutoff(cr, state, ir, pme_lb->setup[pme_lb->cur].rlist);
if (!OK)
{
/* For some reason the chosen cut-off is incompatible with DD.
ic->rcoulomb = set->rcut_coulomb;
ic->rlist = set->rlist;
- ic->rlistlong = set->rlistlong;
- ir->nstcalclr = set->nstcalclr;
ic->ewaldcoeff_q = set->ewaldcoeff_q;
/* TODO: centralize the code that sets the potentials shifts */
if (ic->coulomb_modifier == eintmodPOTSHIFT)
* This also ensures that we won't disable the currently
* optimal setting during a second round of PME balancing.
*/
- set_dd_dlb_max_cutoff(cr, fr->ic->rlistlong);
+ set_dd_dlb_max_cutoff(cr, fr->ic->rlist);
}
}
fr->ewaldcoeff_q = fr->ic->ewaldcoeff_q;
fr->ewaldcoeff_lj = fr->ic->ewaldcoeff_lj;
fr->rlist = fr->ic->rlist;
- fr->rlistlong = fr->ic->rlistlong;
fr->rcoulomb = fr->ic->rcoulomb;
fr->rvdw = fr->ic->rvdw;
return setup->grid[XX]*setup->grid[YY]*setup->grid[ZZ];
}
-/*! \brief Retuern the largest short-range list cut-off radius */
-static real pme_loadbal_rlist(const pme_setup_t *setup)
-{
- /* With the group cut-off scheme we can have twin-range either
- * for Coulomb or for VdW, so we need a check here.
- * With the Verlet cut-off scheme rlist=rlistlong.
- */
- if (setup->rcut_coulomb > setup->rlist)
- {
- return setup->rlistlong;
- }
- else
- {
- return setup->rlist;
- }
-}
-
/*! \brief Print one load-balancing setting */
static void print_pme_loadbal_setting(FILE *fplog,
const char *name,
fprintf(fplog,
" %-7s %6.3f nm %6.3f nm %3d %3d %3d %5.3f nm %5.3f nm\n",
name,
- setup->rcut_coulomb, pme_loadbal_rlist(setup),
+ setup->rcut_coulomb, setup->rlist,
setup->grid[XX], setup->grid[YY], setup->grid[ZZ],
setup->spacing, 1/setup->ewaldcoeff_q);
}
double pp_ratio, grid_ratio;
real pp_ratio_temporary;
- pp_ratio_temporary = pme_loadbal_rlist(&pme_lb->setup[pme_lb->cur])/pme_loadbal_rlist(&pme_lb->setup[0]);
+ pp_ratio_temporary = pme_lb->setup[pme_lb->cur].rlist / pme_lb->setup[0].rlist;
pp_ratio = std::pow(static_cast<double>(pp_ratio_temporary), 3.0);
grid_ratio = pme_grid_points(&pme_lb->setup[pme_lb->cur])/
(double)pme_grid_points(&pme_lb->setup[0]);
tpxv_CompElPolyatomicIonsAndMultipleIonTypes, /**< CompEl now can handle polyatomic ions and more than two types of ions */
tpxv_RemoveAdress, /**< removed support for AdResS */
tpxv_PullCoordNGroup, /**< add ngroup to pull coord */
+ tpxv_RemoveTwinRange, /**< removed support for twin-range interactions */
tpxv_Count /**< the total number of tpxv versions */
};
{ 72, F_NPSOLVATION },
{ 41, F_LJC14_Q },
{ 41, F_LJC_PAIRS_NB },
- { 32, F_BHAM_LR },
+ { 32, F_BHAM_LR_NOLONGERUSED },
{ 32, F_RF_EXCL },
{ 32, F_COUL_RECIP },
{ 93, F_LJ_RECIP },
ir->verletbuf_tol = 0;
}
gmx_fio_do_real(fio, ir->rlist);
- if (file_version >= 67)
+ if (file_version >= 67 && file_version < tpxv_RemoveTwinRange)
{
- gmx_fio_do_real(fio, ir->rlistlong);
+ if (bRead)
+ {
+ // Reading such a file version could invoke the twin-range
+ // scheme, about which mdrun should give a fatal error.
+ real dummy_rlistlong = -1;
+ gmx_fio_do_real(fio, dummy_rlistlong);
+
+ if (ir->rlist > 0 && (dummy_rlistlong == 0 || dummy_rlistlong > ir->rlist))
+ {
+ // Get mdrun to issue an error (regardless of
+ // ir->cutoff_scheme).
+ ir->useTwinRange = true;
+ }
+ else
+ {
+ // grompp used to set rlistlong actively. Users were
+ // probably also confused and set rlistlong == rlist.
+ // However, in all remaining cases, it is safe to let
+ // mdrun proceed normally.
+ ir->useTwinRange = false;
+ }
+ }
}
- if (file_version >= 82 && file_version != 90)
+ else
{
- gmx_fio_do_int(fio, ir->nstcalclr);
+ // No need to read or write anything
+ ir->useTwinRange = false;
}
- else
+ if (file_version >= 82 && file_version != 90)
{
- /* Calculate at NS steps */
- ir->nstcalclr = ir->nstlist;
+ // Multiple time-stepping is no longer enabled, but the old
+ // support required the twin-range scheme, for which mdrun
+ // already emits a fatal error.
+ int dummy_nstcalclr = -1;
+ gmx_fio_do_int(fio, dummy_nstcalclr);
}
gmx_fio_do_int(fio, ir->coulombtype);
if (file_version < 32 && ir->coulombtype == eelRF)
}
gmx_fio_do_real(fio, ir->rvdw_switch);
gmx_fio_do_real(fio, ir->rvdw);
- if (file_version < 67)
- {
- ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
- }
gmx_fio_do_int(fio, ir->eDispCorr);
gmx_fio_do_real(fio, ir->epsilon_r);
if (file_version >= 37)
PS("periodic-molecules", EBOOL(ir->bPeriodicMols));
PR("verlet-buffer-tolerance", ir->verletbuf_tol);
PR("rlist", ir->rlist);
- PR("rlistlong", ir->rlistlong);
- PR("nstcalclr", ir->nstcalclr);
/* Options for electrostatics and VdW */
PS("coulombtype", EELTYPE(ir->coulombtype));
static gmx_bool bMeanEmtx = TRUE;
static int skip = 0, nlevels = 20;
static real cutmax = 1e20, cutmin = -1e20, reftemp = 300.0;
- static gmx_bool bCoulSR = TRUE, bCoulLR = FALSE, bCoul14 = FALSE;
- static gmx_bool bLJSR = TRUE, bLJLR = FALSE, bLJ14 = FALSE, bBhamSR = FALSE, bBhamLR = FALSE,
+ static gmx_bool bCoulSR = TRUE, bCoul14 = FALSE;
+ static gmx_bool bLJSR = TRUE, bLJ14 = FALSE, bBhamSR = FALSE,
bFree = TRUE;
t_pargs pa[] = {
{ "-sum", FALSE, etBOOL, {&bSum},
{ "-max", FALSE, etREAL, {&cutmax}, "max value for energies"},
{ "-min", FALSE, etREAL, {&cutmin}, "min value for energies"},
{ "-coulsr", FALSE, etBOOL, {&bCoulSR}, "extract Coulomb SR energies"},
- { "-coullr", FALSE, etBOOL, {&bCoulLR}, "extract Coulomb LR energies"},
{ "-coul14", FALSE, etBOOL, {&bCoul14}, "extract Coulomb 1-4 energies"},
{ "-ljsr", FALSE, etBOOL, {&bLJSR}, "extract Lennard-Jones SR energies"},
- { "-ljlr", FALSE, etBOOL, {&bLJLR}, "extract Lennard-Jones LR energies"},
{ "-lj14", FALSE, etBOOL, {&bLJ14}, "extract Lennard-Jones 1-4 energies"},
{ "-bhamsr", FALSE, etBOOL, {&bBhamSR}, "extract Buckingham SR energies"},
- { "-bhamlr", FALSE, etBOOL, {&bBhamLR}, "extract Buckingham LR energies"},
{ "-free", FALSE, etBOOL, {&bFree}, "calculate free energy"},
{ "-temp", FALSE, etREAL, {&reftemp},
"reference temperature for free energy calculation"}
egrp_use[egCOULSR] = bCoulSR;
egrp_use[egLJSR] = bLJSR;
egrp_use[egBHAMSR] = bBhamSR;
- egrp_use[egCOULLR] = bCoulLR;
- egrp_use[egLJLR] = bLJLR;
- egrp_use[egBHAMLR] = bBhamLR;
egrp_use[egCOUL14] = bCoul14;
egrp_use[egLJ14] = bLJ14;
egrp_use[egTotal] = TRUE;
typedef struct
{
- int nr_inputfiles; /* The number of tpr and mdp input files */
- gmx_int64_t orig_sim_steps; /* Number of steps to be done in the real simulation */
- gmx_int64_t orig_init_step; /* Init step for the real simulation */
- real *rcoulomb; /* The coulomb radii [0...nr_inputfiles] */
- real *rvdw; /* The vdW radii */
- real *rlist; /* Neighbourlist cutoff radius */
- real *rlistlong;
+ int nr_inputfiles; /* The number of tpr and mdp input files */
+ gmx_int64_t orig_sim_steps; /* Number of steps to be done in the real simulation */
+ gmx_int64_t orig_init_step; /* Init step for the real simulation */
+ real *rcoulomb; /* The coulomb radii [0...nr_inputfiles] */
+ real *rvdw; /* The vdW radii */
+ real *rlist; /* Neighbourlist cutoff radius */
int *nkx, *nky, *nkz;
real *fsx, *fsy, *fsz; /* Fourierspacing in x,y,z dimension */
} t_inputinfo;
{
fprintf(fp, " rlist : %f nm\n", ir->rlist);
}
- if (ir->rlistlong != max_cutoff(ir->rvdw, ir->rcoulomb))
- {
- fprintf(fp, " rlistlong : %f nm\n", ir->rlistlong);
- }
/* Print a descriptive line about the tpr settings tested */
fprintf(fp, "\nWill try these real/reciprocal workload settings:\n");
{
fprintf(fp, " rlist");
}
- if (ir->rlistlong != max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb)) )
- {
- fprintf(fp, " rlistlong");
- }
fprintf(fp, " tpr file\n");
/* Loop to create the requested number of tpr input files */
/* Adjust other radii since various conditions need to be fulfilled */
if (eelPME == ir->coulombtype)
{
- /* plain PME, rcoulomb must be equal to rlist */
+ /* plain PME, rcoulomb must be equal to rlist TODO only in the group scheme? */
ir->rlist = ir->rcoulomb;
}
else
ir->rvdw = std::max(info->rvdw[0], ir->rlist);
}
}
-
- ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
-
} /* end of "if (j != 0)" */
/* for j==0: Save the original settings
info->nky[j] = ir->nky;
info->nkz[j] = ir->nkz;
info->rlist[j] = ir->rlist;
- info->rlistlong[j] = ir->rlistlong;
info->fsx[j] = fac*fourierspacing;
info->fsy[j] = fac*fourierspacing;
info->fsz[j] = fac*fourierspacing;
{
fprintf(fp, "%10f", ir->rlist);
}
- if (info->rlistlong[0] != max_cutoff(info->rlist[0], max_cutoff(info->rvdw[0], info->rcoulomb[0])) )
- {
- fprintf(fp, "%10f", ir->rlistlong);
- }
fprintf(fp, " %-14s\n", fn_bench_tprs[j]);
/* Make it clear to the user that some additional settings were modified */
if (!gmx_within_tol(ir->rvdw, info->rvdw[0], GMX_REAL_EPS)
- || !gmx_within_tol(ir->rlistlong, info->rlistlong[0], GMX_REAL_EPS) )
+ || !gmx_within_tol(ir->rlist, info->rlist[0], GMX_REAL_EPS) )
{
bNote = TRUE;
}
snew(info->rcoulomb, ntprs);
snew(info->rvdw, ntprs);
snew(info->rlist, ntprs);
- snew(info->rlistlong, ntprs);
snew(info->nkx, ntprs);
snew(info->nky, ntprs);
snew(info->nkz, ntprs);
def_bondedz ("LJC_NB", "LJC Pairs NB", 2, 4, 0, eNR_NB14, unimplemented ),
def_nb ("LJ_SR", "LJ (SR)", 2, 2 ),
def_nb ("BHAM", "Buck.ham (SR)", 2, 3 ),
- def_nofc ("LJ_LR", "LJ (LR)" ),
- def_nofc ("BHAM_LR", "Buck.ham (LR)" ),
+ def_nofc ("LJ_LR", "LJ (unused)" ),
+ def_nofc ("BHAM_LR", "B.ham (unused)" ),
def_nofc ("DISPCORR", "Disper. corr." ),
def_nofc ("COUL_SR", "Coulomb (SR)" ),
- def_nofc ("COUL_LR", "Coulomb (LR)" ),
+ def_nofc ("COUL_LR", "Coul (unused)" ),
def_nofc ("RF_EXCL", "RF excl." ),
def_nofc ("COUL_RECIP", "Coul. recip." ),
def_nofc ("LJ_RECIP", "LJ recip." ),
}
void do_nonbonded(t_forcerec *fr,
- rvec x[], rvec f_shortrange[], rvec f_longrange[], t_mdatoms *mdatoms, t_blocka *excl,
+ rvec x[], rvec f_shortrange[], t_mdatoms *mdatoms, t_blocka *excl,
gmx_grppairener_t *grppener,
t_nrnb *nrnb, real *lambda, real *dvdl,
int nls, int eNL, int flags)
{
t_nblist * nlist;
- int n, n0, n1, i, i0, i1, range;
+ int n, n0, n1, i, i0, i1;
t_nblists * nblists;
nb_kernel_data_t kernel_data;
nb_kernel_t * kernelptr = NULL;
kernel_data.table_vdw = nblists->table_vdw;
kernel_data.table_elec_vdw = nblists->table_elec_vdw;
- for (range = 0; range < 2; range++)
{
- /* Are we doing short/long-range? */
- if (range == 0)
{
/* Short-range */
if (!(flags & GMX_NONBONDED_DO_SR))
nlist = nblists->nlist_sr;
f = f_shortrange;
}
- else
- {
- /* Long-range */
- if (!(flags & GMX_NONBONDED_DO_LR))
- {
- continue;
- }
- kernel_data.energygrp_elec = grppener->ener[egCOULLR];
- kernel_data.energygrp_vdw = grppener->ener[fr->bBHAM ? egBHAMLR : egLJLR];
- kernel_data.energygrp_polarization = grppener->ener[egGB];
- nlist = nblists->nlist_lr;
- f = f_longrange;
- }
for (i = i0; (i < i1); i++)
{
-#define GMX_NONBONDED_DO_LR (1<<0)
#define GMX_NONBONDED_DO_FORCE (1<<1)
#define GMX_NONBONDED_DO_SHIFTFORCE (1<<2)
#define GMX_NONBONDED_DO_FOREIGNLAMBDA (1<<3)
void
do_nonbonded(t_forcerec *fr,
- rvec x[], rvec f_shortrange[], rvec f_longrange[], t_mdatoms *md, t_blocka *excl,
+ rvec x[], rvec f_shortrange[], t_mdatoms *md, t_blocka *excl,
gmx_grppairener_t *grppener,
t_nrnb *nrnb, real *lambda, real dvdlambda[],
int nls, int eNL, int flags);
nshells++;
}
}
- if (inputrecTwinRange(ir) && (nshells > 0))
- {
- snprintf(warn_buf, STRLEN,
- "The combination of using shells and a twin-range cut-off is not supported");
- warning_error(wi, warn_buf);
- }
if ((nshells > 0) && (ir->nstcalcenergy != 1))
{
set_warning_line(wi, "unknown", -1);
printf("Calculated rlist for %dx%d atom pair-list as %.3f nm, buffer size %.3f nm\n",
1, 1, rlist_1x1, rlist_1x1-std::max(ir->rvdw, ir->rcoulomb));
- ir->rlistlong = ir->rlist;
printf("Set rlist, assuming %dx%d atom pair-list, to %.3f nm, buffer size %.3f nm\n",
ls.cluster_size_i, ls.cluster_size_j,
ir->rlist, ir->rlist-std::max(ir->rvdw, ir->rcoulomb));
printf("Note that mdrun will redetermine rlist based on the actual pair-list setup\n");
- if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
+ if (sqr(ir->rlist) >= max_cutoff2(ir->ePBC, box))
{
- gmx_fatal(FARGS, "The pair-list cut-off (%g nm) is longer than half the shortest box vector or longer than the smallest box diagonal element (%g nm). Increase the box size or decrease nstlist or increase verlet-buffer-tolerance.", ir->rlistlong, std::sqrt(max_cutoff2(ir->ePBC, box)));
+ gmx_fatal(FARGS, "The pair-list cut-off (%g nm) is longer than half the shortest box vector or longer than the smallest box diagonal element (%g nm). Increase the box size or decrease nstlist or increase verlet-buffer-tolerance.", ir->rlist, std::sqrt(max_cutoff2(ir->ePBC, box)));
}
}
"release when all interaction forms are supported for the verlet scheme. The verlet "
"scheme already scales better, and it is compatible with GPUs and other accelerators.");
- /* BASIC CUT-OFF STUFF */
- if (ir->rlist == 0 ||
- !((ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > ir->rlist) ||
- (ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > ir->rlist)))
+ if (ir->rlist > 0 && ir->rlist < ir->rcoulomb)
{
- /* No switched potential and/or no twin-range:
- * we can set the long-range cut-off to the maximum of the other cut-offs.
- */
- ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
+ gmx_fatal(FARGS, "rcoulomb must not be greater than rlist (twin-range schemes are not supported)");
}
- else if (ir->rlistlong < 0)
+ if (ir->rlist > 0 && ir->rlist < ir->rvdw)
{
- ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
- sprintf(warn_buf, "rlistlong was not set, setting it to %g (no buffer)",
- ir->rlistlong);
- warning(wi, warn_buf);
+ gmx_fatal(FARGS, "rvdw must not be greater than rlist (twin-range schemes are not supported)");
}
- if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
+
+ if (ir->rlist == 0 && ir->ePBC != epbcNONE)
{
warning_error(wi, "Can not have an infinite cut-off with PBC");
}
- if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
- {
- warning_error(wi, "rlistlong can not be shorter than rlist");
- }
- if (inputrecTwinRange(ir) && ir->nstlist == 0)
- {
- warning_error(wi, "Can not have nstlist == 0 with twin-range interactions");
- }
- if (inputrecTwinRange(ir) && EI_VV(ir->eI))
- {
- sprintf(warn_buf, "Twin-range interactions are not supported with integrator %s.", ei_names[ir->eI]);
- warning_error(wi, warn_buf);
- }
- }
-
- if (ir->rlistlong == ir->rlist)
- {
- ir->nstcalclr = 0;
- }
- else if (ir->rlistlong > ir->rlist && ir->nstcalclr == 0)
- {
- warning_error(wi, "With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
}
if (ir->cutoff_scheme == ecutsVERLET)
}
}
}
-
- /* No twin-range calculations with Verlet lists */
- ir->rlistlong = ir->rlist;
- }
-
- if (ir->nstcalclr == -1)
- {
- /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
- ir->nstcalclr = ir->nstlist;
- }
- else if (ir->nstcalclr > 0)
- {
- if (ir->nstlist > 0 && (ir->nstlist % ir->nstcalclr != 0))
- {
- warning_error(wi, "nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
- }
- }
- else if (ir->nstcalclr < -1)
- {
- warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
- }
-
- if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rlist && ir->nstcalclr > 1)
- {
- warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
}
/* GENERAL INTEGRATOR STUFF */
ir->nstpcouple = ir_optimal_nstpcouple(ir);
}
}
- if (inputrecTwinRange(ir))
- {
- check_nst("nstcalclr", ir->nstcalclr,
- "nstcalcenergy", &ir->nstcalcenergy, wi);
- if (ir->epc != epcNO)
- {
- check_nst("nstlist", ir->nstlist,
- "nstpcouple", &ir->nstpcouple, wi);
- }
- }
if (ir->nstcalcenergy > 0)
{
CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
}
}
-
- if (inputrecTwinRange(ir))
- {
- sprintf(err_buf, "nstdhdl must be divisible by nstcalclr");
- CHECK(ir->fepvals->nstdhdl > 0 &&
- ir->fepvals->nstdhdl % ir->nstcalclr != 0);
-
- if (ir->efep == efepEXPANDED)
- {
- sprintf(err_buf, "nstexpanded must be divisible by nstcalclr");
- CHECK(ir->expandedvals->nstexpanded % ir->nstcalclr != 0);
- }
- }
}
if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
warning_note(wi, "Tumbling and or flying ice-cubes: We are not removing rotation around center of mass in a non-periodic system. You should probably set comm_mode = ANGULAR.");
}
- sprintf(err_buf, "Twin-range neighbour searching (NS) with simple NS"
- " algorithm not implemented");
- CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
- && (ir->ns_type == ensSIMPLE));
-
/* TEMPERATURE COUPLING */
if (ir->etc == etcYES)
{
if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
{
sprintf(err_buf,
- "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
- "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
+ "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist.\n"
+ "For optimal energy conservation,consider using\n"
"a potential modifier.", eel_names[ir->coulombtype]);
- if (ir->nstcalclr == 1)
- {
- CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
- }
- else
- {
- CHECK(ir->rcoulomb != ir->rlist);
- }
+ CHECK(ir->rcoulomb != ir->rlist);
}
}
}
"between neighborsearch steps");
}
- if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlistlong <= ir->rcoulomb)
+ if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlist <= ir->rcoulomb)
{
- sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
- inputrecTwinRange(ir) ? "rlistlong" : "rlist");
+ sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rcoulomb.");
warning_note(wi, warn_buf);
}
- if (ir_vdw_switched(ir) && (ir->rlistlong <= ir->rvdw))
+ if (ir_vdw_switched(ir) && (ir->rlist <= ir->rvdw))
{
- sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
- inputrecTwinRange(ir) ? "rlistlong" : "rlist");
+ sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rvdw.");
warning_note(wi, warn_buf);
}
}
REM_TYPE("adress_tf_grp_names");
REM_TYPE("adress_cg_grp_names");
REM_TYPE("adress_do_hybridpairs");
+ REM_TYPE("rlistlong");
+ REM_TYPE("nstcalclr");
/* replace the following commands with the clearer new versions*/
REPL_TYPE("unconstrained-start", "continuation");
CTYPE ("nblist cut-off");
RTYPE ("rlist", ir->rlist, 1.0);
CTYPE ("long-range cut-off for switched potentials");
- RTYPE ("rlistlong", ir->rlistlong, -1);
- ITYPE ("nstcalclr", ir->nstcalclr, -1);
/* Electrostatics */
CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
warninp_t wi)
{
real min_size;
- gmx_bool bTWIN;
char warn_buf[STRLEN];
const char *ptr;
ir->shake_tol);
warning_error(wi, warn_buf);
}
-
- if (inputrecTwinRange(ir) && ir->nstlist > 1)
- {
- sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
- if (ir->epc == epcNO)
- {
- warning(wi, warn_buf);
- }
- else
- {
- warning_error(wi, warn_buf);
- }
- }
}
if ( (ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
{
warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
}
- bTWIN = (ir->rlistlong > ir->rlist);
if (ir->ns_type == ensGRID)
{
- if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
+ if (sqr(ir->rlist) >= max_cutoff2(ir->ePBC, box))
{
- sprintf(warn_buf, "ERROR: The cut-off length is longer than half the shortest box vector or longer than the smallest box diagonal element. Increase the box size or decrease %s.\n",
- bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
+ sprintf(warn_buf, "ERROR: The cut-off length is longer than half the shortest box vector or longer than the smallest box diagonal element. Increase the box size or decrease rlist.\n");
warning_error(wi, warn_buf);
}
}
else
{
min_size = std::min(box[XX][XX], std::min(box[YY][YY], box[ZZ][ZZ]));
- if (2*ir->rlistlong >= min_size)
+ if (2*ir->rlist >= min_size)
{
sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
warning_error(wi, warn_buf);
* not be zero at the cut-off.
*/
if (ir_vdw_is_zero_at_cutoff(ir) &&
- rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
+ rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
{
sprintf(warn_buf, "The sum of the two largest charge group "
- "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
+ "radii (%f) is larger than rlist (%f) - rvdw (%f).\n"
"With exact cut-offs, better performance can be "
"obtained with cutoff-scheme = %s, because it "
"does not use charge groups at all.",
rvdw1+rvdw2,
- ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
- ir->rlistlong, ir->rvdw,
+ ir->rlist, ir->rvdw,
ecutscheme_names[ecutsVERLET]);
if (ir_NVE(ir))
{
}
}
if (ir_coulomb_is_zero_at_cutoff(ir) &&
- rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
+ rcoul1 + rcoul2 > ir->rlist - ir->rcoulomb)
{
- sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
+ sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rcoulomb (%f).\n"
"With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
rcoul1+rcoul2,
- ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
- ir->rlistlong, ir->rcoulomb,
+ ir->rlist, ir->rcoulomb,
ecutscheme_names[ecutsVERLET]);
if (ir_NVE(ir))
{
ene->E_dihe = (float) enerd->term[F_PDIHS ];
ene->E_impr = (float) enerd->term[F_IDIHS ];
ene->E_vdw = (float) enerd->term[F_LJ ];
- ene->E_coul = (float) (enerd->term[F_COUL_SR] + enerd->term[F_COUL_LR]);
+ ene->E_coul = (float) enerd->term[F_COUL_SR];
}
}
#else
t_mdatoms *md,
t_commrec *cr,
t_nrnb *nrnb,
- gmx_bool bFillGrid,
- gmx_bool bDoLongRangeNS)
+ gmx_bool bFillGrid)
{
int nsearch;
init_neighbor_list(fp, fr, md->homenr);
}
- if (fr->bTwinRange)
- {
- fr->nlr = 0;
- }
-
nsearch = search_neighbours(fp, fr, box, top, groups, cr, nrnb, md,
- bFillGrid, bDoLongRangeNS);
+ bFillGrid);
if (debug)
{
fprintf(debug, "nsearch = %d\n", nsearch);
t_mdatoms *md,
rvec x[], history_t *hist,
rvec f[],
- rvec f_longrange[],
gmx_enerdata_t *enerd,
t_fcdata *fcd,
gmx_localtop_t *top,
{
donb_flags |= GMX_NONBONDED_DO_POTENTIAL;
}
- if (flags & GMX_FORCE_DO_LR)
- {
- donb_flags |= GMX_NONBONDED_DO_LR;
- }
wallcycle_sub_start(wcycle, ewcsNONBONDED);
- do_nonbonded(fr, x, f, f_longrange, md, excl,
+ do_nonbonded(fr, x, f, md, excl,
&enerd->grpp, nrnb,
lambda, dvdl_nb, -1, -1, donb_flags);
lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]);
}
reset_foreign_enerdata(enerd);
- do_nonbonded(fr, x, f, f_longrange, md, excl,
+ do_nonbonded(fr, x, f, md, excl,
&(enerd->foreign_grpp), nrnb,
lam_i, dvdl_dum, -1, -1,
(donb_flags & ~GMX_NONBONDED_DO_FORCE) | GMX_NONBONDED_DO_FOREIGNLAMBDA);
epot[F_LJ] = sum_v(grpp->nener, grpp->ener[egLJSR]);
epot[F_LJ14] = sum_v(grpp->nener, grpp->ener[egLJ14]);
epot[F_COUL14] = sum_v(grpp->nener, grpp->ener[egCOUL14]);
- epot[F_COUL_LR] = sum_v(grpp->nener, grpp->ener[egCOULLR]);
- epot[F_LJ_LR] = sum_v(grpp->nener, grpp->ener[egLJLR]);
/* We have already added 1-2,1-3, and 1-4 terms to F_GBPOL */
epot[F_GBPOL] += sum_v(grpp->nener, grpp->ener[egGB]);
* and has been added earlier
*/
epot[F_BHAM] = sum_v(grpp->nener, grpp->ener[egBHAMSR]);
- epot[F_BHAM_LR] = sum_v(grpp->nener, grpp->ener[egBHAMLR]);
epot[F_EPOT] = 0;
for (i = 0; (i < F_EPOT); i++)
* For the constraints this is not exact, but we have no other option
* without literally changing the lengths and reevaluating the energies at each step.
* (try to remedy this post 4.6 - MRS)
- * For the non-bonded LR term we assume that the soft-core (if present)
- * no longer affects the energy beyond the short-range cut-off,
- * which is a very good approximation (except for exotic settings).
- * (investigate how to overcome this post 4.6 - MRS)
*/
if (fepvals->separate_dvdl[efptBONDED])
{
}
}
-void reset_enerdata(t_forcerec *fr, gmx_bool bNS,
- gmx_enerdata_t *enerd,
- gmx_bool bMaster)
+void reset_enerdata(gmx_enerdata_t *enerd)
{
- gmx_bool bKeepLR;
int i, j;
- /* First reset all energy components, except for the long range terms
- * on the master at non neighbor search steps, since the long range
- * terms have already been summed at the last neighbor search step.
- */
- bKeepLR = (fr->bTwinRange && !bNS);
+ /* First reset all energy components. */
for (i = 0; (i < egNR); i++)
{
- if (!(bKeepLR && bMaster && (i == egCOULLR || i == egLJLR)))
+ for (j = 0; (j < enerd->grpp.nener); j++)
{
- for (j = 0; (j < enerd->grpp.nener); j++)
- {
- enerd->grpp.ener[i][j] = 0.0;
- }
+ enerd->grpp.ener[i][j] = 0.0;
}
}
for (i = 0; i < efptNR; i++)
{
enerd->term[i] = 0.0;
}
- /* Initialize the dVdlambda term with the long range contribution */
- /* Initialize the dvdl term with the long range contribution */
enerd->term[F_DVDL] = 0.0;
enerd->term[F_DVDL_COUL] = 0.0;
enerd->term[F_DVDL_VDW] = 0.0;
void reset_foreign_enerdata(gmx_enerdata_t *enerd);
/* Resets only the foreign energy data */
-void reset_enerdata(t_forcerec *fr, gmx_bool bNS,
- gmx_enerdata_t *enerd,
- gmx_bool bMaster);
-/* Resets the energy data, if bNS=TRUE also zeros the long-range part */
+void reset_enerdata(gmx_enerdata_t *enerd);
+/* Resets the energy data */
void sum_epot(gmx_grppairener_t *grpp, real *epot);
/* Locally sum the non-bonded potential energy terms */
t_mdatoms *md,
t_commrec *cr,
t_nrnb *nrnb,
- gmx_bool bFillGrid,
- gmx_bool bDoLongRangeNS);
+ gmx_bool bFillGrid);
/* Call the neighborsearcher */
void do_force_lowlevel(t_forcerec *fr,
rvec x[],
history_t *hist,
rvec f_shortrange[],
- rvec f_longrange[],
gmx_enerdata_t *enerd,
t_fcdata *fcd,
gmx_localtop_t *top,
#define GMX_FORCE_DYNAMICBOX (1<<1)
/* Do neighbor searching */
#define GMX_FORCE_NS (1<<2)
-/* Update long-range neighborlists */
-#define GMX_FORCE_LRNS (1<<3)
/* Calculate listed energies/forces (e.g. bonds, restraints, 1-4, FEP non-bonded) */
#define GMX_FORCE_LISTED (1<<4)
-/* Store long-range forces in a separate array */
-#define GMX_FORCE_SEPLRF (1<<5)
/* Calculate non-bonded energies/forces */
#define GMX_FORCE_NONBONDED (1<<6)
/* Calculate forces (not only energies) */
#define GMX_FORCE_ENERGY (1<<9)
/* Calculate dHdl */
#define GMX_FORCE_DHDL (1<<10)
-/* Calculate long-range energies/forces */
-#define GMX_FORCE_DO_LR (1<<11)
/* Normally one want all energy terms and forces */
#define GMX_FORCE_ALLFORCES (GMX_FORCE_LISTED | GMX_FORCE_NONBONDED | GMX_FORCE_FORCES)
#include "nbnxn_gpu_jit_support.h"
const char *egrp_nm[egNR+1] = {
- "Coul-SR", "LJ-SR", "Buck-SR", "Coul-LR", "LJ-LR", "Buck-LR",
+ "Coul-SR", "LJ-SR", "Buck-SR",
"Coul-14", "LJ-14", NULL
};
if (fr->natoms_force_constr > fr->nalloc_force)
{
fr->nalloc_force = over_alloc_dd(fr->natoms_force_constr);
-
- if (fr->bTwinRange)
- {
- srenew(fr->f_twin, fr->nalloc_force);
- }
}
if (fr->bF_NoVirSum)
snew_aligned(ic->tabq_coul_V, 16, 32);
ic->rlist = fr->rlist;
- ic->rlistlong = fr->rlistlong;
/* Lennard-Jones */
ic->vdwtype = fr->vdwtype;
{
gmx_fatal(FARGS, "AdResS simulations are no longer supported");
}
-
+ if (ir->useTwinRange)
+ {
+ gmx_fatal(FARGS, "Twin-range simulations are no longer supported");
+ }
/* Copy the user determined parameters */
fr->userint1 = ir->userint1;
fr->userint2 = ir->userint2;
copy_rvec(ir->posres_com, fr->posres_com);
copy_rvec(ir->posres_comB, fr->posres_comB);
fr->rlist = cutoff_inf(ir->rlist);
- fr->rlistlong = cutoff_inf(ir->rlistlong);
fr->eeltype = ir->coulombtype;
fr->vdwtype = ir->vdwtype;
fr->ljpme_combination_rule = ir->ljpme_combination_rule;
fr->rcoulomb = cutoff_inf(ir->rcoulomb);
fr->rcoulomb_switch = ir->rcoulomb_switch;
- fr->bTwinRange = fr->rlistlong > fr->rlist;
fr->bEwald = (EEL_PME(fr->eeltype) || fr->eeltype == eelEWALD);
fr->reppow = mtop->ffparams.reppow;
* in that case grompp should already have checked that we do not need
* normal tables and we only generate tables for 1-4 interactions.
*/
- rtab = ir->rlistlong + ir->tabext;
+ rtab = ir->rlist + ir->tabext;
if (bMakeTables)
{
pr_real(fp, fr->rcoulomb);
pr_real(fp, fr->fudgeQQ);
pr_bool(fp, fr->bGrid);
- pr_bool(fp, fr->bTwinRange);
/*pr_int(fp,fr->cg0);
pr_int(fp,fr->hcg);*/
for (i = 0; i < fr->nnblists; i++)
return nstglobalcomm;
}
-void check_ir_old_tpx_versions(t_commrec *cr, FILE *fplog,
- t_inputrec *ir, gmx_mtop_t *mtop)
-{
- /* Check required for old tpx files */
- if (inputrecTwinRange(ir) && ir->nstlist > 1 &&
- ir->nstcalcenergy % ir->nstlist != 0)
- {
- md_print_warn(cr, fplog, "Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies\n");
-
- if (gmx_mtop_ftype_count(mtop, F_CONSTR) +
- gmx_mtop_ftype_count(mtop, F_CONSTRNC) > 0 &&
- ir->eConstrAlg == econtSHAKE)
- {
- md_print_warn(cr, fplog, "With twin-range cut-off's and SHAKE the virial and pressure are incorrect\n");
- if (ir->epc != epcNO)
- {
- gmx_fatal(FARGS, "Can not do pressure coupling with twin-range cut-off's and SHAKE");
- }
- }
- check_nst_param(fplog, cr, "nstlist", ir->nstlist,
- "nstcalcenergy", &ir->nstcalcenergy);
- if (ir->epc != epcNO)
- {
- check_nst_param(fplog, cr, "nstlist", ir->nstlist,
- "nstpcouple", &ir->nstpcouple);
- }
- check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
- "nstenergy", &ir->nstenergy);
- check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
- "nstlog", &ir->nstlog);
- if (ir->efep != efepNO)
- {
- check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
- "nstdhdl", &ir->fepvals->nstdhdl);
- }
- }
-
- if (EI_VV(ir->eI) && inputrecTwinRange(ir) && ir->nstlist > 1)
- {
- gmx_fatal(FARGS, "Twin-range multiple time stepping does not work with integrator %s.", ei_names[ir->eI]);
- }
-}
-
void rerun_parallel_comm(t_commrec *cr, t_trxframe *fr,
gmx_bool *bNotLastFrame)
{
{
md->bEner[i] = ir->bQMMM;
}
- else if (i == F_COUL_LR)
- {
- md->bEner[i] = (ir->rcoulomb > ir->rlist);
- }
- else if (i == F_LJ_LR)
- {
- md->bEner[i] = (!bBHAM && ir->rvdw > ir->rlist);
- }
- else if (i == F_BHAM_LR)
- {
- md->bEner[i] = (bBHAM && ir->rvdw > ir->rlist);
- }
else if (i == F_RF_EXCL)
{
md->bEner[i] = (EEL_RF(ir->coulombtype) && ir->cutoff_scheme == ecutsGROUP);
md->bEInd[egCOULSR] = TRUE;
md->bEInd[egLJSR ] = TRUE;
- if (ir->rcoulomb > ir->rlist)
- {
- md->bEInd[egCOULLR] = TRUE;
- }
- if (!bBHAM)
- {
- if (ir->rvdw > ir->rlist)
- {
- md->bEInd[egLJLR] = TRUE;
- }
- }
- else
+ if (bBHAM)
{
md->bEInd[egLJSR] = FALSE;
md->bEInd[egBHAMSR] = TRUE;
- if (ir->rvdw > ir->rlist)
- {
- md->bEInd[egBHAMLR] = TRUE;
- }
}
if (b14)
{
void PrintFreeEnergyInfoToFile(FILE *outfile, t_lambda *fep, t_expanded *expand, t_simtemp *simtemp, df_history_t *dfhist,
int fep_state, int frequency, gmx_int64_t step);
-/* check the version */
-void check_ir_old_tpx_versions(t_commrec *cr, FILE *fplog,
- t_inputrec *ir, gmx_mtop_t *mtop);
-
/* Allocate and initialize node-local state entries. */
void set_state_entries(t_state *state, const t_inputrec *ir);
ems->s.lambda, graph, fr, vsite, mu_tot, t, NULL, NULL, TRUE,
GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES |
GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY |
- (bNS ? GMX_FORCE_NS | GMX_FORCE_DO_LR : 0));
+ (bNS ? GMX_FORCE_NS : 0));
/* Clear the unused shake virial and pressure */
clear_mat(shake_vir);
}
-static void init_nblist(FILE *log, t_nblist *nl_sr, t_nblist *nl_lr,
- int maxsr, int maxlr,
+static void init_nblist(FILE *log, t_nblist *nl_sr,
+ int maxsr,
int ivdw, int ivdwmod,
int ielec, int ielecmod,
int igeometry, int type,
{
t_nblist *nl;
int homenr;
- int i;
- for (i = 0; (i < 2); i++)
{
- nl = (i == 0) ? nl_sr : nl_lr;
- homenr = (i == 0) ? maxsr : maxlr;
+ nl = nl_sr;
+ homenr = maxsr;
if (nl == NULL)
{
- continue;
+ return;
}
}
/* This will also set the simd_padding_width field */
- gmx_nonbonded_set_kernel_pointers( (i == 0) ? log : NULL, nl, bElecAndVdwSwitchDiffers);
+ gmx_nonbonded_set_kernel_pointers(log, nl, bElecAndVdwSwitchDiffers);
/* maxnri is influenced by the number of shifts (maximum is 8)
* and the number of energy groups.
if (debug)
{
- fprintf(debug, "Initiating neighbourlist (ielec=%d, ivdw=%d, type=%d) for %s interactions,\nwith %d SR, %d LR atoms.\n",
- nl->ielec, nl->ivdw, nl->type, gmx_nblist_geometry_names[nl->igeometry], maxsr, maxlr);
+ fprintf(debug, "Initiating neighbourlist (ielec=%d, ivdw=%d, type=%d) for %s interactions,\nwith %d SR atoms.\n",
+ nl->ielec, nl->ivdw, nl->type, gmx_nblist_geometry_names[nl->igeometry], maxsr);
}
}
}
void init_neighbor_list(FILE *log, t_forcerec *fr, int homenr)
{
- /* Make maxlr tunable! (does not seem to be a big difference though)
- * This parameter determines the number of i particles in a long range
- * neighbourlist. Too few means many function calls, too many means
- * cache trashing.
- */
- int maxsr, maxsr_wat, maxlr, maxlr_wat;
+ int maxsr, maxsr_wat;
int ielec, ivdw, ielecmod, ivdwmod, type;
int igeometry_def, igeometry_w, igeometry_ww;
int i;
* The numbers seem very accurate, but they are uncritical.
*/
maxsr_wat = std::min(fr->nWatMol, (homenr+2)/3);
- if (fr->bTwinRange)
- {
- maxlr = 50;
- maxlr_wat = std::min(maxsr_wat, maxlr);
- }
- else
- {
- maxlr = maxlr_wat = 0;
- }
/* Determine the values for ielec/ivdw. */
ielec = fr->nbkernel_elec_interaction;
{
nbl = &(fr->nblists[i]);
- init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ], &nbl->nlist_lr[eNL_VDWQQ],
- maxsr, maxlr, ivdw, ivdwmod, ielec, ielecmod, igeometry_def, type, bElecAndVdwSwitchDiffers);
- init_nblist(log, &nbl->nlist_sr[eNL_VDW], &nbl->nlist_lr[eNL_VDW],
- maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, igeometry_def, type, bElecAndVdwSwitchDiffers);
- init_nblist(log, &nbl->nlist_sr[eNL_QQ], &nbl->nlist_lr[eNL_QQ],
- maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_def, type, bElecAndVdwSwitchDiffers);
- init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_WATER], &nbl->nlist_lr[eNL_VDWQQ_WATER],
- maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_w, type, bElecAndVdwSwitchDiffers);
- init_nblist(log, &nbl->nlist_sr[eNL_QQ_WATER], &nbl->nlist_lr[eNL_QQ_WATER],
- maxsr_wat, maxlr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_w, type, bElecAndVdwSwitchDiffers);
- init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_WATERWATER], &nbl->nlist_lr[eNL_VDWQQ_WATERWATER],
- maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_ww, type, bElecAndVdwSwitchDiffers);
- init_nblist(log, &nbl->nlist_sr[eNL_QQ_WATERWATER], &nbl->nlist_lr[eNL_QQ_WATERWATER],
- maxsr_wat, maxlr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_ww, type, bElecAndVdwSwitchDiffers);
+ init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ],
+ maxsr, ivdw, ivdwmod, ielec, ielecmod, igeometry_def, type, bElecAndVdwSwitchDiffers);
+ init_nblist(log, &nbl->nlist_sr[eNL_VDW],
+ maxsr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, igeometry_def, type, bElecAndVdwSwitchDiffers);
+ init_nblist(log, &nbl->nlist_sr[eNL_QQ],
+ maxsr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_def, type, bElecAndVdwSwitchDiffers);
+ init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_WATER],
+ maxsr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_w, type, bElecAndVdwSwitchDiffers);
+ init_nblist(log, &nbl->nlist_sr[eNL_QQ_WATER],
+ maxsr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_w, type, bElecAndVdwSwitchDiffers);
+ init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_WATERWATER],
+ maxsr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_ww, type, bElecAndVdwSwitchDiffers);
+ init_nblist(log, &nbl->nlist_sr[eNL_QQ_WATERWATER],
+ maxsr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_ww, type, bElecAndVdwSwitchDiffers);
/* Did we get the solvent loops so we can use optimized water kernels? */
if (nbl->nlist_sr[eNL_VDWQQ_WATER].kernelptr_vf == NULL
if (fr->efep != efepNO)
{
- init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_FREE], &nbl->nlist_lr[eNL_VDWQQ_FREE],
- maxsr, maxlr, ivdw, ivdwmod, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
- init_nblist(log, &nbl->nlist_sr[eNL_VDW_FREE], &nbl->nlist_lr[eNL_VDW_FREE],
- maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
- init_nblist(log, &nbl->nlist_sr[eNL_QQ_FREE], &nbl->nlist_lr[eNL_QQ_FREE],
- maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
+ init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_FREE],
+ maxsr, ivdw, ivdwmod, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
+ init_nblist(log, &nbl->nlist_sr[eNL_VDW_FREE],
+ maxsr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
+ init_nblist(log, &nbl->nlist_sr[eNL_QQ_FREE],
+ maxsr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
}
}
/* QMMM MM list */
{
snew(fr->QMMMlist, 1);
}
- init_nblist(log, fr->QMMMlist, NULL,
- maxsr, maxlr, 0, 0, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_STANDARD, bElecAndVdwSwitchDiffers);
+ init_nblist(log, fr->QMMMlist,
+ maxsr, 0, 0, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_STANDARD, bElecAndVdwSwitchDiffers);
}
if (log != NULL)
}
}
-static void reset_neighbor_lists(t_forcerec *fr, gmx_bool bResetSR, gmx_bool bResetLR)
+static void reset_neighbor_lists(t_forcerec *fr)
{
int n, i;
{
for (i = 0; i < eNL_NR; i++)
{
- if (bResetSR)
- {
- reset_nblist( &(fr->nblists[n].nlist_sr[i]) );
- }
- if (bResetLR)
- {
- reset_nblist( &(fr->nblists[n].nlist_lr[i]) );
- }
+ reset_nblist( &(fr->nblists[n].nlist_sr[i]) );
}
}
}
for (i = 0; (i < eNL_NR); i++)
{
close_nblist(&(fr->nblists[n].nlist_sr[i]));
- close_nblist(&(fr->nblists[n].nlist_lr[i]));
}
}
}
-static gmx_inline void add_j_to_nblist(t_nblist *nlist, int j_atom, gmx_bool bLR)
+static gmx_inline void add_j_to_nblist(t_nblist *nlist, int j_atom)
{
int nrj = nlist->nrj;
if (gmx_debug_at)
{
- fprintf(debug, "Increasing %s nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
- bLR ? "LR" : "SR", nlist->ielec, nlist->ivdw, nlist->type, nlist->igeometry, nlist->maxnrj);
+ fprintf(debug, "Increasing SR nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
+ nlist->ielec, nlist->ivdw, nlist->type, nlist->igeometry, nlist->maxnrj);
}
srenew(nlist->jjnr, nlist->maxnrj);
static gmx_inline void add_j_to_nblist_cg(t_nblist *nlist,
int j_start, int j_end,
- t_excl *bexcl, gmx_bool i_is_j,
- gmx_bool bLR)
+ t_excl *bexcl, gmx_bool i_is_j)
{
int nrj = nlist->nrj;
int j;
nlist->maxnrj = over_alloc_small(nlist->nrj + 1);
if (gmx_debug_at)
{
- fprintf(debug, "Increasing %s nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
- bLR ? "LR" : "SR", nlist->ielec, nlist->ivdw, nlist->type, nlist->igeometry, nlist->maxnrj);
+ fprintf(debug, "Increasing SR nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
+ nlist->ielec, nlist->ivdw, nlist->type, nlist->igeometry, nlist->maxnrj);
}
srenew(nlist->jjnr, nlist->maxnrj);
t_excl bExcl[],
int shift,
t_forcerec * fr,
- gmx_bool bLR,
gmx_bool bDoVdW,
gmx_bool bDoCoul,
int solvent_opt);
t_excl bExcl[],
int shift,
t_forcerec * fr,
- gmx_bool bLR,
gmx_bool bDoVdW,
gmx_bool bDoCoul,
int solvent_opt)
{
nbl_ind = fr->gid2nblists[GID(igid, jgid, ngid)];
}
- if (bLR)
- {
- nlist = fr->nblists[nbl_ind].nlist_lr;
- }
- else
- {
- nlist = fr->nblists[nbl_ind].nlist_sr;
- }
+ nlist = fr->nblists[nbl_ind].nlist_sr;
if (iwater != esolNO)
{
if (!bDoCoul)
{
/* VdW only - only first atoms in each water interact */
- add_j_to_nblist(vdw, jj0, bLR);
+ add_j_to_nblist(vdw, jj0);
}
else
{
/* One entry for the entire water-water interaction */
if (!bDoVdW)
{
- add_j_to_nblist(coul_ww, jj0, bLR);
+ add_j_to_nblist(coul_ww, jj0);
}
else
{
- add_j_to_nblist(vdwc_ww, jj0, bLR);
+ add_j_to_nblist(vdwc_ww, jj0);
}
}
}
if (!bDoCoul)
{
/* VdW only - only first atoms in each water interact */
- add_j_to_nblist(vdw, jj0, bLR);
+ add_j_to_nblist(vdw, jj0);
}
else
{
/* One entry for the entire water-water interaction */
if (!bDoVdW)
{
- add_j_to_nblist(coul_ww, jj0, bLR);
+ add_j_to_nblist(coul_ww, jj0);
}
else
{
- add_j_to_nblist(vdwc_ww, jj0, bLR);
+ add_j_to_nblist(vdwc_ww, jj0);
}
}
}
{
if (charge[jj] != 0)
{
- add_j_to_nblist(coul, jj, bLR);
+ add_j_to_nblist(coul, jj);
}
}
}
{
if (bHaveVdW[type[jj]])
{
- add_j_to_nblist(vdw, jj, bLR);
+ add_j_to_nblist(vdw, jj);
}
}
}
{
if (charge[jj] != 0)
{
- add_j_to_nblist(vdwc, jj, bLR);
+ add_j_to_nblist(vdwc, jj);
}
else
{
- add_j_to_nblist(vdw, jj, bLR);
+ add_j_to_nblist(vdw, jj);
}
}
else if (charge[jj] != 0)
{
- add_j_to_nblist(coul, jj, bLR);
+ add_j_to_nblist(coul, jj);
}
}
}
{
if (charge[jj] != 0)
{
- add_j_to_nblist(coul, jj, bLR);
+ add_j_to_nblist(coul, jj);
}
}
else if (!bDoCoul_i)
{
if (bHaveVdW[type[jj]])
{
- add_j_to_nblist(vdw, jj, bLR);
+ add_j_to_nblist(vdw, jj);
}
}
else
{
if (charge[jj] != 0)
{
- add_j_to_nblist(vdwc, jj, bLR);
+ add_j_to_nblist(vdwc, jj);
}
else
{
- add_j_to_nblist(vdw, jj, bLR);
+ add_j_to_nblist(vdw, jj);
}
}
else if (charge[jj] != 0)
{
- add_j_to_nblist(coul, jj, bLR);
+ add_j_to_nblist(coul, jj);
}
}
}
{
if (charge[jj] != 0 || chargeB[jj] != 0)
{
- add_j_to_nblist(coul_free, jj, bLR);
+ add_j_to_nblist(coul_free, jj);
}
}
else if (!bDoCoul_i)
{
if (bHaveVdW[type[jj]] || bHaveVdW[typeB[jj]])
{
- add_j_to_nblist(vdw_free, jj, bLR);
+ add_j_to_nblist(vdw_free, jj);
}
}
else
{
if (charge[jj] != 0 || chargeB[jj] != 0)
{
- add_j_to_nblist(vdwc_free, jj, bLR);
+ add_j_to_nblist(vdwc_free, jj);
}
else
{
- add_j_to_nblist(vdw_free, jj, bLR);
+ add_j_to_nblist(vdw_free, jj);
}
}
else if (charge[jj] != 0 || chargeB[jj] != 0)
{
- add_j_to_nblist(coul_free, jj, bLR);
+ add_j_to_nblist(coul_free, jj);
}
}
}
/* This is done whether or not bWater is set */
if (charge[jj] != 0)
{
- add_j_to_nblist(coul, jj, bLR);
+ add_j_to_nblist(coul, jj);
}
}
else if (!bDoCoul_i_sol)
{
if (bHaveVdW[type[jj]])
{
- add_j_to_nblist(vdw, jj, bLR);
+ add_j_to_nblist(vdw, jj);
}
}
else
{
if (charge[jj] != 0)
{
- add_j_to_nblist(vdwc, jj, bLR);
+ add_j_to_nblist(vdwc, jj);
}
else
{
- add_j_to_nblist(vdw, jj, bLR);
+ add_j_to_nblist(vdw, jj);
}
}
else if (charge[jj] != 0)
{
- add_j_to_nblist(coul, jj, bLR);
+ add_j_to_nblist(coul, jj);
}
}
}
t_excl bExcl[],
int shift,
t_forcerec * fr,
- gmx_bool bLR,
gmx_bool gmx_unused bDoVdW,
gmx_bool gmx_unused bDoCoul,
int gmx_unused solvent_opt)
bNotEx = NOTEXCL(bExcl, i, jj);
if (bNotEx)
{
- add_j_to_nblist(coul, jj, bLR);
+ add_j_to_nblist(coul, jj);
}
}
}
t_excl bExcl[],
int shift,
t_forcerec * fr,
- gmx_bool bLR,
gmx_bool gmx_unused bDoVdW,
gmx_bool gmx_unused bDoCoul,
int gmx_unused solvent_opt)
{
nbl_ind = fr->gid2nblists[gid];
}
- if (bLR)
- {
- vdwc = &fr->nblists[nbl_ind].nlist_lr[eNL_VDWQQ];
- }
- else
- {
- vdwc = &fr->nblists[nbl_ind].nlist_sr[eNL_VDWQQ];
- }
+ vdwc = &fr->nblists[nbl_ind].nlist_sr[eNL_VDWQQ];
/* Make a new neighbor list for charge group icg.
* Currently simply one neighbor list is made with LJ and Coulomb.
/* Here we add the j charge group jcg to the list,
* exclusions are also added to the list.
*/
- add_j_to_nblist_cg(vdwc, index[jcg], index[jcg+1], bExcl, icg == jcg, bLR);
+ add_j_to_nblist_cg(vdwc, index[jcg], index[jcg+1], bExcl, icg == jcg);
}
}
if (nsbuf->nj + nrj > MAX_CG)
{
put_in_list(bHaveVdW, ngid, md, icg, jgid, nsbuf->ncg, nsbuf->jcg,
- cgs->index, bexcl, shift, fr, FALSE, TRUE, TRUE, fr->solvent_opt);
+ cgs->index, bexcl, shift, fr, TRUE, TRUE, fr->solvent_opt);
/* Reset buffer contents */
nsbuf->ncg = nsbuf->nj = 0;
}
if (nsbuf->ncg > 0)
{
put_in_list(bHaveVdW, ngid, md, icg, nn, nsbuf->ncg, nsbuf->jcg,
- cgs->index, bexcl, k, fr, FALSE, TRUE, TRUE, fr->solvent_opt);
+ cgs->index, bexcl, k, fr, TRUE, TRUE, fr->solvent_opt);
nsbuf->ncg = nsbuf->nj = 0;
}
}
****************************************************/
-static void get_cutoff2(t_forcerec *fr, gmx_bool bDoLongRange,
- real *rvdw2, real *rcoul2,
- real *rs2, real *rm2, real *rl2)
+static void get_cutoff2(t_forcerec *fr, real *rs2)
{
*rs2 = sqr(fr->rlist);
-
- if (bDoLongRange && fr->bTwinRange)
- {
- /* With plain cut-off or RF we need to make the list exactly
- * up to the cut-off and the cut-off's can be different,
- * so we can not simply set them to rlistlong.
- * To keep this code compatible with (exotic) old cases,
- * we also create lists up to rvdw/rcoulomb for PME and Ewald.
- * The interaction check should correspond to:
- * !ir_vdw/coulomb_might_be_zero_at_cutoff from inputrec.c.
- */
- if (((fr->vdwtype == evdwCUT || fr->vdwtype == evdwPME) &&
- fr->vdw_modifier == eintmodNONE) ||
- fr->rvdw <= fr->rlist)
- {
- *rvdw2 = sqr(fr->rvdw);
- }
- else
- {
- *rvdw2 = sqr(fr->rlistlong);
- }
- if (((fr->eeltype == eelCUT ||
- (EEL_RF(fr->eeltype) && fr->eeltype != eelRF_ZERO) ||
- fr->eeltype == eelPME ||
- fr->eeltype == eelEWALD) &&
- fr->coulomb_modifier == eintmodNONE) ||
- fr->rcoulomb <= fr->rlist)
- {
- *rcoul2 = sqr(fr->rcoulomb);
- }
- else
- {
- *rcoul2 = sqr(fr->rlistlong);
- }
- }
- else
- {
- /* Workaround for a gcc -O3 or -ffast-math problem */
- *rvdw2 = *rs2;
- *rcoul2 = *rs2;
- }
- *rm2 = std::min(*rvdw2, *rcoul2);
- *rl2 = std::max(*rvdw2, *rcoul2);
}
static void init_nsgrid_lists(t_forcerec *fr, int ngid, gmx_ns_t *ns)
{
- real rvdw2, rcoul2, rs2, rm2, rl2;
+ real rs2;
int j;
- get_cutoff2(fr, TRUE, &rvdw2, &rcoul2, &rs2, &rm2, &rl2);
+ get_cutoff2(fr, &rs2);
/* Short range buffers */
snew(ns->nl_sr, ngid);
/* Counters */
snew(ns->nsr, ngid);
- snew(ns->nlr_ljc, ngid);
- snew(ns->nlr_one, ngid);
-
- /* Always allocate both list types, since rcoulomb might now change with PME load balancing */
- /* Long range VdW and Coul buffers */
- snew(ns->nl_lr_ljc, ngid);
- /* Long range VdW or Coul only buffers */
- snew(ns->nl_lr_one, ngid);
for (j = 0; (j < ngid); j++)
{
snew(ns->nl_sr[j], MAX_CG);
- snew(ns->nl_lr_ljc[j], MAX_CG);
- snew(ns->nl_lr_one[j], MAX_CG);
}
if (debug)
{
fprintf(debug,
- "ns5_core: rs2 = %g, rm2 = %g, rl2 = %g (nm^2)\n",
- rs2, rm2, rl2);
+ "ns5_core: rs2 = %g (nm^2)\n",
+ rs2);
}
}
t_mdatoms *md,
put_in_list_t *put_in_list,
gmx_bool bHaveVdW[],
- gmx_bool bDoLongRange, gmx_bool bMakeQMMMnblist)
+ gmx_bool bMakeQMMMnblist)
{
gmx_ns_t *ns;
- int **nl_lr_ljc, **nl_lr_one, **nl_sr;
- int *nlr_ljc, *nlr_one, *nsr;
+ int **nl_sr;
+ int *nsr;
gmx_domdec_t *dd;
t_block *cgs = &(top->cgs);
int *cginfo = fr->cginfo;
int cg0, cg1, icg = -1, cgsnr, i0, igid, naaj, max_jcg;
int jcg0, jcg1, jjcg, cgj0, jgid;
int *grida, *gridnra, *gridind;
- gmx_bool rvdw_lt_rcoul, rcoul_lt_rvdw;
rvec *cgcm, grid_offset;
- real r2, rs2, rvdw2, rcoul2, rm2, rl2, XI, YI, ZI, tmp1, tmp2;
+ real r2, rs2, XI, YI, ZI, tmp1, tmp2;
int *i_egp_flags;
gmx_bool bDomDec, bTriclinicX, bTriclinicY;
ivec ncpddc;
cgsnr = cgs->nr;
- get_cutoff2(fr, bDoLongRange, &rvdw2, &rcoul2, &rs2, &rm2, &rl2);
-
- rvdw_lt_rcoul = (rvdw2 >= rcoul2);
- rcoul_lt_rvdw = (rcoul2 >= rvdw2);
-
- if (bMakeQMMMnblist)
- {
- rm2 = rl2;
- rs2 = rl2;
- }
+ get_cutoff2(fr, &rs2);
nl_sr = ns->nl_sr;
nsr = ns->nsr;
- nl_lr_ljc = ns->nl_lr_ljc;
- nl_lr_one = ns->nl_lr_one;
- nlr_ljc = ns->nlr_ljc;
- nlr_one = ns->nlr_one;
/* Unpack arrays */
cgcm = fr->cg_cm;
else
{
if (d == XX &&
- box[XX][XX] - fabs(box[YY][XX]) - fabs(box[ZZ][XX]) < std::sqrt(rl2))
+ box[XX][XX] - fabs(box[YY][XX]) - fabs(box[ZZ][XX]) < std::sqrt(rs2))
{
shp[d] = 2;
}
ZI = cgcm[icg][ZZ]+tz*box[ZZ][ZZ];
/* Calculate range of cells in Z direction that have the shift tz */
zgi = cell_z + tz*Nz;
- get_dx_dd(Nz, gridz, rl2, zgi, ZI-grid_offset[ZZ],
+ get_dx_dd(Nz, gridz, rs2, zgi, ZI-grid_offset[ZZ],
ncpddc[ZZ], sh0[ZZ], sh1[ZZ], &dz0, &dz1, dcz2);
if (dz0 > dz1)
{
{
ygi = cell_y + ty*Ny;
}
- get_dx_dd(Ny, gridy, rl2, ygi, YI-grid_offset[YY],
+ get_dx_dd(Ny, gridy, rs2, ygi, YI-grid_offset[YY],
ncpddc[YY], sh0[YY], sh1[YY], &dy0, &dy1, dcy2);
if (dy0 > dy1)
{
{
xgi = cell_x + tx*Nx;
}
- get_dx_dd(Nx, gridx, rl2, xgi, XI-grid_offset[XX],
+ get_dx_dd(Nx, gridx, rs2, xgi, XI-grid_offset[XX],
ncpddc[XX], sh0[XX], sh1[XX], &dx0, &dx1, dcx2);
if (dx0 > dx1)
{
for (nn = 0; (nn < ngid); nn++)
{
nsr[nn] = 0;
- nlr_ljc[nn] = 0;
- nlr_one[nn] = 0;
}
#ifdef NS5DB
fprintf(log, "shift: %2d, dx0,1: %2d,%2d, dy0,1: %2d,%2d, dz0,1: %2d,%2d\n",
#endif
for (dx = dx0; (dx <= dx1); dx++)
{
- tmp1 = rl2 - dcx2[dx];
+ tmp1 = rs2 - dcx2[dx];
for (dy = dy0; (dy <= dy1); dy++)
{
tmp2 = tmp1 - dcy2[dy];
(jjcg < max_jcg))
{
r2 = calc_dx2(XI, YI, ZI, cgcm[jjcg]);
- if (r2 < rl2)
+ if (r2 < rs2)
{
/* jgid = gid[cgsatoms[cgsindex[jjcg]]]; */
jgid = GET_CGINFO_GID(cginfo[jjcg]);
/* check energy group exclusions */
if (!(i_egp_flags[jgid] & EGP_EXCL))
{
- if (r2 < rs2)
- {
- if (nsr[jgid] >= MAX_CG)
- {
- /* Add to short-range list */
- put_in_list(bHaveVdW, ngid, md, icg, jgid,
- nsr[jgid], nl_sr[jgid],
- cgs->index, /* cgsatoms, */ bexcl,
- shift, fr, FALSE, TRUE, TRUE, fr->solvent_opt);
- nsr[jgid] = 0;
- }
- nl_sr[jgid][nsr[jgid]++] = jjcg;
- }
- else if (r2 < rm2)
- {
- if (nlr_ljc[jgid] >= MAX_CG)
- {
- /* Add to LJ+coulomb long-range list */
- put_in_list(bHaveVdW, ngid, md, icg, jgid,
- nlr_ljc[jgid], nl_lr_ljc[jgid], top->cgs.index,
- bexcl, shift, fr, TRUE, TRUE, TRUE, fr->solvent_opt);
- nlr_ljc[jgid] = 0;
- }
- nl_lr_ljc[jgid][nlr_ljc[jgid]++] = jjcg;
- }
- else
+ if (nsr[jgid] >= MAX_CG)
{
- if (nlr_one[jgid] >= MAX_CG)
- {
- /* Add to long-range list with only coul, or only LJ */
- put_in_list(bHaveVdW, ngid, md, icg, jgid,
- nlr_one[jgid], nl_lr_one[jgid], top->cgs.index,
- bexcl, shift, fr, TRUE, rvdw_lt_rcoul, rcoul_lt_rvdw, fr->solvent_opt);
- nlr_one[jgid] = 0;
- }
- nl_lr_one[jgid][nlr_one[jgid]++] = jjcg;
+ /* Add to short-range list */
+ put_in_list(bHaveVdW, ngid, md, icg, jgid,
+ nsr[jgid], nl_sr[jgid],
+ cgs->index, /* cgsatoms, */ bexcl,
+ shift, fr, TRUE, TRUE, fr->solvent_opt);
+ nsr[jgid] = 0;
}
+ nl_sr[jgid][nsr[jgid]++] = jjcg;
}
}
nns++;
{
put_in_list(bHaveVdW, ngid, md, icg, nn, nsr[nn], nl_sr[nn],
cgs->index, /* cgsatoms, */ bexcl,
- shift, fr, FALSE, TRUE, TRUE, fr->solvent_opt);
- }
-
- if (nlr_ljc[nn] > 0)
- {
- put_in_list(bHaveVdW, ngid, md, icg, nn, nlr_ljc[nn],
- nl_lr_ljc[nn], top->cgs.index,
- bexcl, shift, fr, TRUE, TRUE, TRUE, fr->solvent_opt);
- }
-
- if (nlr_one[nn] > 0)
- {
- put_in_list(bHaveVdW, ngid, md, icg, nn, nlr_one[nn],
- nl_lr_one[nn], top->cgs.index,
- bexcl, shift, fr, TRUE, rvdw_lt_rcoul, rcoul_lt_rvdw, fr->solvent_opt);
+ shift, fr, TRUE, TRUE, fr->solvent_opt);
}
}
}
gmx_groups_t *groups,
t_commrec *cr,
t_nrnb *nrnb, t_mdatoms *md,
- gmx_bool bFillGrid,
- gmx_bool bDoLongRangeNS)
+ gmx_bool bFillGrid)
{
t_block *cgs = &(top->cgs);
rvec box_size, grid_x0, grid_x1;
if (fr->ePBC != epbcNONE)
{
- if (sqr(fr->rlistlong) >= max_cutoff2(fr->ePBC, box))
+ if (sqr(fr->rlist) >= max_cutoff2(fr->ePBC, box))
{
gmx_fatal(FARGS, "One of the box vectors has become shorter than twice the cut-off length or box_yy-|box_zy| or box_zz has become smaller than the cut-off.");
}
if (!bGrid)
{
min_size = std::min(box_size[XX], std::min(box_size[YY], box_size[ZZ]));
- if (2*fr->rlistlong >= min_size)
+ if (2*fr->rlist >= min_size)
{
gmx_fatal(FARGS, "One of the box diagonal elements has become smaller than twice the cut-off length.");
}
}
/* Reset the neighbourlists */
- reset_neighbor_lists(fr, TRUE, TRUE);
+ reset_neighbor_lists(fr);
if (bGrid && bFillGrid)
{
cgs->nr, fr->cg_cm, grid_x0, grid_x1, &grid_dens);
grid_first(log, grid, NULL, NULL, box, grid_x0, grid_x1,
- fr->rlistlong, grid_dens);
+ fr->rlist, grid_dens);
}
start = 0;
nsearch = nsgrid_core(cr, fr, box, ngid, top,
grid, ns->bexcl, ns->bExcludeAlleg,
md, put_in_list, ns->bHaveVdW,
- bDoLongRangeNS, FALSE);
+ FALSE);
/* neighbour searching withouth QMMM! QM atoms have zero charge in
* the classical calculation. The charge-charge interaction
nsearch += nsgrid_core(cr, fr, box, ngid, top,
grid, ns->bexcl, ns->bExcludeAlleg,
md, put_in_list_qmmm, ns->bHaveVdW,
- bDoLongRangeNS, TRUE);
+ TRUE);
}
}
else
#endif
inc_nrnb(nrnb, eNR_NS, nsearch);
- /* inc_nrnb(nrnb,eNR_LR,fr->nlr); */
return nsearch;
}
gmx_groups_t *groups,
t_commrec *cr,
t_nrnb *nrnb, t_mdatoms *md,
- gmx_bool bFillGrid,
- gmx_bool bDoLongRangeNS);
+ gmx_bool bFillGrid);
/* Debugging routines from wnblist.c */
void grid_first(FILE *fplog, t_grid *grid,
gmx_domdec_t *dd, const gmx_ddbox_t *ddbox,
matrix box, rvec izones_x0, rvec izones_x1,
- real rlistlong, real grid_density)
+ real rlist, real grid_density)
{
int i, m;
- set_grid_sizes(box, izones_x0, izones_x1, rlistlong, dd, ddbox, grid, grid_density);
+ set_grid_sizes(box, izones_x0, izones_x1, rlist, dd, ddbox, grid, grid_density);
grid->ncells = grid->n[XX]*grid->n[YY]*grid->n[ZZ];
{
donb_flags |= GMX_NONBONDED_DO_POTENTIAL;
}
- if (flags & GMX_FORCE_DO_LR)
- {
- donb_flags |= GMX_NONBONDED_DO_LR;
- }
kernel_data.flags = donb_flags;
kernel_data.lambda = lambda;
int start, homenr;
double mu[2*DIM];
gmx_bool bStateChanged, bNS, bFillGrid, bCalcCGCM;
- gmx_bool bDoLongRange, bDoForces, bSepLRF, bUseGPU, bUseOrEmulGPU;
+ gmx_bool bDoForces, bUseGPU, bUseOrEmulGPU;
gmx_bool bDiffKernels = FALSE;
rvec vzero, box_diag;
float cycles_pme, cycles_force, cycles_wait_gpu;
bNS = (flags & GMX_FORCE_NS) && (fr->bAllvsAll == FALSE);
bFillGrid = (bNS && bStateChanged);
bCalcCGCM = (bFillGrid && !DOMAINDECOMP(cr));
- bDoLongRange = (fr->bTwinRange && bNS && (flags & GMX_FORCE_DO_LR));
bDoForces = (flags & GMX_FORCE_FORCES);
- bSepLRF = (bDoLongRange && bDoForces && (flags & GMX_FORCE_SEPLRF));
bUseGPU = fr->nbv->bUseGPU;
bUseOrEmulGPU = bUseGPU || (nbv->grp[0].kernel_type == nbnxnk8x8x8_PlainC);
}
/* Reset energies */
- reset_enerdata(fr, bNS, enerd, MASTER(cr));
+ reset_enerdata(enerd);
clear_rvecs(SHIFTS, fr->fshift);
if (DOMAINDECOMP(cr) && !(cr->duty & DUTY_PME))
/* Clear the short- and long-range forces */
clear_rvecs(fr->natoms_force_constr, f);
- if (bSepLRF && do_per_step(step, inputrec->nstcalclr))
- {
- clear_rvecs(fr->natoms_force_constr, fr->f_twin);
- }
clear_rvec(fr->vir_diag_posres);
}
/* Compute the bonded and non-bonded energies and optionally forces */
do_force_lowlevel(fr, inputrec, &(top->idef),
cr, nrnb, wcycle, mdatoms,
- x, hist, f, bSepLRF ? fr->f_twin : f, enerd, fcd, top, fr->born,
+ x, hist, f, enerd, fcd, top, fr->born,
bBornRadii, box,
inputrec->fepvals, lambda, graph, &(top->excls), fr->mu_tot,
flags, &cycles_pme);
- if (bSepLRF)
- {
- if (do_per_step(step, inputrec->nstcalclr))
- {
- /* Add the long range forces to the short range forces */
- for (i = 0; i < fr->natoms_force_constr; i++)
- {
- rvec_add(fr->f_twin[i], f[i], f[i]);
- }
- }
- }
-
cycles_force += wallcycle_stop(wcycle, ewcFORCE);
if (ed)
/* Communicate the forces */
wallcycle_start(wcycle, ewcMOVEF);
dd_move_f(cr->dd, f, fr->fshift);
- if (bSepLRF)
- {
- /* We should not update the shift forces here,
- * since f_twin is already included in f.
- */
- dd_move_f(cr->dd, fr->f_twin, NULL);
- }
wallcycle_stop(wcycle, ewcMOVEF);
}
spread_vsite_f(vsite, x, f, fr->fshift, FALSE, NULL, nrnb,
&top->idef, fr->ePBC, fr->bMolPBC, graph, box, cr);
wallcycle_stop(wcycle, ewcVSITESPREAD);
-
- if (bSepLRF)
- {
- wallcycle_start(wcycle, ewcVSITESPREAD);
- spread_vsite_f(vsite, x, fr->f_twin, NULL, FALSE, NULL,
- nrnb,
- &top->idef, fr->ePBC, fr->bMolPBC, graph, box, cr);
- wallcycle_stop(wcycle, ewcVSITESPREAD);
- }
}
if (flags & GMX_FORCE_VIRIAL)
int start, homenr;
double mu[2*DIM];
gmx_bool bStateChanged, bNS, bFillGrid, bCalcCGCM;
- gmx_bool bDoLongRangeNS, bDoForces, bSepLRF;
+ gmx_bool bDoForces;
float cycles_pme, cycles_force;
start = 0;
bStateChanged = (flags & GMX_FORCE_STATECHANGED);
bNS = (flags & GMX_FORCE_NS) && (fr->bAllvsAll == FALSE);
- /* Should we update the long-range neighborlists at this step? */
- bDoLongRangeNS = fr->bTwinRange && bNS;
/* Should we perform the long-range nonbonded evaluation inside the neighborsearching? */
bFillGrid = (bNS && bStateChanged);
bCalcCGCM = (bFillGrid && !DOMAINDECOMP(cr));
bDoForces = (flags & GMX_FORCE_FORCES);
- bSepLRF = ((inputrec->nstcalclr > 1) && bDoForces &&
- (flags & GMX_FORCE_SEPLRF) && (flags & GMX_FORCE_DO_LR));
if (bStateChanged)
{
}
/* Reset energies */
- reset_enerdata(fr, bNS, enerd, MASTER(cr));
+ reset_enerdata(enerd);
clear_rvecs(SHIFTS, fr->fshift);
if (bNS)
/* Do the actual neighbour searching */
ns(fplog, fr, box,
groups, top, mdatoms,
- cr, nrnb, bFillGrid,
- bDoLongRangeNS);
+ cr, nrnb, bFillGrid);
wallcycle_stop(wcycle, ewcNS);
}
/* Clear the short- and long-range forces */
clear_rvecs(fr->natoms_force_constr, f);
- if (bSepLRF && do_per_step(step, inputrec->nstcalclr))
- {
- clear_rvecs(fr->natoms_force_constr, fr->f_twin);
- }
clear_rvec(fr->vir_diag_posres);
}
/* Compute the bonded and non-bonded energies and optionally forces */
do_force_lowlevel(fr, inputrec, &(top->idef),
cr, nrnb, wcycle, mdatoms,
- x, hist, f, bSepLRF ? fr->f_twin : f, enerd, fcd, top, fr->born,
+ x, hist, f, enerd, fcd, top, fr->born,
bBornRadii, box,
inputrec->fepvals, lambda,
graph, &(top->excls), fr->mu_tot,
flags,
&cycles_pme);
- if (bSepLRF)
- {
- if (do_per_step(step, inputrec->nstcalclr))
- {
- /* Add the long range forces to the short range forces */
- for (i = 0; i < fr->natoms_force_constr; i++)
- {
- rvec_add(fr->f_twin[i], f[i], f[i]);
- }
- }
- }
-
cycles_force = wallcycle_stop(wcycle, ewcFORCE);
if (ed)
{
dd_move_f(cr->dd, fr->f_novirsum, NULL);
}
- if (bSepLRF)
- {
- /* We should not update the shift forces here,
- * since f_twin is already included in f.
- */
- dd_move_f(cr->dd, fr->f_twin, NULL);
- }
wallcycle_stop(wcycle, ewcMOVEF);
}
spread_vsite_f(vsite, x, f, fr->fshift, FALSE, NULL, nrnb,
&top->idef, fr->ePBC, fr->bMolPBC, graph, box, cr);
wallcycle_stop(wcycle, ewcVSITESPREAD);
-
- if (bSepLRF)
- {
- wallcycle_start(wcycle, ewcVSITESPREAD);
- spread_vsite_f(vsite, x, fr->f_twin, NULL, FALSE, NULL,
- nrnb,
- &top->idef, fr->ePBC, fr->bMolPBC, graph, box, cr);
- wallcycle_stop(wcycle, ewcVSITESPREAD);
- }
}
if (flags & GMX_FORCE_VIRIAL)
#include "gromacs/topology/mtop_util.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
//! Global max algorithm
{
fprintf(debug, "TPI cg %d, atoms %d-%d\n", cg_tp, a_tp0, a_tp1);
}
- if (a_tp1 - a_tp0 > 1 &&
- (inputrec->rlist < inputrec->rcoulomb ||
- inputrec->rlist < inputrec->rvdw))
- {
- gmx_fatal(FARGS, "Can not do TPI for multi-atom molecule with a twin-range cut-off");
- }
+
+ GMX_RELEASE_ASSERT(inputrec->rcoulomb <= inputrec->rlist && inputrec->rvdw <= inputrec->rlist, "Twin-range interactions are not supported with TPI");
+
snew(x_mol, a_tp1-a_tp0);
bDispCorr = (inputrec->eDispCorr != edispcNO);
* and the RF exclusion terms of the inserted molecule occur
* within a single charge group we can pass NULL for the graph.
* This also avoids shifts that would move charge groups
- * out of the box.
- *
- * Some checks above ensure than we can not have
- * twin-range interactions together with nstlist > 1,
- * therefore we do not need to remember the LR energies.
- */
+ * out of the box. */
/* Make do_force do a single node force calculation */
cr->nnodes = 1;
do_force(fplog, cr, inputrec,
state_global->lambda,
NULL, fr, NULL, mu_tot, t, NULL, NULL, FALSE,
GMX_FORCE_NONBONDED | GMX_FORCE_ENERGY |
- (bNS ? GMX_FORCE_DYNAMICBOX | GMX_FORCE_NS | GMX_FORCE_DO_LR : 0) |
+ (bNS ? GMX_FORCE_DYNAMICBOX | GMX_FORCE_NS : 0) |
(bStateChanged ? GMX_FORCE_STATECHANGED : 0));
cr->nnodes = nnodes;
bStateChanged = FALSE;
for (i = 0; i < ngid; i++)
{
sum_UgembU[e++] +=
- (enerd->grpp.ener[egBHAMSR][GID(i, gid_tp, ngid)] +
- enerd->grpp.ener[egBHAMLR][GID(i, gid_tp, ngid)])*embU;
+ enerd->grpp.ener[egBHAMSR][GID(i, gid_tp, ngid)]*embU;
}
}
else
for (i = 0; i < ngid; i++)
{
sum_UgembU[e++] +=
- (enerd->grpp.ener[egLJSR][GID(i, gid_tp, ngid)] +
- enerd->grpp.ener[egLJLR][GID(i, gid_tp, ngid)])*embU;
+ enerd->grpp.ener[egLJSR][GID(i, gid_tp, ngid)]*embU;
}
}
if (bDispCorr)
{
for (i = 0; i < ngid; i++)
{
- sum_UgembU[e++] +=
- (enerd->grpp.ener[egCOULSR][GID(i, gid_tp, ngid)] +
- enerd->grpp.ener[egCOULLR][GID(i, gid_tp, ngid)])*embU;
+ sum_UgembU[e++] += enerd->grpp.ener[egCOULSR][GID(i, gid_tp, ngid)] * embU;
}
if (bRFExcl)
{
return upd->xp;
}
-static void combine_forces(gmx_update_t upd,
- int nstcalclr,
- gmx_constr_t constr,
- t_inputrec *ir, t_mdatoms *md, t_idef *idef,
- t_commrec *cr,
- gmx_int64_t step,
- t_state *state, gmx_bool bMolPBC,
- int start, int nrend,
- rvec f[], rvec f_lr[],
- tensor *vir_lr_constr,
- t_nrnb *nrnb)
-{
- int i, d;
-
- /* f contains the short-range forces + the long range forces
- * which are stored separately in f_lr.
- */
-
- if (constr != NULL && vir_lr_constr != NULL &&
- !(ir->eConstrAlg == econtSHAKE && ir->epc == epcNO))
- {
- /* We need to constrain the LR forces separately,
- * because due to the different pre-factor for the SR and LR
- * forces in the update algorithm, we have to correct
- * the constraint virial for the nstcalclr-1 extra f_lr.
- * Constrain only the additional LR part of the force.
- */
- rvec *xp;
- real fac;
- int gf = 0;
-
- xp = get_xprime(state, upd);
-
- fac = (nstcalclr - 1)*ir->delta_t*ir->delta_t;
-
- for (i = 0; i < md->homenr; i++)
- {
- if (md->cFREEZE != NULL)
- {
- gf = md->cFREEZE[i];
- }
- for (d = 0; d < DIM; d++)
- {
- if ((md->ptype[i] != eptVSite) &&
- (md->ptype[i] != eptShell) &&
- !ir->opts.nFreeze[gf][d])
- {
- xp[i][d] = state->x[i][d] + fac*f_lr[i][d]*md->invmass[i];
- }
- else
- {
- xp[i][d] = state->x[i][d];
- }
- }
- }
- constrain(NULL, FALSE, FALSE, constr, idef, ir, cr, step, 0, 1.0, md,
- state->x, xp, xp, bMolPBC, state->box, state->lambda[efptBONDED], NULL,
- NULL, vir_lr_constr, nrnb, econqForce);
- }
-
- /* Add nstcalclr-1 times the LR force to the sum of both forces
- * and store the result in forces_lr.
- */
- for (i = start; i < nrend; i++)
- {
- for (d = 0; d < DIM; d++)
- {
- f_lr[i][d] = f[i][d] + (nstcalclr - 1)*f_lr[i][d];
- }
- }
-}
-
void update_constraints(FILE *fplog,
gmx_int64_t step,
real *dvdlambda, /* the contribution to be added to the bonded interactions */
t_inputrec *inputrec, /* input record and box stuff */
t_mdatoms *md,
t_state *state,
- gmx_bool bMolPBC,
rvec *f, /* forces on home particles */
- gmx_bool bDoLR,
- rvec *f_lr,
- tensor *vir_lr_constr,
t_fcdata *fcd,
gmx_ekindata_t *ekind,
matrix M,
gmx_bool bInitStep,
int UpdatePart,
t_commrec *cr, /* these shouldn't be here -- need to think about it */
- t_nrnb *nrnb,
- gmx_constr_t constr,
- t_idef *idef)
+ gmx_constr_t constr)
{
gmx_bool bNH, bPR, bDoConstr = FALSE;
double dt, alpha;
- rvec *force;
int start, homenr, nrend;
rvec *xprime;
int nth, th;
bNH = inputrec->etc == etcNOSEHOOVER;
bPR = ((inputrec->epc == epcPARRINELLORAHMAN) || (inputrec->epc == epcMTTK));
- if (bDoLR && inputrec->nstcalclr > 1)
- {
- GMX_RELEASE_ASSERT(!EI_VV(inputrec->eI), "The twin-range setup is not supported by velocity-Verlet integrators");
- /* Store the total force + nstcalclr-1 times the LR force
- * in forces_lr, so it can be used in a normal update algorithm
- * to produce twin time stepping.
- */
- combine_forces(upd,
- inputrec->nstcalclr, constr, inputrec, md, idef, cr,
- step, state, bMolPBC,
- start, nrend, f, f_lr, vir_lr_constr, nrnb);
- force = f_lr;
- }
- else
- {
- force = f;
- }
-
/* ############# START The update of velocities and positions ######### */
where();
dump_it_all(fplog, "Before update",
- state->natoms, state->x, xprime, state->v, force);
+ state->natoms, state->x, xprime, state->v, f);
if (inputrec->eI == eiSD2)
{
inputrec->opts.nFreeze,
md->invmass, md->ptype,
md->cFREEZE, md->cACC, md->cTC,
- state->x, xprime, state->v, force, M,
+ state->x, xprime, state->v, f, M,
bNH, bPR);
}
else
do_update_visc(start_th, end_th, dt,
ekind->tcstat, state->nosehoover_vxi,
md->invmass, md->ptype,
- md->cTC, state->x, xprime, state->v, force, M,
+ md->cTC, state->x, xprime, state->v, f, M,
state->box,
ekind->cosacc.cos_accel,
ekind->cosacc.vcos,
inputrec->opts.acc, inputrec->opts.nFreeze,
md->invmass, md->ptype,
md->cFREEZE, md->cACC, md->cTC,
- state->x, xprime, state->v, force,
+ state->x, xprime, state->v, f,
inputrec->opts.ngtc, inputrec->opts.ref_t,
bDoConstr, TRUE,
step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
inputrec->opts.acc, inputrec->opts.nFreeze,
md->invmass, md->ptype,
md->cFREEZE, md->cACC, md->cTC,
- state->x, xprime, state->v, force, state->sd_X,
+ state->x, xprime, state->v, f, state->sd_X,
inputrec->opts.tau_t,
TRUE, step, inputrec->ld_seed,
DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
do_update_bd(start_th, end_th, dt,
inputrec->opts.nFreeze, md->invmass, md->ptype,
md->cFREEZE, md->cTC,
- state->x, xprime, state->v, force,
+ state->x, xprime, state->v, f,
inputrec->bd_fric,
upd->sd->bd_rf,
step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
inputrec->opts.acc, inputrec->opts.nFreeze,
md->invmass, md->ptype,
md->cFREEZE, md->cACC,
- state->v, force,
+ state->v, f,
(bNH || bPR), state->veta, alpha);
break;
case etrtPOSITION:
t_inputrec *inputrec, /* input record and box stuff */
t_mdatoms *md,
t_state *state,
- gmx_bool bMolPBC,
rvec *f, /* forces on home particles */
- gmx_bool bDoLR,
- rvec *f_lr,
- tensor *vir_lr_constr,
t_fcdata *fcd,
gmx_ekindata_t *ekind,
matrix M,
gmx_bool bInitStep,
int bUpdatePart,
t_commrec *cr, /* these shouldn't be here -- need to think about it */
- t_nrnb *nrnb,
- gmx_constr *constr,
- t_idef *idef);
+ gmx_constr *constr);
/* Return TRUE if OK, FALSE in case of Shake Error */
};
enum {
- egCOULSR, egLJSR, egBHAMSR, egCOULLR, egLJLR, egBHAMLR,
+ egCOULSR, egLJSR, egBHAMSR,
egCOUL14, egLJ14, egGB, egNR
};
extern const char *egrp_nm[egNR+1];
/* Cut-Off stuff.
* Infinite cut-off's will be GMX_CUTOFF_INF (unlike in t_inputrec: 0).
*/
- real rlist, rlistlong;
+ real rlist;
/* Dielectric constant resp. multiplication factor for charges */
real zsquare, temp;
/* The allocation size of vectors of size natoms_force */
int nalloc_force;
- /* Twin Range stuff, f_twin has size natoms_force */
- gmx_bool bTwinRange;
- int nlr;
- rvec *f_twin;
- /* Constraint virial correction for multiple time
- stepping. Supported for the group scheme when not using
- velocity-Verlet integrators. */
- tensor vir_twin_constr;
-
/* Forces that should not enter into the virial summation:
* PPPM/PME/Ewald/posres
*/
(ir->ewald_geometry == eewg3DC || ir->epsilon_surface != 0));
}
-gmx_bool inputrecTwinRange(const t_inputrec *ir)
-{
- return (ir->rlist > 0 && (ir->rlistlong == 0 || ir->rlistlong > ir->rlist));
-}
-
gmx_bool inputrecElecField(const t_inputrec *ir)
{
return (ir->ex[XX].n > 0 || ir->ex[YY].n > 0 || ir->ex[ZZ].n > 0);
int andersen_seed; /* Random seed for Andersen thermostat (obsolete) */
real verletbuf_tol; /* Per atom pair energy drift tolerance (kJ/mol/ps/atom) for list buffer */
real rlist; /* short range pairlist cut-off (nm) */
- real rlistlong; /* long range pairlist cut-off (nm) */
- int nstcalclr; /* Frequency of evaluating direct space long-range interactions */
real rtpi; /* Radius for test particle insertion */
int coulombtype; /* Type of electrostatics treatment */
int coulomb_modifier; /* Modify the Coulomb interaction */
real scalefactor; /* factor for scaling the MM charges in QM calc.*/
/* Fields for removed features go here (better caching) */
- gmx_bool bAdress;
+ gmx_bool bAdress; // Whether AdResS is enabled - always false if a valid .tpr was read
+ gmx_bool useTwinRange; // Whether twin-range scheme is active - always false if a valid .tpr was read
} t_inputrec;
int ir_optimal_nstcalcenergy(const t_inputrec *ir);
/* Cut-off */
real rlist;
- real rlistlong;
/* PME/Ewald */
real ewaldcoeff_q;
cmp_rvec(fp, "inputrec->posres_comB", -1, ir1->posres_comB, ir2->posres_comB, ftol, abstol);
cmp_real(fp, "inputrec->verletbuf_tol", -1, ir1->verletbuf_tol, ir2->verletbuf_tol, ftol, abstol);
cmp_real(fp, "inputrec->rlist", -1, ir1->rlist, ir2->rlist, ftol, abstol);
- cmp_real(fp, "inputrec->rlistlong", -1, ir1->rlistlong, ir2->rlistlong, ftol, abstol);
- cmp_int(fp, "inputrec->nstcalclr", -1, ir1->nstcalclr, ir2->nstcalclr);
cmp_real(fp, "inputrec->rtpi", -1, ir1->rtpi, ir2->rtpi, ftol, abstol);
cmp_int(fp, "inputrec->coulombtype", -1, ir1->coulombtype, ir2->coulombtype);
cmp_int(fp, "inputrec->coulomb_modifier", -1, ir1->coulomb_modifier, ir2->coulomb_modifier);
F_LJC_PAIRS_NB,
F_LJ,
F_BHAM,
- F_LJ_LR,
- F_BHAM_LR,
+ F_LJ_LR_NOLONGERUSED,
+ F_BHAM_LR_NOLONGERUSED,
F_DISPCORR,
F_COUL_SR,
- F_COUL_LR,
+ F_COUL_LR_NOLONGERUSED,
F_RF_EXCL,
F_COUL_RECIP,
F_LJ_RECIP,
gmx_bool bConverged = TRUE, bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
gmx_bool bResetCountersHalfMaxH = FALSE;
gmx_bool bTemp, bPres, bTrotter;
- gmx_bool bUpdateDoLR;
real dvdl_constr;
rvec *cbuf = NULL;
int cbuf_nalloc = 0;
nstglobalcomm = 1;
}
- check_ir_old_tpx_versions(cr, fplog, ir, top_global);
-
nstglobalcomm = check_nstglobalcomm(fplog, cr, nstglobalcomm, ir);
bGStatEveryStep = (nstglobalcomm == 1);
}
}
- if (inputrecTwinRange(ir) && repl_ex_nst % ir->nstcalclr != 0)
- {
- /* We should exchange at nstcalclr steps to get correct integration */
- gmx_fatal(FARGS, "The replica exchange period (%d) is not divisible by nstcalclr (%d)", repl_ex_nst, ir->nstcalclr);
- }
-
if (ir->efep != efepNO)
{
/* Set free energy calculation frequency as the greatest common
- * denominator of nstdhdl and repl_ex_nst.
- * Check for nstcalclr with twin-range, since we need the long-range
- * contribution to the free-energy at the correct (nstcalclr) steps.
- */
+ * denominator of nstdhdl and repl_ex_nst. */
nstfep = ir->fepvals->nstdhdl;
if (ir->bExpanded)
{
- if (inputrecTwinRange(ir) &&
- ir->expandedvals->nstexpanded % ir->nstcalclr != 0)
- {
- gmx_fatal(FARGS, "nstexpanded should be divisible by nstcalclr");
- }
nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
}
if (repl_ex_nst > 0)
{
nstfep = gmx_greatest_common_divisor(repl_ex_nst, nstfep);
}
- /* We checked divisibility of repl_ex_nst and nstcalclr above */
- if (inputrecTwinRange(ir) && nstfep % ir->nstcalclr != 0)
- {
- gmx_incons("nstfep not divisible by nstcalclr");
- }
}
/* Be REALLY careful about what flags you set here. You CANNOT assume
{
gmx_fatal(FARGS, "Rerun trajectory frame step %d time %f does not contain a box, while pbc is used", rerun_fr.step, rerun_fr.time);
}
- if (max_cutoff2(ir->ePBC, rerun_fr.box) < sqr(fr->rlistlong))
+ if (max_cutoff2(ir->ePBC, rerun_fr.box) < sqr(fr->rlist))
{
gmx_fatal(FARGS, "Rerun trajectory frame step %d time %f has too small box dimensions", rerun_fr.step, rerun_fr.time);
}
}
else
{
- /* Determine whether or not to do Neighbour Searching and LR */
+ /* Determine whether or not to do Neighbour Searching */
bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
}
force_flags = (GMX_FORCE_STATECHANGED |
((inputrecDynamicBox(ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
GMX_FORCE_ALLFORCES |
- GMX_FORCE_SEPLRF |
(bCalcVir ? GMX_FORCE_VIRIAL : 0) |
(bCalcEner ? GMX_FORCE_ENERGY : 0) |
(bDoFEP ? GMX_FORCE_DHDL : 0)
);
- if (fr->bTwinRange)
- {
- if (do_per_step(step, ir->nstcalclr))
- {
- force_flags |= GMX_FORCE_DO_LR;
- }
- }
-
if (shellfc)
{
/* Now is the time to relax the shells */
trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
}
- update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC,
- f, FALSE, NULL, NULL, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, f, fcd,
ekind, M, upd, bInitStep, etrtVELOCITY1,
- cr, nrnb, constr, &top->idef);
+ cr, constr);
if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */
{
if (EI_VV(ir->eI))
{
/* velocity half-step update */
- update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
- FALSE, NULL, NULL, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, f, fcd,
ekind, M, upd, FALSE, etrtVELOCITY2,
- cr, nrnb, constr, &top->idef);
+ cr, constr);
}
/* Above, initialize just copies ekinh into ekin,
}
copy_rvecn(state->x, cbuf, 0, state->natoms);
}
- bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr) && !EI_VV(ir->eI));
- update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
- bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
- ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
+ update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ ekind, M, upd, bInitStep, etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state,
cr, nrnb, wcycle, upd, constr,
FALSE, bCalcVir);
- if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
- {
- /* Correct the virial for multiple time stepping */
- m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
- }
-
if (ir->eI == eiVVAK)
{
/* erase F_EKIN and F_TEMP here? */
/* now we know the scaling, we can compute the positions again again */
copy_rvecn(cbuf, state->x, 0, state->natoms);
- update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
- FALSE, NULL, NULL, fcd,
- ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
+ update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ ekind, M, upd, bInitStep, etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
/* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
fprintf(fp, "%s\n\n", buf);
}
ir->rlist = rlist_new;
- ir->rlistlong = rlist_new;
}
}
ls.cluster_size_i, ls.cluster_size_j);
}
ir->rlist = rlist_new;
- ir->rlistlong = rlist_new;
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff