% LocalWords: pymbar multinode subensemble Monte solute subst groupconcept GPU
% LocalWords: dodecahedron octahedron dodecahedra equilibration usinggroups nm
% LocalWords: topologies rlistlong CUDA GPUs rcoulomb SIMD BlueGene FPUs erfc
-% LocalWords: cutoffschemesupport unbuffered bondeds AdResS OpenMP ewald rtol
+% LocalWords: cutoffschemesupport unbuffered bondeds OpenMP ewald rtol
% LocalWords: verletdrift peptide RMS rescaling ergodicity ergodic discretized
% LocalWords: isothermal compressibility isotropically anisotropic iteratively
% LocalWords: incompressible integrations translational biomolecules NMA PCA
year = 2011
}
-@Article{DelleSite2007,
- title = {Some fundamental problems for an energy-conserving adaptive-resolution molecular dynamics scheme},
- volume = {76},
- journal = BTpre,
- author = {L. {Delle Site}},
- year = {2007}
-}
-
-@article{Praprotnik2005,
- title = {Adaptive resolution molecular-dynamics simulation: Changing the degrees of freedom on the fly},
- volume = {123},
- journal = BTJcp,
- author = {M. Praprotnik and L. {Delle Site} and K. Kremer},
- year = {2005},
- pages = {224106}
-}
-
-@article{Praprotnik2008,
- title = {Multiscale Simulation of Soft Matter: From Scale Bridging to Adaptive Resolution},
- volume = {59},
- journal = BTarpc,
- author = {Matej Praprotnik and Luigi {Delle Site} and Kurt Kremer},
- year = {2008},
- pages = {545--571}
-}
-
-@article{Praprotnik2011,
- title = {Statistical physics problems in adaptive resolution computer simulations of complex fluids},
- volume = {145},
- journal = {J. Stat. Phys.},
- author = {Matej Praprotnik and Simon Poblete and Kurt Kremer},
- year = {2011},
- pages = {946--966}
-}
-
-@article{Junghans2010,
- title = {A reference implementation of the adaptive resolution scheme in {ESPResSo}},
- volume = {181},
- journal = BTcpc,
- author = {Junghans, Christoph and Poblete, S.},
- year = {2010},
- pages = {1449--1454}
-}
-
-@article{Fritsch2012,
- author={S. Fritsch and C. Junghans and K. Kremer},
- title={Structure formation of toluene around C60: Implementation of the Adaptive Resolution Scheme (AdResS) into GROMACS},
- journal= BTjctc,
- year= {2012},
- volume = {8},
- pages = {398--403},
- OPTurl = {http://dx.doi.org/10.1063/1.463137}
-}
-
@article{Ruehle2009,
title = {Versatile {Object-Oriented} Toolkit for {Coarse-Graining} Applications},
volume = {5},
pages = {3211--3223}
}
-@article{Poblete2010,
- title = {Coupling different levels of resolution in molecular simulations},
- volume = {132},
- journal = BTjcp,
- author = {Poblete, S. and Praprotnik, Matej and Kremer, Kurt and {Delle Site}, Luigi},
- year = {2010},
- pages = {114101}
-}
-
-@article{Fritsch2012b,
- title = {Adaptive resolution molecular dynamics simulation through coupling to an internal particle reservoir},
- journal = BTprl,
- author = {Fritsch, S. and Poblete, S. and Junghans, C. and Ciccottii, G. and {Delle Site}, L. and Kremer, K.},
- year = {2012},
- volume = {108},
- pages = {170602}
-},
-
@article{wang_jpcb10,
Author = {Wang, Z.-J. and Deserno, M.},
Title = {A Systematically Coarse-Grained Solvent-Free Model for
What is needed as well is a transition state optimizer.
-\section{\normindex{Adaptive Resolution Scheme}}
-\newcommand{\adress}{AdResS}
-The adaptive resolution scheme~\cite{Praprotnik2005,Praprotnik2008} (\seeindex{\adress}{Adaptive Resolution Scheme}) couples two systems with different resolutions by a force interpolation scheme.
-In contrast to the mixed Quantum-Classical simulation techniques of the previous section, the number of high resolution particles is not fixed, but can vary over the simulation time.
-
-Below we discuss {\adress} for a double resolution (atomistic and coarse grained) representation of the same system. See \figref{adress} for illustration.
-The details of implementation described in this section were published in~\cite{Junghans2010,Fritsch2012}.
-
-\begin{figure}
-\begin{center}
-\includegraphics[width=0.5\textwidth]{plots/adress}
-\caption{A schematic illustration of the {\adress} method for water.}
-\label{fig:adress}
-\end{center}
-\end{figure}
-Every molecule needs a well-defined mapping point (usually the center of mass)
-but any other linear combination of particle coordinates is also sufficient.
-In the topology the mapping point is defined by a virtual site. The forces in the coarse-grained region are functions of the mapping point positions only.
-In this implementation molecules are modeled by charge groups or sets of charge groups, which actually allows one to have multiple mapping points per molecule. This can be useful for bigger molecules like polymers. In that case one has to also extend the AdResS description to bonded interactions~\cite{Praprotnik2011}, which will be implemented into \gromacs in one of the future versions.
-
-The force between two molecules is given by~\cite{Praprotnik2005}
-\footnote{Note that the equation obeys Newton's third law, which is not the case for other interpolation schemes~\cite{DelleSite2007}.}:
-\begin{equation}
-\vec{F}_{\alpha\beta}=w_\alpha w_\beta \vec{F}^\mathrm{ex,mol}_{\alpha\beta} + \left[1-w_\alpha w_\beta\right] \vec{F}^\mathrm{cg,mol}_{\alpha\beta}~,
-\label{eqn:interpolation}
-\end{equation}
-where $\alpha$ and $\beta$ label the two molecules and $w_\alpha$, $w_\beta$ are the adaptive weights of the two molecules.
-
-The first part, which represents the explicit interaction of the molecules, can be written as:
-\begin{equation}
-\vec{F}^\mathrm{ex,mol}_{\alpha\beta}=\sum_{i\in\alpha}\sum_{j\in\beta} \vec{F}^\mathrm{ex}_{ij}~,
-\end{equation}
-where $\vec{F}^\mathrm{ex}_{ij}$ is the force between the $i$th atom in $\alpha$th molecule and the $j$th atom in the $\beta$th molecule, which is given by an explicit force field.
-The second part of \eqnref{interpolation} comes from the coarse-grained interaction of the molecules.
-In \gromacs a slightly extended case is implemented:
-\begin{equation}
- \vec{F}_{\alpha\beta}=\sum_{i\in\alpha}\sum_{j\in\beta} w_i w_j \vec{F}^\mathrm{ex}_{ij} + \left[1-w_\alpha w_\beta\right] \vec{F}^\mathrm{cg,mol}_{\alpha\beta}~,
-\label{eqn:interpolation2}
-\end{equation}
-where $w_i$ and $w_j$ are atom-wise weights, which are determined by the {\tt adress-site} option. For {\tt adress-site} being the center of mass, atom $i$ has the weight of the center of mass of its \emph{charge group}.
-The weight $w_\alpha$ of molecule $\alpha$ is determined by the position of coarse-grained particle, which is constructed as a virtual site from the atomistic particles as specified in the topology.
-This extension allows one to perform all kind of AdResS variations, but the common case can be recovered by using a center of mass virtual site in the topology, {\tt adress-site=COM} and putting all atoms (except the virtual site representing the coarse-grained interaction) of a molecule into one charge group.
-For big molecules, it is sometimes useful to use an atom-based weight, which can be either be achieved by setting {\tt adress-site=atomperatom} or putting every atom into a separate charge group (the center of mass of a charge group with one atom is the atom itself).
-
-The coarse-grained force field $\vec{F}^\mathrm{cg}$ is usually derived from the atomistic system by structure-based coarse-graining (see \secref{cg-forcefields}). To specify which atoms belong to a coarse-grained representation, energy groups are used.
-Each coarse-grained interaction has to be associated with a specific energy group, which is why the virtual sites representing the coarse-grained interactions also have to be in different charge groups. The energy groups which are treated as coarse-grained interactions are then listed in {\tt adress_cg_grp_names}.
-The most important element of this interpolation (see \eqnref{interpolation} and \eqnref{interpolation2}) is the adaptive weighting function (for illustration see \figref{adress}):
-\begin{equation}
-w(x)=
-\left\{\begin{array}{c@{\;:\;}l}
-1&\mathrm{atomistic/explicit\;region}\\
-0<w<1&\mathrm{hybrid\;region}\\
-0&\mathrm{coarse-grained\;region}
-\end{array}\right.~,
-\label{equ:weighting}
-\end{equation}
-which has a value between 0 and 1.
-This definition of $w$ gives a purely explicit force in the explicit region and a purely coarse-grained force in the coarse-grained region,
-so essentially \eqnref{interpolation} only the hybrid region has mixed interactions which would not appear in a standard simulation.
-In {\gromacs}, a $\cos^2$-like function is implemented as a weighting function:
-\begin{equation}
-w(x)=
-\left\{\begin{array}{c@{\;:\;}r@{x}l}
-0&&>d_\mathrm{ex}+d_\mathrm{hy}\\
-\cos^2\left(\frac{\pi}{2d_\mathrm{hy}}(x-d_\mathrm{ex})\right)&d_\mathrm{ex}+d_\mathrm{hy}>&>d_\mathrm{ex}\\
-1&d_\mathrm{ex}>&
-\end{array}\right.~,
-\label{equ:wf}
-\end{equation}
-where $d_\mathrm{ex}$ and $d_\mathrm{hy}$ are the sizes of the explicit and the hybrid region, respectively.
-Depending on the physical interest of the research, other functions could be implemented as long as the following boundary conditions are fulfilled:
-The function is 1) continuous, 2) monotonic and 3) has zero derivatives at the boundaries.
-Spherical and one-dimensional splitting of the simulation box has been implemented ({\tt adress-type} option)
-and depending on this, the distance $x$ to the center of the explicit region is calculated as follows:
-\begin{equation}
-x=
-\left\{
-\begin{array}{c@{\;:\;}l}
- |(\vec{R}_\alpha-\vec{R}_\mathrm{ct})\cdot\hat{e}|&\mathrm{splitting\;in\;}\hat{e}\mathrm{\;direction}\\
-|\vec{R}_\alpha-\vec{R}_\mathrm{ct}|&\mathrm{spherical\;splitting}
-\end{array}
-\right.~,
-\end{equation}
-where $\vec{R}_\mathrm{ct}$ is the center of the explicit zone (defined by {\tt adress-reference-coords} option). $\vec{R}_\alpha$ is the mapping point of the $\alpha$th molecule. For the center of mass mapping, it is given by:
-\begin{equation}
-R_\alpha=\frac{\sum_{i\in\alpha}m_i r_i}{\sum_{i\in\alpha}m_i}
-\label{equ:com-def}
-\end{equation}
-Note that the value of the weighting function depends exclusively on the mapping of the molecule.
-
-The interpolation of forces (see \eqnref{interpolation2}) can produce inhomogeneities in the density and affect the structure of the system in the hybrid region.
-
-One way of reducing the density inhomogeneities is by the application of the so-called thermodynamic force (TF)~\cite{Poblete2010}.
-Such a force consists of a space-dependent external field applied in the hybrid region on the coarse-grained site of each molecule. It can be specified for each of the species of the system.
-The TF compensates the pressure profile~\cite{Fritsch2012b} that emerges under a homogeneous density profile. Therefore, it can correct the local density inhomogeneities in the hybrid region and it also allows the coupling of atomistic and coarse-grained representations which by construction have different pressures at the target density.
-The field can be determined by an iterative procedure, which is described in detail in the \href{http://code.google.com/p/votca/downloads/list?&q=manual}{manual} of the \normindex{VOTCA package}~\cite{ruehle2009}. Setting the {\tt adress-interface-correction} to {\tt thermoforce} enables the TF correction and\newline{\tt adress-tf-grp-names} defines the energy groups to act on.
-
-\subsection{Example: Adaptive resolution simulation of water}\label{subsec:adressexample}
-In this section the set up of an adaptive resolution simulation coupling atomistic SPC ~\cite{Berendsen81} water to its coarse-grained representation will be explained (as used in \cite{Fritsch2012b}).
-The following steps are required to setup the simulation:
-\begin{itemize}
-\item Perform a reference all-atom simulation
-\item Create a coarse-grained representation and save it as tabulated interaction function
-\item Create a hybrid topology for the SPC water
-\item Modify the atomistic coordinate file to include the coarse grained representation
-\item Define the geometry of the adaptive simulation in the grompp input file
-\item Create an index file
-\end{itemize}
-The coarse-grained representation of the interaction is stored as tabulated interaction function see \ssecref{userpot}. The convention is to use the $C^{(12)}$ columns with the $C^{(12)}$- coefficient set to 1. All other columns should be zero. The VOTCA manual has detailed instructions and a tutorial for SPC water on how to coarse-grain the interaction using various techniques.
-Here we named the coarse grained interaction CG, so the corresponding tabulated file is {\tt table_CG_CG.xvg}. To create the topology one can start from the atomistic topology file (e.g. share/gromacs/top/oplsaa.ff/spc.itp), we are assuming rigid water here. In the VOTCA tutorial the file is named {\tt hybrid_spc.itp}.
-The only difference to the atomistic topology is the addition of a coarse-grained virtual site:
-{\small
-\begin{verbatim}
-[ moleculetype ]
-; molname nrexcl
-SOL 2
-
-[ atoms ]
-; nr type resnr residue atom cgnr charge mass
- 1 opls_116 1 SOL OW 1 -0.82
- 2 opls_117 1 SOL HW1 1 0.41
- 3 opls_117 1 SOL HW2 1 0.41
- 4 CG 1 SOL CG 2 0
-
-
-[ settles ]
-; OW funct doh dhh
-1 1 0.1 0.16330
-
-[ exclusions ]
-1 2 3
-2 1 3
-3 1 2
-
-[ virtual_sites3 ]
-; Site from funct a d
-4 1 2 3 1 0.05595E+00 0.05595E+00
-\end{verbatim}}
-The virtual site type 3 with the specified coefficients places the virtual site in the center of mass of the molecule (for larger molecules virtual_sitesn has to be used).
-We now need to include our modified water model in the topology file and define the type {\tt CG}. In {\tt topol.top}:
-{\small
-\begin{verbatim}
-#include "ffoplsaa.itp"
-
-[ atomtypes ]
-;name mass charge ptype sigma epsilon
- CG 0.00000 0.0000 V 1 0.25
-
-#include "hybrid_spc.itp"
-
-[ system ]
-Adaptive water
-[ molecules ]
-SOL 8507
-\end{verbatim}}
-The $\sigma$ and $\epsilon$ values correspond to $C_6=1$ and $C_{12}=1$ and thus the table file should contain the coarse-grained interaction in either the $C_6$ or $C_{12}$ column. In the example the OPLS force field is used where $\sigma$ and $\epsilon$ are specified.
-Note that for force fields which define atomtypes directly in terms of $C_6$ and $C_{12}$, one can simply set $C_6=0$ and $C_{12}=1$. See section \ssecref{userpot} for more details on tabulated interactions. Since now the water molecule has a virtual site the coordinate file also needs to include that.
-{\small
-\begin{verbatim}
-adaptive water coordinates
-34028
- 1SOL OW 1 0.283 0.886 0.647
- 1SOL HW1 2 0.359 0.884 0.711
- 1SOL HW2 3 0.308 0.938 0.566
- 1SOL CG 4 0.289 0.889 0.646
- 1SOL OW 5 1.848 0.918 0.082
- 1SOL HW1 6 1.760 0.930 0.129
- 1SOL HW2 7 1.921 0.912 0.150
- 1SOL CG 8 1.847 0.918 0.088
- (...)
-\end{verbatim}}
-This file can be created manually or using the VOTCA tool {\tt csg_map } with the {\tt --hybrid} option.\\
-In the grompp input file the AdResS feature needs to be enabled and the geometry defined.
-{\small
-\begin{verbatim}
-(...)
-; AdResS relevant options
-energygrps = CG
-energygrp_table = CG CG
-
-; Method for doing Van der Waals
-vdw-type = user
-
-adress = yes
-adress_type = xsplit
-adress_ex_width = 1.5
-adress_hy_width = 1.5
-adress_interface_correction = off
-adress_reference_coords = 8 0 0
-adress_cg_grp_names = CG
-\end{verbatim}}
-
-Here we are defining an energy group {\tt CG} which consists of the coarse-grained virtual site.
-As discussed above, the coarse-grained interaction is usually tabulated. This requires the {\tt vdw-type} parameter to be set to {\tt user}. In the case where multi-component systems are coarse-grained, an energy group has to be defined for each component. Note that all the energy groups defining coarse-grained representations have to be listed again in {\tt adress_cg_grp_names} to distinguish them from regular energy groups.\\
-The index file has to be updated to have a group CG which includes all the coarse-grained virtual sites. This can be done easily using the {\tt make_ndx} tool of gromacs.
-
\section{Using VMD plug-ins for trajectory file I/O}
\index{VMD plug-ins}\index{trajectory file}{\gromacs} tools are able
to use the plug-ins found in an existing installation of
% LocalWords: nstrout rotangles rotslabs rottorque RMSD rmsdfit excl NH amine
% LocalWords: positionrestraint es SH phenylalanine solvated sh nanometer QM
% LocalWords: Lennard Buckingham UK Svensson ab vsitetop co UHF MP interatomic
-% LocalWords: AdResS adress cg grp coords TF VOTCA thermoforce tf SPC userpot
+% LocalWords: cg grp coords SPC userpot
% LocalWords: itp sitesn atomtypes ff csg ndx Tesla CHARMM AA gauss
% LocalWords: CMAP nocmap fators Monte performant lib LD DIR llllcc
% LocalWords: CMake homepage DEV overclocking GT dlopen vmd VMDDIR
free energy perturbed non-bondeds yes yes
energy group contributions yes only on CPU
energy group exclusions yes no
-AdResS multi-scale yes no
OpenMP multi-threading only PME all
native GPU support no yes
Coulomb PME yes yes
calculations are done.
-Adaptive Resolution Simulation
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-.. mdp:: adress
-
- (no)
- Decide whether the AdResS feature is turned on.
-
-.. mdp:: adress-type
-
- .. mdp-value:: Off
-
- Do an AdResS simulation with weight equal 1, which is equivalent
- to an explicit (normal) MD simulation. The difference to
- disabled AdResS is that the AdResS variables are still read-in
- and hence are defined.
-
- .. mdp-value:: Constant
-
- Do an AdResS simulation with a constant weight,
- :mdp:`adress-const-wf` defines the value of the weight
-
- .. mdp-value:: XSplit
-
- Do an AdResS simulation with simulation box split in
- x-direction, so basically the weight is only a function of the x
- coordinate and all distances are measured using the x coordinate
- only.
-
- .. mdp-value:: Sphere
-
- Do an AdResS simulation with spherical explicit zone.
-
-.. mdp:: adress-const-wf
-
- (1)
- Provides the weight for a constant weight simulation
- (:mdp:`adress-type` =Constant)
-
-.. mdp:: adress-ex-width
-
- (0)
- Width of the explicit zone, measured from
- :mdp:`adress-reference-coords`.
-
-.. mdp:: adress-hy-width
-
- (0)
- Width of the hybrid zone.
-
-.. mdp:: adress-reference-coords
-
- (0,0,0)
- Position of the center of the explicit zone. Periodic boundary
- conditions apply for measuring the distance from it.
-
-.. mdp:: adress-cg-grp-names
-
- The names of the coarse-grained energy groups. All other energy
- groups are considered explicit and their interactions will be
- automatically excluded with the coarse-grained groups.
-
-.. mdp:: adress-site
-
- The mapping point from which the weight is calculated.
-
- .. mdp-value:: COM
-
- The weight is calculated from the center of mass of each charge group.
-
- .. mdp-value:: COG
-
- The weight is calculated from the center of geometry of each charge group.
-
- .. mdp-value:: Atom
-
- The weight is calculated from the position of 1st atom of each charge group.
-
- .. mdp-value:: AtomPerAtom
-
- The weight is calculated from the position of each individual atom.
-
-.. mdp:: adress-interface-correction
-
- .. mdp-value:: Off
-
- Do not apply any interface correction.
-
- .. mdp-value:: thermoforce
-
- Apply thermodynamic force interface correction. The table can be
- specified using the ``-tabletf`` option of :ref:`gmx mdrun`. The
- table should contain the potential and force (acting on
- molecules) as function of the distance from
- :mdp:`adress-reference-coords`.
-
-.. mdp:: adress-tf-grp-names
-
- The names of the energy groups to which the thermoforce is applied
- if enabled in :mdp:`adress-interface-correction`. If no group is
- given the default table is applied.
-
-.. mdp:: adress-ex-forcecap
-
- (0)
- Cap the force in the hybrid region, useful for big molecules. 0
- disables force capping.
-
-
Computational Electrophysiology
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Use these options to switch on and control ion/water position exchanges in "Computational
reals and groups to your subroutine. Check the inputrec definition
in ``src/gromacs/legacyheaders/types/inputrec.h``
+Removed features
+^^^^^^^^^^^^^^^^
+
+This feature has been removed from |Gromacs|, but so that old
+:ref:`mdp` and :ref:`tpr` files cannot be mistakenly misused, we still
+parse this option. :ref:`gmx grompp` and :ref:`gmx mdrun` will issue a
+fatal error if this is set.
+
+.. mdp:: adress
+
+ (no)
+
.. _reference manual: gmx-manual-parent-dir_
"Aldert van Buuren",
"Rudi van Drunen",
"Anton Feenstra",
- "Sebastian Fritsch",
"Gerrit Groenhof",
"Christoph Junghans",
"Anca Hamuraru",
"Driving Forces for Adsorption of Amphiphilic Peptides to Air-Water Interface",
"J. Phys. Chem. B",
114, 2010, "11093" },
- { "Fritsch12",
- "S. Fritsch, C. Junghans and K. Kremer",
- "Adaptive molecular simulation study on structure formation of toluene around C60 using Gromacs",
- "J. Chem. Theo. Comp.",
- 8, 2012, "398" },
- { "Junghans10",
- "C. Junghans and S. Poblete",
- "A reference implementation of the adaptive resolution scheme in ESPResSo",
- "Comp. Phys. Comm.",
- 181, 2010, "1449" },
{ "Wang2010",
"H. Wang, F. Dommert, C.Holm",
"Optimizing working parameters of the smooth particle mesh Ewald algorithm in terms of accuracy and efficiency",
#include <cstring>
#include <algorithm>
+#include <vector>
#include "gromacs/fileio/filenm.h"
#include "gromacs/fileio/gmxfio.h"
tpxv_IntermolecularBondeds, /**< permit inter-molecular bonded interactions in the topology */
tpxv_CompElWithSwapLayerOffset, /**< added parameters for improved CompEl setups */
tpxv_CompElPolyatomicIonsAndMultipleIonTypes, /**< CompEl now can handle polyatomic ions and more than two types of ions */
+ tpxv_RemoveAdress, /**< removed support for AdResS */
tpxv_Count /**< the total number of tpxv versions */
};
gmx_fio_do_real(fio, ir->userreal3);
gmx_fio_do_real(fio, ir->userreal4);
- /* AdResS stuff */
- if (file_version >= 77)
+ /* AdResS is removed, but we need to be able to read old files,
+ and let mdrun refuse to run them */
+ if (file_version >= 77 && file_version < tpxv_RemoveAdress)
{
gmx_fio_do_gmx_bool(fio, ir->bAdress);
if (ir->bAdress)
{
- if (bRead)
- {
- snew(ir->adress, 1);
- }
- gmx_fio_do_int(fio, ir->adress->type);
- gmx_fio_do_real(fio, ir->adress->const_wf);
- gmx_fio_do_real(fio, ir->adress->ex_width);
- gmx_fio_do_real(fio, ir->adress->hy_width);
- gmx_fio_do_int(fio, ir->adress->icor);
- gmx_fio_do_int(fio, ir->adress->site);
- gmx_fio_do_rvec(fio, ir->adress->refs);
- gmx_fio_do_int(fio, ir->adress->n_tf_grps);
- gmx_fio_do_real(fio, ir->adress->ex_forcecap);
- gmx_fio_do_int(fio, ir->adress->n_energy_grps);
- gmx_fio_do_int(fio, ir->adress->do_hybridpairs);
+ int idum, numThermoForceGroups, numEnergyGroups;
+ real rdum;
+ rvec rvecdum;
+ gmx_fio_do_int(fio, idum);
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_int(fio, idum);
+ gmx_fio_do_int(fio, idum);
+ gmx_fio_do_rvec(fio, rvecdum);
+ gmx_fio_do_int(fio, numThermoForceGroups);
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_int(fio, numEnergyGroups);
+ gmx_fio_do_int(fio, idum);
- if (bRead)
- {
- snew(ir->adress->tf_table_index, ir->adress->n_tf_grps);
- }
- if (ir->adress->n_tf_grps > 0)
- {
- gmx_fio_ndo_int(fio, ir->adress->tf_table_index, ir->adress->n_tf_grps);
- }
- if (bRead)
+ if (numThermoForceGroups > 0)
{
- snew(ir->adress->group_explicit, ir->adress->n_energy_grps);
+ std::vector<int> idumn(numThermoForceGroups);
+ gmx_fio_ndo_int(fio, idumn.data(), idumn.size());
}
- if (ir->adress->n_energy_grps > 0)
+ if (numEnergyGroups > 0)
{
- gmx_fio_ndo_int(fio, ir->adress->group_explicit, ir->adress->n_energy_grps);
+ std::vector<int> idumn(numEnergyGroups);
+ gmx_fio_ndo_int(fio, idumn.data(), idumn.size());
}
}
}
pr_swap(fp, indent, ir->swap);
}
- /* AdResS PARAMETERS */
- PS("adress", EBOOL(ir->bAdress));
- if (ir->bAdress)
- {
- PS("adress-type", EADRESSTYPE(ir->adress->type));
- PR("adress-const-wf", ir->adress->const_wf);
- PR("adress-ex-width", ir->adress->ex_width);
- PR("adress-hy-width", ir->adress->hy_width);
- PR("adress-ex-forcecap", ir->adress->ex_forcecap);
- PS("adress-interface-correction", EADRESSICTYPE(ir->adress->icor));
- PS("adress-site", EADRESSSITETYPE(ir->adress->site));
- pr_rvec(fp, indent, "adress-reference-coords", ir->adress->refs, DIM, TRUE);
- PS("adress-do-hybridpairs", EBOOL(ir->adress->do_hybridpairs));
- }
-
/* USER-DEFINED THINGIES */
PI("userint1", ir->userint1);
PI("userint2", ir->userint2);
{ efXVG, "-dhdl", "dhdl", ffOPTWR },
{ efXVG, "-field", "field", ffOPTWR },
{ efXVG, "-table", "table", ffOPTRD },
- { efXVG, "-tabletf", "tabletf", ffOPTRD },
{ efXVG, "-tablep", "tablep", ffOPTRD },
{ efXVG, "-tableb", "table", ffOPTRD },
{ efTRX, "-rerun", "rerun", ffOPTRD },
"multiple_entries", "no", "use_last", NULL
};
-const char *eAdresstype_names[eAdressNR+1] = {
- "off", "constant", "xsplit", "sphere", NULL
-};
-
-const char *eAdressICtype_names[eAdressICNR+1] = {
- "off", "thermoforce", NULL
-};
-
-const char *eAdressSITEtype_names[eAdressSITENR+1] = {
- "com", "cog", "atom", "atomperatom", NULL
-};
-
const char *gmx_nblist_geometry_names[GMX_NBLIST_GEOMETRY_NR+1] = {
"Particle-Particle", "Water3-Particle", "Water3-Water3", "Water4-Particle", "Water4-Water4", "CG-CG", NULL
};
const char *gmx_nblist_interaction_names[GMX_NBLIST_INTERACTION_NR+1] = {
- "Standard", "Free_Energy", "Adress", NULL
+ "Standard", "Free_Energy", NULL
};
const char *gmx_nbkernel_elec_names[GMX_NBKERNEL_ELEC_NR+1] =
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2009 Christoph Junghans, Brad Lambeth.
- * Copyright (c) 2011 Christoph Junghans, Sebastian Fritsch.
- * Copyright (c) 2011,2012,2013,2014,2015, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include "nb_generic_adress.h"
-
-#include <cmath>
-
-#include "gromacs/gmxlib/nonbonded/nb_kernel.h"
-#include "gromacs/gmxlib/nonbonded/nonbonded.h"
-#include "gromacs/legacyheaders/nrnb.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/utility/fatalerror.h"
-
-#define ALMOST_ZERO 1e-30
-#define ALMOST_ONE 1-(1e-30)
-void
-gmx_nb_generic_adress_kernel(t_nblist * nlist,
- rvec * xx,
- rvec * ff,
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- nb_kernel_data_t * kernel_data,
- t_nrnb * nrnb)
-{
- int ntype, table_nelements, ielec, ivdw;
- real facel;
- int n, ii, is3, ii3, k, nj0, nj1, jnr, j3, ggid, nnn, n0;
- real shX, shY, shZ;
- real fscal, felec, fvdw, velec, vvdw, tx, ty, tz;
- real rinvsq;
- real iq;
- real qq, vctot;
- int nti, nvdwparam;
- int tj;
- real rt, r, eps, eps2, Y, F, Geps, Heps2, VV, FF, Fp, fijD, fijR;
- real rinvsix;
- real vvdwtot;
- real vvdw_rep, vvdw_disp;
- real ix, iy, iz, fix, fiy, fiz;
- real jx, jy, jz;
- real dx, dy, dz, rsq, rinv;
- real c6, c12, cexp1, cexp2, br;
- real * charge;
- real * shiftvec;
- real * vdwparam;
- int * type;
- real * fshift;
- real * velecgrp;
- real * vvdwgrp;
- real tabscale;
- real * VFtab;
- real * x;
- real * f;
- int ewitab;
- real ewtabscale, eweps, ewrt, ewtabhalfspace;
- real * ewtab;
- real rcoulomb2, rvdw, rvdw2, sh_dispersion, sh_repulsion;
- real rcutoff, rcutoff2;
- real d, d2, sw, dsw, rinvcorr;
- real elec_swV3, elec_swV4, elec_swV5, elec_swF2, elec_swF3, elec_swF4;
- real vdw_swV3, vdw_swV4, vdw_swV5, vdw_swF2, vdw_swF3, vdw_swF4;
- gmx_bool bExactElecCutoff, bExactVdwCutoff, bExactCutoff;
-
- real * wf;
- real weight_cg1;
- real weight_cg2;
- real weight_product;
- real hybscal; /* the multiplicator to the force for hybrid interactions*/
- real force_cap;
- gmx_bool bCG;
- int egp_nr;
-
- wf = mdatoms->wf;
-
- force_cap = fr->adress_ex_forcecap;
-
- x = xx[0];
- f = ff[0];
- ielec = nlist->ielec;
- ivdw = nlist->ivdw;
-
- fshift = fr->fshift[0];
- velecgrp = kernel_data->energygrp_elec;
- vvdwgrp = kernel_data->energygrp_vdw;
- tabscale = kernel_data->table_elec_vdw->scale;
- VFtab = kernel_data->table_elec_vdw->data;
-
- ewtab = fr->ic->tabq_coul_FDV0;
- ewtabscale = fr->ic->tabq_scale;
- ewtabhalfspace = 0.5/ewtabscale;
-
- rcoulomb2 = fr->rcoulomb*fr->rcoulomb;
- rvdw = fr->rvdw;
- rvdw2 = rvdw*rvdw;
- sh_dispersion = fr->ic->dispersion_shift.cpot;
- sh_repulsion = fr->ic->repulsion_shift.cpot;
-
- if (fr->coulomb_modifier == eintmodPOTSWITCH)
- {
- d = fr->rcoulomb-fr->rcoulomb_switch;
- elec_swV3 = -10.0/(d*d*d);
- elec_swV4 = 15.0/(d*d*d*d);
- elec_swV5 = -6.0/(d*d*d*d*d);
- elec_swF2 = -30.0/(d*d*d);
- elec_swF3 = 60.0/(d*d*d*d);
- elec_swF4 = -30.0/(d*d*d*d*d);
- }
- else
- {
- /* Avoid warnings from stupid compilers (looking at you, Clang!) */
- elec_swV3 = elec_swV4 = elec_swV5 = elec_swF2 = elec_swF3 = elec_swF4 = 0.0;
- }
- if (fr->vdw_modifier == eintmodPOTSWITCH)
- {
- d = fr->rvdw-fr->rvdw_switch;
- vdw_swV3 = -10.0/(d*d*d);
- vdw_swV4 = 15.0/(d*d*d*d);
- vdw_swV5 = -6.0/(d*d*d*d*d);
- vdw_swF2 = -30.0/(d*d*d);
- vdw_swF3 = 60.0/(d*d*d*d);
- vdw_swF4 = -30.0/(d*d*d*d*d);
- }
- else
- {
- /* Avoid warnings from stupid compilers (looking at you, Clang!) */
- vdw_swV3 = vdw_swV4 = vdw_swV5 = vdw_swF2 = vdw_swF3 = vdw_swF4 = 0.0;
- }
-
- bExactElecCutoff = (fr->coulomb_modifier != eintmodNONE) || fr->eeltype == eelRF_ZERO;
- bExactVdwCutoff = (fr->vdw_modifier != eintmodNONE);
- bExactCutoff = bExactElecCutoff || bExactVdwCutoff;
-
- if (bExactCutoff)
- {
- rcutoff = ( fr->rcoulomb > fr->rvdw ) ? fr->rcoulomb : fr->rvdw;
- rcutoff2 = rcutoff*rcutoff;
- }
- else
- {
- /* Fix warnings for stupid compilers */
- rcutoff2 = 1e30;
- }
-
- /* avoid compiler warnings for cases that cannot happen */
- nnn = 0;
- eps = 0.0;
- eps2 = 0.0;
-
- /* 3 VdW parameters for buckingham, otherwise 2 */
- nvdwparam = (ivdw == GMX_NBKERNEL_VDW_BUCKINGHAM) ? 3 : 2;
- table_nelements = 12;
-
- charge = mdatoms->chargeA;
- type = mdatoms->typeA;
- facel = fr->epsfac;
- shiftvec = fr->shift_vec[0];
- vdwparam = fr->nbfp;
- ntype = fr->ntype;
-
- for (n = 0; (n < nlist->nri); n++)
- {
- is3 = 3*nlist->shift[n];
- shX = shiftvec[is3];
- shY = shiftvec[is3+1];
- shZ = shiftvec[is3+2];
- nj0 = nlist->jindex[n];
- nj1 = nlist->jindex[n+1];
- ii = nlist->iinr[n];
- ii3 = 3*ii;
- ix = shX + x[ii3+0];
- iy = shY + x[ii3+1];
- iz = shZ + x[ii3+2];
- iq = facel*charge[ii];
- nti = nvdwparam*ntype*type[ii];
- vctot = 0;
- vvdwtot = 0;
- fix = 0;
- fiy = 0;
- fiz = 0;
-
- /* We need to find out if this i atom is part of an
- all-atom or CG energy group */
- egp_nr = mdatoms->cENER[ii];
- bCG = !fr->adress_group_explicit[egp_nr];
-
- weight_cg1 = wf[ii];
-
- if ((!bCG) && weight_cg1 < ALMOST_ZERO)
- {
- continue;
- }
-
- for (k = nj0; (k < nj1); k++)
- {
- jnr = nlist->jjnr[k];
- weight_cg2 = wf[jnr];
- weight_product = weight_cg1*weight_cg2;
-
- if (weight_product < ALMOST_ZERO)
- {
- /* if it's a explicit loop, skip this atom */
- if (!bCG)
- {
- continue;
- }
- else /* if it's a coarse grained loop, include this atom */
- {
- hybscal = 1.0;
- }
- }
- else if (weight_product >= ALMOST_ONE)
- {
-
- /* if it's a explicit loop, include this atom */
- if (!bCG)
- {
- hybscal = 1.0;
- }
- else /* if it's a coarse grained loop, skip this atom */
- {
- continue;
- }
- }
- /* both have double identity, get hybrid scaling factor */
- else
- {
- hybscal = weight_product;
-
- if (bCG)
- {
- hybscal = 1.0 - hybscal;
- }
- }
-
- j3 = 3*jnr;
- jx = x[j3+0];
- jy = x[j3+1];
- jz = x[j3+2];
- dx = ix - jx;
- dy = iy - jy;
- dz = iz - jz;
- rsq = dx*dx+dy*dy+dz*dz;
- rinv = gmx_invsqrt(rsq);
- rinvsq = rinv*rinv;
- felec = 0;
- fvdw = 0;
- velec = 0;
- vvdw = 0;
-
- if (bExactCutoff && rsq > rcutoff2)
- {
- continue;
- }
-
- if (ielec == GMX_NBKERNEL_ELEC_CUBICSPLINETABLE || ivdw == GMX_NBKERNEL_VDW_CUBICSPLINETABLE)
- {
- r = rsq*rinv;
- rt = r*tabscale;
- n0 = rt;
- eps = rt-n0;
- eps2 = eps*eps;
- nnn = table_nelements*n0;
- }
-
- /* Coulomb interaction. ielec==0 means no interaction */
- if (ielec != GMX_NBKERNEL_ELEC_NONE)
- {
- qq = iq*charge[jnr];
-
- switch (ielec)
- {
- case GMX_NBKERNEL_ELEC_NONE:
- break;
-
- case GMX_NBKERNEL_ELEC_COULOMB:
- /* Vanilla cutoff coulomb */
- velec = qq*rinv;
- felec = velec*rinvsq;
- break;
-
- case GMX_NBKERNEL_ELEC_REACTIONFIELD:
- /* Reaction-field */
- velec = qq*(rinv+fr->k_rf*rsq-fr->c_rf);
- felec = qq*(rinv*rinvsq-2.0*fr->k_rf);
- break;
-
- case GMX_NBKERNEL_ELEC_CUBICSPLINETABLE:
- /* Tabulated coulomb */
- Y = VFtab[nnn];
- F = VFtab[nnn+1];
- Geps = eps*VFtab[nnn+2];
- Heps2 = eps2*VFtab[nnn+3];
- Fp = F+Geps+Heps2;
- VV = Y+eps*Fp;
- FF = Fp+Geps+2.0*Heps2;
- velec = qq*VV;
- felec = -qq*FF*tabscale*rinv;
- break;
-
- case GMX_NBKERNEL_ELEC_GENERALIZEDBORN:
- /* GB */
- gmx_fatal(FARGS, "Death & horror! GB generic interaction not implemented.\n");
- break;
-
- case GMX_NBKERNEL_ELEC_EWALD:
- ewrt = rsq*rinv*ewtabscale;
- ewitab = ewrt;
- eweps = ewrt-ewitab;
- ewitab = 4*ewitab;
- felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
- rinvcorr = (fr->coulomb_modifier == eintmodPOTSHIFT) ? rinv-fr->ic->sh_ewald : rinv;
- velec = qq*(rinvcorr-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
- felec = qq*rinv*(rinvsq-felec);
- break;
-
- default:
- gmx_fatal(FARGS, "Death & horror! No generic coulomb interaction for ielec=%d.\n", ielec);
- break;
- }
- if (fr->coulomb_modifier == eintmodPOTSWITCH)
- {
- d = rsq*rinv-fr->rcoulomb_switch;
- d = (d > 0.0) ? d : 0.0;
- d2 = d*d;
- sw = 1.0+d2*d*(elec_swV3+d*(elec_swV4+d*elec_swV5));
- dsw = d2*(elec_swF2+d*(elec_swF3+d*elec_swF4));
- /* Apply switch function. Note that felec=f/r since it will be multiplied
- * by the i-j displacement vector. This means felec'=f'/r=-(v*sw)'/r=
- * -(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=felec*sw-v*dsw/r
- */
- felec = felec*sw - rinv*velec*dsw;
- /* Once we have used velec to update felec we can modify velec too */
- velec *= sw;
- }
- if (bExactElecCutoff)
- {
- felec = (rsq <= rcoulomb2) ? felec : 0.0;
- velec = (rsq <= rcoulomb2) ? velec : 0.0;
- }
- vctot += velec;
- } /* End of coulomb interactions */
-
-
- /* VdW interaction. ivdw==0 means no interaction */
- if (ivdw != GMX_NBKERNEL_VDW_NONE)
- {
- tj = nti+nvdwparam*type[jnr];
-
- switch (ivdw)
- {
- case GMX_NBKERNEL_VDW_NONE:
- break;
-
- case GMX_NBKERNEL_VDW_LENNARDJONES:
- /* Vanilla Lennard-Jones cutoff */
- c6 = vdwparam[tj];
- c12 = vdwparam[tj+1];
- rinvsix = rinvsq*rinvsq*rinvsq;
- vvdw_disp = c6*rinvsix;
- vvdw_rep = c12*rinvsix*rinvsix;
- fvdw = (vvdw_rep-vvdw_disp)*rinvsq;
- if (fr->vdw_modifier == eintmodPOTSHIFT)
- {
- vvdw = (vvdw_rep + c12*sh_repulsion)/12.0 - (vvdw_disp + c6*sh_dispersion)/6.0;
- }
- else
- {
- vvdw = vvdw_rep/12.0-vvdw_disp/6.0;
- }
- break;
-
- case GMX_NBKERNEL_VDW_BUCKINGHAM:
- /* Buckingham */
- c6 = vdwparam[tj];
- cexp1 = vdwparam[tj+1];
- cexp2 = vdwparam[tj+2];
-
- rinvsix = rinvsq*rinvsq*rinvsq;
- vvdw_disp = c6*rinvsix;
- br = cexp2*rsq*rinv;
- vvdw_rep = cexp1*std::exp(-br);
- fvdw = (br*vvdw_rep-vvdw_disp)*rinvsq;
- if (fr->vdw_modifier == eintmodPOTSHIFT)
- {
- vvdw = (vvdw_rep-cexp1*std::exp(-cexp2*rvdw)) - (vvdw_disp + c6*sh_dispersion)/6.0;
- }
- else
- {
- vvdw = vvdw_rep-vvdw_disp/6.0;
- }
- break;
-
- case GMX_NBKERNEL_VDW_CUBICSPLINETABLE:
- /* Tabulated VdW */
- c6 = vdwparam[tj];
- c12 = vdwparam[tj+1];
- Y = VFtab[nnn+4];
- F = VFtab[nnn+5];
- Geps = eps*VFtab[nnn+6];
- Heps2 = eps2*VFtab[nnn+7];
- Fp = F+Geps+Heps2;
- VV = Y+eps*Fp;
- FF = Fp+Geps+2.0*Heps2;
- vvdw_disp = c6*VV;
- fijD = c6*FF;
- Y = VFtab[nnn+8];
- F = VFtab[nnn+9];
- Geps = eps*VFtab[nnn+10];
- Heps2 = eps2*VFtab[nnn+11];
- Fp = F+Geps+Heps2;
- VV = Y+eps*Fp;
- FF = Fp+Geps+2.0*Heps2;
- vvdw_rep = c12*VV;
- fijR = c12*FF;
- fvdw = -(fijD+fijR)*tabscale*rinv;
- vvdw = vvdw_disp + vvdw_rep;
- break;
-
- default:
- gmx_fatal(FARGS, "Death & horror! No generic VdW interaction for ivdw=%d.\n", ivdw);
- break;
- }
- if (fr->vdw_modifier == eintmodPOTSWITCH)
- {
- d = rsq*rinv-fr->rvdw_switch;
- d = (d > 0.0) ? d : 0.0;
- d2 = d*d;
- sw = 1.0+d2*d*(vdw_swV3+d*(vdw_swV4+d*vdw_swV5));
- dsw = d2*(vdw_swF2+d*(vdw_swF3+d*vdw_swF4));
- /* See coulomb interaction for the force-switch formula */
- fvdw = fvdw*sw - rinv*vvdw*dsw;
- vvdw *= sw;
- }
- if (bExactVdwCutoff)
- {
- fvdw = (rsq <= rvdw2) ? fvdw : 0.0;
- vvdw = (rsq <= rvdw2) ? vvdw : 0.0;
- }
- vvdwtot += vvdw;
- } /* end VdW interactions */
-
- fscal = felec+fvdw;
-
- if (!bCG && force_cap > 0 && (std::abs(fscal) > force_cap))
- {
- fscal = force_cap*fscal/std::abs(fscal);
- }
-
- fscal *= hybscal;
-
- tx = fscal*dx;
- ty = fscal*dy;
- tz = fscal*dz;
- fix = fix + tx;
- fiy = fiy + ty;
- fiz = fiz + tz;
- f[j3+0] = f[j3+0] - tx;
- f[j3+1] = f[j3+1] - ty;
- f[j3+2] = f[j3+2] - tz;
- }
-
- f[ii3+0] = f[ii3+0] + fix;
- f[ii3+1] = f[ii3+1] + fiy;
- f[ii3+2] = f[ii3+2] + fiz;
- fshift[is3] = fshift[is3]+fix;
- fshift[is3+1] = fshift[is3+1]+fiy;
- fshift[is3+2] = fshift[is3+2]+fiz;
- ggid = nlist->gid[n];
- velecgrp[ggid] += vctot;
- vvdwgrp[ggid] += vvdwtot;
- }
- /* Estimate flops, average for generic adress kernel:
- * 14 flops per outer iteration
- * 54 flops per inner iteration
- */
- inc_nrnb(nrnb, eNR_NBKERNEL_GENERIC_ADRESS, nlist->nri*14 + nlist->jindex[n]*54);
-}
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2009 Christoph Junghans, Brad Lambeth.
- * Copyright (c) 2011 Christoph Junghans, Sebastian Fritsch.
- * Copyright (c) 2011,2012,2014,2015, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-
-#ifndef _nb_generic_adress_h_
-#define _nb_generic_adress_h_
-
-#include "gromacs/gmxlib/nonbonded/nb_kernel.h"
-#include "gromacs/legacyheaders/types/mdatom.h"
-#include "gromacs/legacyheaders/types/nblist.h"
-#include "gromacs/legacyheaders/types/nrnb.h"
-#include "gromacs/math/vectypes.h"
-
-void
-gmx_nb_generic_adress_kernel(t_nblist * nlist,
- rvec * xx,
- rvec * ff,
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- nb_kernel_data_t * kernel_data,
- t_nrnb * nrnb);
-
-#endif
+++ /dev/null
-#!/usr/bin/python
-#
-# This file is part of the GROMACS molecular simulation package.
-#
-# Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
-# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
-# and including many others, as listed in the AUTHORS file in the
-# top-level source directory and at http://www.gromacs.org.
-#
-# GROMACS is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public License
-# as published by the Free Software Foundation; either version 2.1
-# of the License, or (at your option) any later version.
-#
-# GROMACS is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with GROMACS; if not, see
-# http://www.gnu.org/licenses, or write to the Free Software Foundation,
-# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-#
-# If you want to redistribute modifications to GROMACS, please
-# consider that scientific software is very special. Version
-# control is crucial - bugs must be traceable. We will be happy to
-# consider code for inclusion in the official distribution, but
-# derived work must not be called official GROMACS. Details are found
-# in the README & COPYING files - if they are missing, get the
-# official version at http://www.gromacs.org.
-#
-# To help us fund GROMACS development, we humbly ask that you cite
-# the research papers on the package. Check out http://www.gromacs.org.
-
-import sys
-import os
-sys.path.append ( "../preprocessor" )
-from gmxpreprocess import gmxpreprocess
-
-# "The happiest programs are programs that write other programs."
-#
-#
-# This script controls the generation of Gromacs nonbonded kernels.
-#
-# We no longer generate kernels on-the-fly, so this file is not run
-# during a Gromacs compile - only when we need to update the kernels (=rarely).
-#
-# To maximize performance, each combination of interactions in Gromacs
-# has a separate nonbonded kernel without conditionals in the code.
-# To avoid writing hundreds of different routines for each architecture,
-# we instead use a custom preprocessor so we can encode the conditionals
-# and expand for-loops (e.g, for water-water interactions)
-# from a general kernel template. While that file will contain quite a
-# few preprocessor directives, it is still an order of magnitude easier
-# to maintain than ~200 different kernels (not to mention it avoids bugs).
-#
-# To actually generate the kernels, this program iteratively calls the
-# preprocessor with different define settings corresponding to all
-# combinations of coulomb/van-der-Waals/geometry options.
-#
-# A main goal in the design was to make this new generator _general_. For
-# this reason we have used a lot of different fields to identify a particular
-# kernel and interaction. Basically, each kernel will have a name like
-#
-# nbkernel_ElecXX_VdwYY_GeomZZ_VF_QQ()
-#
-# Where XX/YY/ZZ/VF are strings to identify what the kernel computes.
-#
-# Elec/Vdw describe the type of interaction for electrostatics and van der Waals.
-# The geometry settings correspond e.g. to water-water or water-particle kernels,
-# and finally the VF setting is V,F,or VF depending on whether we calculate
-# only the potential, only the force, or both of them. The final string (QQ)
-# is the architecture/language/optimization of the kernel.
-#
-Arch = 'c'
-
-# Explanation of the 'properties':
-#
-# It is cheap to compute r^2, and the kernels require various other functions of r for
-# different kinds of interaction. Depending on the needs of the kernel and the available
-# processor instructions, this will be done in different ways.
-#
-# 'rinv' means we need 1/r, which is calculated as 1/sqrt(r^2).
-# 'rinvsq' means we need 1/(r*r). This is calculated as rinv*rinv if we already did rinv, otherwise 1/r^2.
-# 'r' is similarly calculated as r^2*rinv when needed
-# 'table' means the interaction is tabulated, in which case we will calculate a table index before the interaction
-# 'shift' means the interaction will be modified by a constant to make it zero at the cutoff.
-# 'cutoff' means the interaction is set to 0.0 outside the cutoff
-#
-
-FileHeader = \
-'/*\n' \
-' * This file is part of the GROMACS molecular simulation package.\n' \
-' *\n' \
-' * Copyright (c) 2012, by the GROMACS development team, led by\n' \
-' * David van der Spoel, Berk Hess, Erik Lindahl, and including many\n' \
-' * others, as listed in the AUTHORS file in the top-level source\n' \
-' * directory and at http://www.gromacs.org.\n' \
-' *\n' \
-' * GROMACS is free software; you can redistribute it and/or\n' \
-' * modify it under the terms of the GNU Lesser General Public License\n' \
-' * as published by the Free Software Foundation; either version 2.1\n' \
-' * of the License, or (at your option) any later version.\n' \
-' *\n' \
-' * GROMACS is distributed in the hope that it will be useful,\n' \
-' * but WITHOUT ANY WARRANTY; without even the implied warranty of\n' \
-' * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU\n' \
-' * Lesser General Public License for more details.\n' \
-' *\n' \
-' * You should have received a copy of the GNU Lesser General Public\n' \
-' * License along with GROMACS; if not, see\n' \
-' * http://www.gnu.org/licenses, or write to the Free Software Foundation,\n' \
-' * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.\n' \
-' *\n' \
-' * If you want to redistribute modifications to GROMACS, please\n' \
-' * consider that scientific software is very special. Version\n' \
-' * control is crucial - bugs must be traceable. We will be happy to\n' \
-' * consider code for inclusion in the official distribution, but\n' \
-' * derived work must not be called official GROMACS. Details are found\n' \
-' * in the README & COPYING files - if they are missing, get the\n' \
-' * official version at http://www.gromacs.org.\n' \
-' *\n' \
-' * To help us fund GROMACS development, we humbly ask that you cite\n' \
-' * the research papers on the package. Check out http://www.gromacs.org.\n' \
-' */\n' \
-'/*\n' \
-' * Note: this file was generated by the GROMACS '+Arch+' kernel generator.\n' \
-' */\n'
-
-###############################################
-# ELECTROSTATICS
-# Interactions and flags for them
-###############################################
-ElectrostaticsList = {
- 'None' : [],
- 'Coulomb' : ['rinv','rinvsq'],
- 'ReactionField' : ['rinv','rinvsq'],
- 'GeneralizedBorn' : ['rinv','r'],
- 'CubicSplineTable' : ['rinv','r','table'],
- 'Ewald' : ['rinv','rinvsq','r'],
-}
-
-
-###############################################
-# VAN DER WAALS
-# Interactions and flags for them
-###############################################
-VdwList = {
- 'None' : [],
- 'LennardJones' : ['rinvsq'],
- 'Buckingham' : ['rinv','rinvsq','r'],
- 'CubicSplineTable' : ['rinv','r','table'],
-}
-
-
-###############################################
-# MODIFIERS
-# Different ways to adjust/modify interactions to conserve energy
-###############################################
-ModifierList = {
- 'None' : [],
- 'ExactCutoff' : ['exactcutoff'], # Zero the interaction outside the cutoff, used for reaction-field-zero
- 'PotentialShift' : ['shift','exactcutoff'],
- 'PotentialSwitch' : ['rinv','r','switch','exactcutoff']
-}
-
-
-###############################################
-# GEOMETRY COMBINATIONS
-###############################################
-GeometryNameList = [
- [ 'Particle' , 'Particle' ],
- [ 'Water3' , 'Particle' ],
- [ 'Water3' , 'Water3' ],
- [ 'Water4' , 'Particle' ],
- [ 'Water4' , 'Water4' ]
-]
-
-
-###############################################
-# POTENTIAL / FORCE
-###############################################
-VFList = [
- 'PotentialAndForce',
-# 'Potential', # Not used yet
- 'Force'
-]
-
-
-###############################################
-# GEOMETRY PROPERTIES
-###############################################
-# Dictionaries with lists telling which interactions are present
-# 1,2,3 means particles 1,2,3 (but not 0) have electrostatics!
-GeometryElectrostatics = {
- 'Particle' : [ 0 ],
- 'Particle2' : [ 0 , 1 ],
- 'Particle3' : [ 0 , 1 , 2 ],
- 'Particle4' : [ 0 , 1 , 2 , 3 ],
- 'Water3' : [ 0 , 1 , 2 ],
- 'Water4' : [ 1 , 2 , 3 ]
-}
-
-GeometryVdw = {
- 'Particle' : [ 0 ],
- 'Particle2' : [ 0 , 1 ],
- 'Particle3' : [ 0 , 1 , 2 ],
- 'Particle4' : [ 0 , 1 , 2 , 3 ],
- 'Water3' : [ 0 ],
- 'Water4' : [ 0 ]
-}
-
-
-###############################################
-# ADRESS PROPERTIES
-###############################################
-AdressList = [
- 'CG', 'EX'
-]
-
-# Dictionary to abbreviate all strings (mixed from all the lists)
-Abbreviation = {
- 'None' : 'None',
- 'Coulomb' : 'Coul',
- 'Ewald' : 'Ew',
- 'ReactionField' : 'RF',
- 'GeneralizedBorn' : 'GB',
- 'CubicSplineTable' : 'CSTab',
- 'LennardJones' : 'LJ',
- 'Buckingham' : 'Bham',
- 'PotentialShift' : 'Sh',
- 'PotentialSwitch' : 'Sw',
- 'ExactCutoff' : 'Cut',
- 'PotentialAndForce' : 'VF',
- 'Potential' : 'V',
- 'Force' : 'F',
- 'Water3' : 'W3',
- 'Water4' : 'W4',
- 'Particle' : 'P1',
- 'Particle2' : 'P2',
- 'Particle3' : 'P3',
- 'Particle4' : 'P4'
-}
-
-
-###############################################
-# Functions
-###############################################
-
-# Return a string with the kernel name from current settings
-def MakeKernelFileName(KernelElec,KernelElecMod,KernelVdw,KernelVdwMod,KernelGeom,KernelRes):
- ElecStr = 'Elec' + Abbreviation[KernelElec]
- if(KernelElecMod!='None'):
- ElecStr = ElecStr + Abbreviation[KernelElecMod]
- VdwStr = 'Vdw' + Abbreviation[KernelVdw]
- if(KernelVdwMod!='None'):
- VdwStr = VdwStr + Abbreviation[KernelVdwMod]
- GeomStr = 'Geom' + Abbreviation[KernelGeom[0]] + Abbreviation[KernelGeom[1]]
- return 'nb_kernel_' + ElecStr + '_' + VdwStr + '_' + GeomStr + '_' + KernelRes + '_' + Arch
-
-def MakeKernelName(KernelElec,KernelElecMod,KernelVdw,KernelVdwMod,KernelGeom,KernelVF,KernelRes):
- ElecStr = 'Elec' + Abbreviation[KernelElec]
- if(KernelElecMod!='None'):
- ElecStr = ElecStr + Abbreviation[KernelElecMod]
- VdwStr = 'Vdw' + Abbreviation[KernelVdw]
- if(KernelVdwMod!='None'):
- VdwStr = VdwStr + Abbreviation[KernelVdwMod]
- GeomStr = 'Geom' + Abbreviation[KernelGeom[0]] + Abbreviation[KernelGeom[1]]
- VFStr = Abbreviation[KernelVF]
- return 'nb_kernel_' + ElecStr + '_' + VdwStr + '_' + GeomStr + '_' + VFStr + '_' + KernelRes + '_' + Arch
-
-# Return a string with a declaration to use for the kernel;
-# this will be a sequence of string combinations as well as the actual function name
-# Dont worry about field widths - that is just pretty-printing for the header!
-def MakeKernelDecl(KernelName,KernelElec,KernelElecMod,KernelVdw,KernelVdwMod,KernelGeom,KernelOther,KernelVF):
- KernelStr = '\"'+KernelName+'\"'
- ArchStr = '\"'+Arch+'\"'
- ElecStr = '\"'+KernelElec+'\"'
- ElecModStr = '\"'+KernelElecMod+'\"'
- VdwStr = '\"'+KernelVdw+'\"'
- VdwModStr = '\"'+KernelVdwMod+'\"'
- GeomStr = '\"'+KernelGeom[0]+KernelGeom[1]+'\"'
- OtherStr = '\"'+KernelOther+'\"'
- VFStr = '\"'+KernelVF+'\"'
-
- ThisSpec = ArchStr+', '+ElecStr+', '+ElecModStr+', '+VdwStr+', '+VdwModStr+', '+GeomStr+', '+OtherStr+', '+VFStr
- ThisDecl = ' { '+KernelName+', '+KernelStr+', '+ThisSpec+' }'
- return ThisDecl
-
-
-# Returns 1 if this kernel should be created, 0 if we should skip it
-# This routine is not critical - it is not the end of the world if we create more kernels,
-# but since the number is pretty large we save both space and compile-time by reducing it a bit.
-def KeepKernel(KernelElec,KernelElecMod,KernelVdw,KernelVdwMod,KernelGeom,KernelVF):
-
- # No need for kernels without interactions
- if(KernelElec=='None' and KernelVdw=='None'):
- return 0
-
- # No need for modifiers without interactions
- if((KernelElec=='None' and KernelElecMod!='None') or (KernelVdw=='None' and KernelVdwMod!='None')):
- return 0
-
- # No need for LJ-only water optimization, or water optimization with implicit solvent.
- if('Water' in KernelGeom[0] and (KernelElec=='None' or 'GeneralizedBorn' in KernelElec)):
- return 0
-
- # Non-matching table settings are pointless
- if( ('Table' in KernelElec) and ('Table' in KernelVdw) and KernelElec!=KernelVdw ):
- return 0
-
- # Try to reduce the number of different switch/shift options to get a reasonable number of kernels
- # For electrostatics, reaction-field can use 'exactcutoff', and ewald can use switch or shift.
- if(KernelElecMod=='ExactCutoff' and KernelElec!='ReactionField'):
- return 0
- if(KernelElecMod in ['PotentialShift','PotentialSwitch'] and KernelElec!='Ewald'):
- return 0
- # For Vdw, we support switch and shift for Lennard-Jones/Buckingham
- if((KernelVdwMod=='ExactCutoff') or
- (KernelVdwMod in ['PotentialShift','PotentialSwitch'] and KernelVdw not in ['LennardJones','Buckingham'])):
- return 0
-
- # Choose either switch or shift and don't mix them...
- if((KernelElecMod=='PotentialShift' and KernelVdwMod=='PotentialSwitch') or
- (KernelElecMod=='PotentialSwitch' and KernelVdwMod=='PotentialShift')):
- return 0
-
- # Don't use a Vdw kernel with a modifier if the electrostatics one does not have one
- if(KernelElec!='None' and KernelElecMod=='None' and KernelVdwMod!='None'):
- return 0
-
- # Don't use an electrostatics kernel with a modifier if the vdw one does not have one,
- # unless the electrostatics one is reaction-field with exact cutoff.
- if(KernelVdw!='None' and KernelVdwMod=='None' and KernelElecMod!='None'):
- if(KernelElec=='ReactionField' and KernelVdw!='CubicSplineTable'):
- return 0
- elif(KernelElec!='ReactionField'):
- return 0
-
- return 1
-
-
-
-#
-# The preprocessor will automatically expand the interactions for water and other
-# geometries inside the kernel, but to get this right we need to setup a couple
-# of defines - we do them in a separate routine to keep the main loop clean.
-#
-# While this routine might look a bit complex it is actually quite straightforward,
-# and the best news is that you wont have to modify _anything_ for a new geometry
-# as long as you correctly define its Electrostatics/Vdw geometry in the lists above!
-#
-def SetDefines(KernelElec,KernelElecMod,KernelVdw,KernelVdwMod,KernelGeom,KernelVF,defines):
- # What is the _name_ for the i/j group geometry?
- igeometry = KernelGeom[0]
- jgeometry = KernelGeom[1]
- # define so we can access it in the source when the preprocessor runs
- defines['GEOMETRY_I'] = igeometry
- defines['GEOMETRY_J'] = jgeometry
-
- # For the i/j groups, extract a python list of which sites have electrostatics
- # For SPC/TIP3p this will be [1,1,1], while TIP4p (no elec on first site) will be [0,1,1,1]
- ielec = GeometryElectrostatics[igeometry]
- jelec = GeometryElectrostatics[jgeometry]
- # Zero out the corresponding lists in case we dont do Elec
- if(KernelElec=='None'):
- ielec = []
- jelec = []
-
- # Extract similar interaction lists for Vdw interactions (example for SPC: [1,0,0])
- iVdw = GeometryVdw[igeometry]
- jVdw = GeometryVdw[jgeometry]
-
- # Zero out the corresponding lists in case we dont do Vdw
- if(KernelVdw=='None'):
- iVdw = []
- jVdw = []
-
- # iany[] and jany[] contains lists of the particles actually used (for interactions) in this kernel
- iany = list(set(ielec+iVdw)) # convert to+from set to make elements unique
- jany = list(set(jelec+jVdw))
-
- defines['PARTICLES_ELEC_I'] = ielec
- defines['PARTICLES_ELEC_J'] = jelec
- defines['PARTICLES_VDW_I'] = iVdw
- defines['PARTICLES_VDW_J'] = jVdw
- defines['PARTICLES_I'] = iany
- defines['PARTICLES_J'] = jany
-
- # elecij,Vdwij are sets with pairs of particles for which the corresponding interaction is done
- # (and anyij again corresponds to either electrostatics or Vdw)
- elecij = []
- Vdwij = []
- anyij = []
-
- for i in ielec:
- for j in jelec:
- elecij.append([i,j])
-
- for i in iVdw:
- for j in jVdw:
- Vdwij.append([i,j])
-
- for i in iany:
- for j in jany:
- if [i,j] in elecij or [i,j] in Vdwij:
- anyij.append([i,j])
-
- defines['PAIRS_IJ'] = anyij
-
- # Make an 2d list-of-distance-properties-to-calculate for i,j
- ni = max(iany)+1
- nj = max(jany)+1
- # Each element properties[i][j] is an empty list
- properties = [ [ [] for j in range(0,nj) ] for i in range (0,ni) ]
- # Add properties to each set
- for i in range(0,ni):
- for j in range(0,nj):
- if [i,j] in elecij:
- properties[i][j] = properties[i][j] + ['electrostatics'] + ElectrostaticsList[KernelElec] + ModifierList[KernelElecMod]
- if [i,j] in Vdwij:
- properties[i][j] = properties[i][j] + ['vdw'] + VdwList[KernelVdw] + ModifierList[KernelVdwMod]
- # Add rinv if we need r
- if 'r' in properties[i][j]:
- properties[i][j] = properties[i][j] + ['rinv']
- # Add rsq if we need rinv or rinsq
- if 'rinv' in properties[i][j] or 'rinvsq' in properties[i][j]:
- properties[i][j] = properties[i][j] + ['rsq']
-
- defines['INTERACTION_FLAGS'] = properties
-
-
-
-def PrintStatistics(ratio):
- ratio = 100.0*ratio
- print '\rGenerating %s nonbonded kernels... %5.1f%%' % (Arch,ratio),
- sys.stdout.flush()
-
-
-
-defines = {}
-kerneldecl = []
-
-cnt = 0.0
-nelec = len(ElectrostaticsList)
-nVdw = len(VdwList)
-nmod = len(ModifierList)
-ngeom = len(GeometryNameList)
-
-ntot = nelec*nmod*nVdw*nmod*ngeom
-
-numKernels = 0
-
-fpdecl = open('nb_kernel_' + Arch + '.c','w')
-fpdecl.write( FileHeader )
-fpdecl.write( '#include "gmxpre.h"\n\n' )
-fpdecl.write( '#include "../nb_kernel.h"\n\n' )
-
-for KernelElec in ElectrostaticsList:
- defines['KERNEL_ELEC'] = KernelElec
-
- for KernelElecMod in ModifierList:
- defines['KERNEL_MOD_ELEC'] = KernelElecMod
-
- for KernelVdw in VdwList:
- defines['KERNEL_VDW'] = KernelVdw
-
- for KernelVdwMod in ModifierList:
- defines['KERNEL_MOD_VDW'] = KernelVdwMod
-
- for KernelGeom in GeometryNameList:
-
- cnt += 1
- KernelFilename = MakeKernelFileName(KernelElec,KernelElecMod,KernelVdw,KernelVdwMod,KernelGeom) + '.c'
- fpkernel = open(KernelFilename,'w')
- defines['INCLUDE_HEADER'] = 1 # Include header first time in new file
- DoHeader = 1
-
- for KernelVF in VFList:
-
- KernelName = MakeKernelName(KernelElec,KernelElecMod,KernelVdw,KernelVdwMod,KernelGeom,KernelVF)
-
- defines['KERNEL_NAME'] = KernelName
- defines['KERNEL_VF'] = KernelVF
-
- # Check if this is a valid/sane/usable combination
- if not KeepKernel(KernelElec,KernelElecMod,KernelVdw,KernelVdwMod,KernelGeom,KernelVF):
- continue;
-
- # The overall kernel settings determine what the _kernel_ calculates, but for the water
- # kernels this does not mean that every pairwise interaction has e.g. Vdw interactions.
- # This routine sets defines of what to calculate for each pair of particles in those cases.
- SetDefines(KernelElec,KernelElecMod,KernelVdw,KernelVdwMod,KernelGeom,KernelVF,defines)
-
- if(DoHeader==1):
- fpkernel.write( FileHeader )
-
- gmxpreprocess('nb_kernel_template_' + Arch + '.pre', KernelName+'.tmp' , defines, force=1,contentType='C')
- numKernels = numKernels + 1
-
- defines['INCLUDE_HEADER'] = 0 # Header has been included once now
- DoHeader=0
-
- # Append temp file contents to the common kernelfile
- fptmp = open(KernelName+'.tmp','r')
- fpkernel.writelines(fptmp.readlines())
- fptmp.close()
- os.remove(KernelName+'.tmp')
-
- # Add a declaration for this kernel
- fpdecl.write('nb_kernel_t ' + KernelName + ';\n');
-
- # Add declaration to the buffer
- KernelOther=''
- kerneldecl.append(MakeKernelDecl(KernelName,KernelElec,KernelElecMod,KernelVdw,KernelVdwMod,KernelGeom,KernelOther,KernelVF))
-
- filesize = fpkernel.tell()
- fpkernel.close()
- if(filesize==0):
- os.remove(KernelFilename)
-
- PrintStatistics(cnt/ntot)
- pass
- pass
- pass
- pass
-pass
-
-# Write out the list of settings and corresponding kernels to the declaration file
-fpdecl.write( '\n\n' )
-fpdecl.write( 'nb_kernel_info_t\n' )
-fpdecl.write( ' kernellist_'+Arch+'[] =\n' )
-fpdecl.write( '{\n' )
-for decl in kerneldecl[0:-1]:
- fpdecl.write( decl + ',\n' )
-fpdecl.write( kerneldecl[-1] + '\n' )
-fpdecl.write( '};\n\n' )
-fpdecl.write( 'int\n' )
-fpdecl.write( ' kernellist_'+Arch+'_size = sizeof(kernellist_'+Arch+')/sizeof(kernellist_'+Arch+'[0]);\n')
-fpdecl.close()
+++ /dev/null
-/* #if 0 */
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2012,2014,2015, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#error This file must be processed with the Gromacs pre-preprocessor
-/* #endif */
-/* #if INCLUDE_HEADER */
-#include "gmxpre.h"
-
-#include "config.h"
-
-#include <math.h>
-
-#include "../nb_kernel.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/legacyheaders/nrnb.h"
-/* #endif */
-
-#define ALMOST_ZERO 1e-30
-#define ALMOST_ONE 1-(1e-30)
-
-/* ## List of variables set by the generating script: */
-/* ## */
-/* ## Setttings that apply to the entire kernel: */
-/* ## KERNEL_ELEC: String, choice for electrostatic interactions */
-/* ## KERNEL_VDW: String, choice for van der Waals interactions */
-/* ## KERNEL_NAME: String, name of this kernel */
-/* ## KERNEL_RESOLUTION: String, choice for resoultion */
-/* ## KERNEL_VF: String telling if we calculate potential, force, or both */
-/* ## GEOMETRY_I/GEOMETRY_J: String, name of each geometry, e.g. 'Water3' or '1Particle' */
-/* ## */
-/* ## Setttings that apply to particles in the outer (I) or inner (J) loops: */
-/* ## PARTICLES_I[]/ Arrays with lists of i/j particles to use in kernel. It is */
-/* ## PARTICLES_J[]: just [0] for particle geometry, but can be longer for water */
-/* ## PARTICLES_ELEC_I[]/ Arrays with lists of i/j particle that have electrostatics */
-/* ## PARTICLES_ELEC_J[]: interactions that should be calculated in this kernel. */
-/* ## PARTICLES_VDW_I[]/ Arrays with the list of i/j particle that have VdW */
-/* ## PARTICLES_VDW_J[]: interactions that should be calculated in this kernel. */
-/* ## */
-/* ## Setttings for pairs of interactions (e.g. 2nd i particle against 1st j particle) */
-/* ## PAIRS_IJ[]: Array with (i,j) tuples of pairs for which interactions */
-/* ## should be calculated in this kernel. Zero-charge particles */
-/* ## do not have interactions with particles without vdw, and */
-/* ## Vdw-only interactions are not evaluated in a no-vdw-kernel. */
-/* ## INTERACTION_FLAGS[][]: 2D matrix, dimension e.g. 3*3 for water-water interactions. */
-/* ## For each i-j pair, the element [I][J] is a list of strings */
-/* ## defining properties/flags of this interaction. Examples */
-/* ## include 'electrostatics'/'vdw' if that type of interaction */
-/* ## should be evaluated, 'rsq'/'rinv'/'rinvsq' if those values */
-/* ## are needed, and 'exactcutoff' or 'shift','switch' to */
-/* ## decide if the force/potential should be modified. This way */
-/* ## we only calculate values absolutely needed for each case. */
-
-/* ## Calculate the size and offset for (merged/interleaved) table data */
-
-/* #if ('CubicSplineTable' in [KERNEL_ELEC,KERNEL_VDW]) or KERNEL_VF=='PotentialAndForce' */
-/* #define TABLE_POINT_SIZE 4 */
-/* #else */
-/* #define TABLE_POINT_SIZE 2 */
-/* #endif */
-
-/* #if 'Table' in KERNEL_ELEC and 'Table' in KERNEL_VDW */
-/* #define TABLE_NINTERACTIONS 3 */
-/* #define TABLE_VDW_OFFSET TABLE_POINT_SIZE */
-/* #elif 'Table' in KERNEL_ELEC */
-/* #define TABLE_NINTERACTIONS 1 */
-/* #elif 'Table' in KERNEL_VDW */
-/* #define TABLE_NINTERACTIONS 2 */
-/* #define TABLE_VDW_OFFSET 0 */
-/* #else */
-/* #define TABLE_NINTERACTIONS 0 */
-/* #endif */
-
-/* #if 'Buckingham' in KERNEL_VDW */
-/* #define NVDWPARAM 3 */
-/* #else */
-/* #define NVDWPARAM 2 */
-/* #endif */
-
-/*
- * Gromacs nonbonded kernel: {KERNEL_NAME}
- * Electrostatics interaction: {KERNEL_ELEC}
- * VdW interaction: {KERNEL_VDW}
- * Geometry: {GEOMETRY_I}-{GEOMETRY_J}
- * Calculate force/pot: {KERNEL_VF}
- */
-void
-{KERNEL_NAME}
- (t_nblist * gmx_restrict nlist,
- rvec * gmx_restrict xx,
- rvec * gmx_restrict ff,
- t_forcerec * gmx_restrict fr,
- t_mdatoms * gmx_restrict mdatoms,
- nb_kernel_data_t * gmx_restrict kernel_data,
- t_nrnb * gmx_restrict nrnb)
-{
- /* ## Not all variables are used for all kernels, but any optimizing compiler fixes that, */
- /* ## so there is no point in going to extremes to exclude variables that are not needed. */
- int i_shift_offset,i_coord_offset,j_coord_offset;
- int j_index_start,j_index_end;
- int nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
- real shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
- int *iinr,*jindex,*jjnr,*shiftidx,*gid;
- real *shiftvec,*fshift,*x,*f;
- /* #for I in PARTICLES_I */
- int vdwioffset{I};
- real ix{I},iy{I},iz{I},fix{I},fiy{I},fiz{I},iq{I},isai{I};
- /* #endfor */
- /* #for J in PARTICLES_J */
- int vdwjidx{J};
- real jx{J},jy{J},jz{J},fjx{J},fjy{J},fjz{J},jq{J},isaj{J};
- /* #endfor */
- /* #for I,J in PAIRS_IJ */
- real dx{I}{J},dy{I}{J},dz{I}{J},rsq{I}{J},rinv{I}{J},rinvsq{I}{J},r{I}{J},qq{I}{J},c6_{I}{J},c12_{I}{J},cexp1_{I}{J},cexp2_{I}{J};
- /* #endfor */
- /* #if KERNEL_ELEC != 'None' */
- real velec,felec,velecsum,facel,crf,krf,krf2;
- real *charge;
- /* #endif */
- /* #if 'GeneralizedBorn' in KERNEL_ELEC */
- int gbitab;
- real vgb,fgb,vgbsum,dvdasum,gbscale,gbtabscale,isaprod,gbqqfactor,gbinvepsdiff,dvdaj,gbeps,dvdatmp;
- real *invsqrta,*dvda,*gbtab;
- /* #endif */
- /* #if KERNEL_VDW != 'None' */
- int nvdwtype;
- real rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,br,vvdwexp,sh_vdw_invrcut6;
- int *vdwtype;
- real *vdwparam;
- /* #endif */
- /* #if 'Table' in KERNEL_ELEC or 'GeneralizedBorn' in KERNEL_ELEC or 'Table' in KERNEL_VDW */
- int vfitab;
- real rt,vfeps,vftabscale,Y,F,Geps,Heps2,Fp,VV,FF;
- real *vftab;
- /* #endif */
- /* #if 'Ewald' in KERNEL_ELEC */
- int ewitab;
- real ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace;
- real *ewtab;
- /* #endif */
- /* #if 'PotentialSwitch' in [KERNEL_MOD_ELEC,KERNEL_MOD_VDW] */
- real rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
- /* #endif */
-
- real * wf;
- real weight_cg1;
- real weight_cg2;
- real weight_product;
- real hybscal; /* the multiplicator to the force for hybrid interactions*/
- gmx_bool bHybrid; /*Are we in the hybrid zone ?*/
- real force_cap;
-
- wf = mdatoms->wf;
-
- force_cap = fr->adress_ex_forcecap;
-
- x = xx[0];
- f = ff[0];
-
- nri = nlist->nri;
- iinr = nlist->iinr;
- jindex = nlist->jindex;
- jjnr = nlist->jjnr;
- shiftidx = nlist->shift;
- gid = nlist->gid;
- shiftvec = fr->shift_vec[0];
- fshift = fr->fshift[0];
- /* #if KERNEL_ELEC != 'None' */
- facel = fr->epsfac;
- charge = mdatoms->chargeA;
- /* #if 'ReactionField' in KERNEL_ELEC */
- krf = fr->ic->k_rf;
- krf2 = krf*2.0;
- crf = fr->ic->c_rf;
- /* #endif */
- /* #endif */
- /* #if KERNEL_VDW != 'None' */
- nvdwtype = fr->ntype;
- vdwparam = fr->nbfp;
- vdwtype = mdatoms->typeA;
- /* #endif */
-
- /* #if 'Table' in KERNEL_ELEC and 'Table' in KERNEL_VDW */
- vftab = kernel_data->table_elec_vdw->data;
- vftabscale = kernel_data->table_elec_vdw->scale;
- /* #elif 'Table' in KERNEL_ELEC */
- vftab = kernel_data->table_elec->data;
- vftabscale = kernel_data->table_elec->scale;
- /* #elif 'Table' in KERNEL_VDW */
- vftab = kernel_data->table_vdw->data;
- vftabscale = kernel_data->table_vdw->scale;
- /* #endif */
-
- /* #if 'Ewald' in KERNEL_ELEC */
- sh_ewald = fr->ic->sh_ewald;
- /* #if KERNEL_VF=='Force' and KERNEL_MOD_ELEC!='PotentialSwitch' */
- ewtab = fr->ic->tabq_coul_F;
- ewtabscale = fr->ic->tabq_scale;
- ewtabhalfspace = 0.5/ewtabscale;
- /* #else */
- ewtab = fr->ic->tabq_coul_FDV0;
- ewtabscale = fr->ic->tabq_scale;
- ewtabhalfspace = 0.5/ewtabscale;
- /* #endif */
- /* #endif */
-
- /* #if KERNEL_ELEC=='GeneralizedBorn' */
- invsqrta = fr->invsqrta;
- dvda = fr->dvda;
- gbtabscale = fr->gbtab.scale;
- gbtab = fr->gbtab.data;
- gbinvepsdiff = (1.0/fr->epsilon_r) - (1.0/fr->gb_epsilon_solvent);
- /* #endif */
-
- /* #if 'Water' in GEOMETRY_I */
- /* Setup water-specific parameters */
- inr = nlist->iinr[0];
- /* #for I in PARTICLES_ELEC_I */
- iq{I} = facel*charge[inr+{I}];
- /* #endfor */
- /* #for I in PARTICLES_VDW_I */
- vdwioffset{I} = {NVDWPARAM}*nvdwtype*vdwtype[inr+{I}];
- /* #endfor */
- /* #endif */
-
- /* #if 'Water' in GEOMETRY_J */
- /* #for J in PARTICLES_ELEC_J */
- jq{J} = charge[inr+{J}];
- /* #endfor */
- /* #for J in PARTICLES_VDW_J */
- vdwjidx{J} = {NVDWPARAM}*vdwtype[inr+{J}];
- /* #endfor */
- /* #for I,J in PAIRS_IJ */
- /* #if 'electrostatics' in INTERACTION_FLAGS[I][J] */
- qq{I}{J} = iq{I}*jq{J};
- /* #endif */
- /* #if 'vdw' in INTERACTION_FLAGS[I][J] */
- /* #if 'Buckingham' in KERNEL_VDW */
- c6_{I}{J} = vdwparam[vdwioffset{I}+vdwjidx{J}];
- cexp1_{I}{J} = vdwparam[vdwioffset{I}+vdwjidx{J}+1];
- cexp2_{I}{J} = vdwparam[vdwioffset{I}+vdwjidx{J}+2];
- /* #else */
- c6_{I}{J} = vdwparam[vdwioffset{I}+vdwjidx{J}];
- c12_{I}{J} = vdwparam[vdwioffset{I}+vdwjidx{J}+1];
- /* #endif */
- /* #endif */
- /* #endfor */
- /* #endif */
-
- /* #if KERNEL_MOD_ELEC!='None' or KERNEL_MOD_VDW!='None' */
- /* #if KERNEL_ELEC!='None' */
- /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
- rcutoff = fr->rcoulomb;
- /* #else */
- rcutoff = fr->rvdw;
- /* #endif */
- rcutoff2 = rcutoff*rcutoff;
- /* #endif */
-
- /* #if KERNEL_MOD_VDW=='PotentialShift' */
- sh_vdw_invrcut6 = fr->ic->sh_invrc6;
- rvdw = fr->rvdw;
- /* #endif */
-
- /* #if 'PotentialSwitch' in [KERNEL_MOD_ELEC,KERNEL_MOD_VDW] */
- /* #if KERNEL_MOD_ELEC=='PotentialSwitch' */
- rswitch = fr->rcoulomb_switch;
- /* #else */
- rswitch = fr->rvdw_switch;
- /* #endif */
- /* Setup switch parameters */
- d = rcutoff-rswitch;
- swV3 = -10.0/(d*d*d);
- swV4 = 15.0/(d*d*d*d);
- swV5 = -6.0/(d*d*d*d*d);
- /* #if 'Force' in KERNEL_VF */
- swF2 = -30.0/(d*d*d);
- swF3 = 60.0/(d*d*d*d);
- swF4 = -30.0/(d*d*d*d*d);
- /* #endif */
- /* #endif */
-
- /* ## Keep track of the floating point operations we issue for reporting! */
- /* #define OUTERFLOPS 0 */
- /* #define INNERFLOPS 0 */
- outeriter = 0;
- inneriter = 0;
-
- /* Start outer loop over neighborlists */
- for(iidx=0; iidx<nri; iidx++)
- {
- /* Load shift vector for this list */
- i_shift_offset = DIM*shiftidx[iidx];
- shX = shiftvec[i_shift_offset+XX];
- shY = shiftvec[i_shift_offset+YY];
- shZ = shiftvec[i_shift_offset+ZZ];
-
- /* Load limits for loop over neighbors */
- j_index_start = jindex[iidx];
- j_index_end = jindex[iidx+1];
-
- /* Get outer coordinate index */
- inr = iinr[iidx];
- i_coord_offset = DIM*inr;
-
- /* Load i particle coords and add shift vector */
- /* ## Loop over i particles, but only include ones that we use - skip e.g. vdw-only sites for elec-only kernel */
- /* #for I in PARTICLES_I */
- ix{I} = shX + x[i_coord_offset+DIM*{I}+XX];
- iy{I} = shY + x[i_coord_offset+DIM*{I}+YY];
- iz{I} = shZ + x[i_coord_offset+DIM*{I}+ZZ];
- /* #define OUTERFLOPS OUTERFLOPS+3 */
- /* #endfor */
-
- /* #if 'Force' in KERNEL_VF */
- /* #for I in PARTICLES_I */
- fix{I} = 0.0;
- fiy{I} = 0.0;
- fiz{I} = 0.0;
- /* #endfor */
- /* #endif */
-
- weight_cg1 = wf[inr];
-
- /* ## For water we already preloaded parameters at the start of the kernel */
- /* #if not 'Water' in GEOMETRY_I */
- /* Load parameters for i particles */
- /* #for I in PARTICLES_ELEC_I */
- iq{I} = facel*charge[inr+{I}];
- /* #define OUTERFLOPS OUTERFLOPS+1 */
- /* #if KERNEL_ELEC=='GeneralizedBorn' */
- isai{I} = invsqrta[inr+{I}];
- /* #endif */
- /* #endfor */
- /* #for I in PARTICLES_VDW_I */
- vdwioffset{I} = {NVDWPARAM}*nvdwtype*vdwtype[inr+{I}];
- /* #endfor */
- /* #endif */
-
- /* #if 'Potential' in KERNEL_VF */
- /* Reset potential sums */
- /* #if KERNEL_ELEC != 'None' */
- velecsum = 0.0;
- /* #endif */
- /* #if 'GeneralizedBorn' in KERNEL_ELEC */
- vgbsum = 0.0;
- /* #endif */
- /* #if KERNEL_VDW != 'None' */
- vvdwsum = 0.0;
- /* #endif */
- /* #endif */
- /* #if 'GeneralizedBorn' in KERNEL_ELEC and 'Force' in KERNEL_VF */
- dvdasum = 0.0;
- /* #endif */
-
- /* Start inner kernel loop */
- for(jidx=j_index_start; jidx<j_index_end; jidx++)
- {
- /* Get j neighbor index, and coordinate index */
- jnr = jjnr[jidx];
- weight_cg2 = wf[jnr];
- weight_product = weight_cg1*weight_cg2;
- if (weight_product < ALMOST_ZERO) {
- /* #if KERNEL_RESOLUTION=='CG' */
- hybscal = 1.0;
- /* #else */
- continue;
- /* #endif */
- }
- else if (weight_product >= ALMOST_ONE)
- {
- /* #if KERNEL_RESOLUTION=='CG' */
- continue;
- /* #else */
- hybscal = 1.0;
- /* #endif */
- }
- else
- {
- /* #if KERNEL_RESOLUTION=='CG' */
- hybscal = 1.0 - weight_product;
- /* #else */
- hybscal = weight_product;
- /* #endif */
- }
- j_coord_offset = DIM*jnr;
-
- /* load j atom coordinates */
- /* #for J in PARTICLES_J */
- jx{J} = x[j_coord_offset+DIM*{J}+XX];
- jy{J} = x[j_coord_offset+DIM*{J}+YY];
- jz{J} = x[j_coord_offset+DIM*{J}+ZZ];
- /* #endfor */
-
- /* Calculate displacement vector */
- /* #for I,J in PAIRS_IJ */
- dx{I}{J} = ix{I} - jx{J};
- dy{I}{J} = iy{I} - jy{J};
- dz{I}{J} = iz{I} - jz{J};
- /* #define INNERFLOPS INNERFLOPS+3 */
- /* #endfor */
-
- /* Calculate squared distance and things based on it */
- /* #for I,J in PAIRS_IJ */
- rsq{I}{J} = dx{I}{J}*dx{I}{J}+dy{I}{J}*dy{I}{J}+dz{I}{J}*dz{I}{J};
- /* #define INNERFLOPS INNERFLOPS+5 */
- /* #endfor */
-
- /* #for I,J in PAIRS_IJ */
- /* #if 'rinv' in INTERACTION_FLAGS[I][J] */
- rinv{I}{J} = gmx_invsqrt(rsq{I}{J});
- /* #define INNERFLOPS INNERFLOPS+5 */
- /* #endif */
- /* #endfor */
-
- /* #for I,J in PAIRS_IJ */
- /* #if 'rinvsq' in INTERACTION_FLAGS[I][J] */
- /* # if 'rinv' not in INTERACTION_FLAGS[I][J] */
- rinvsq{I}{J} = 1.0/rsq{I}{J};
- /* #define INNERFLOPS INNERFLOPS+4 */
- /* #else */
- rinvsq{I}{J} = rinv{I}{J}*rinv{I}{J};
- /* #define INNERFLOPS INNERFLOPS+1 */
- /* #endif */
- /* #endif */
- /* #endfor */
-
- /* #if not 'Water' in GEOMETRY_J */
- /* Load parameters for j particles */
- /* #for J in PARTICLES_ELEC_J */
- jq{J} = charge[jnr+{J}];
- /* #if KERNEL_ELEC=='GeneralizedBorn' */
- isaj{J} = invsqrta[jnr+{J}];
- /* #endif */
- /* #endfor */
- /* #for J in PARTICLES_VDW_J */
- vdwjidx{J} = {NVDWPARAM}*vdwtype[jnr+{J}];
- /* #endfor */
- /* #endif */
-
- /* #for I,J in PAIRS_IJ */
-
- /**************************
- * CALCULATE INTERACTIONS *
- **************************/
-
- /* ## Note special check for TIP4P-TIP4P. Since we are cutting of all hydrogen interactions we also cut the LJ-only O-O interaction */
- /* #if 'exactcutoff' in INTERACTION_FLAGS[I][J] or (GEOMETRY_I=='Water4' and GEOMETRY_J=='Water4' and 'exactcutoff' in INTERACTION_FLAGS[1][1]) */
- /* ## We always calculate rinv/rinvsq above to enable pipelineing in compilers (performance tested on x86) */
- if (rsq{I}{J}<rcutoff2)
- {
- /* #if 0 ## this and the next two lines is a hack to maintain auto-indentation in template file */
- }
- /* #endif */
- /* #endif */
-
- /* #if 'r' in INTERACTION_FLAGS[I][J] */
- r{I}{J} = rsq{I}{J}*rinv{I}{J};
- /* #define INNERFLOPS INNERFLOPS+1 */
- /* #endif */
-
- /* ## For water geometries we already loaded parameters at the start of the kernel */
- /* #if not 'Water' in GEOMETRY_J */
- /* #if 'electrostatics' in INTERACTION_FLAGS[I][J] */
- qq{I}{J} = iq{I}*jq{J};
- /* #define INNERFLOPS INNERFLOPS+1 */
- /* #endif */
- /* #if 'vdw' in INTERACTION_FLAGS[I][J] */
- /* #if KERNEL_VDW=='Buckingham' */
- c6_{I}{J} = vdwparam[vdwioffset{I}+vdwjidx{J}];
- cexp1_{I}{J} = vdwparam[vdwioffset{I}+vdwjidx{J}+1];
- cexp2_{I}{J} = vdwparam[vdwioffset{I}+vdwjidx{J}+2];
- /* #else */
- c6_{I}{J} = vdwparam[vdwioffset{I}+vdwjidx{J}];
- c12_{I}{J} = vdwparam[vdwioffset{I}+vdwjidx{J}+1];
- /* #endif */
- /* #endif */
- /* #endif */
-
- /* #if 'table' in INTERACTION_FLAGS[I][J] */
- /* Calculate table index by multiplying r with table scale and truncate to integer */
- rt = r{I}{J}*vftabscale;
- vfitab = rt;
- vfeps = rt-vfitab;
- vfitab = {TABLE_NINTERACTIONS}*{TABLE_POINT_SIZE}*vfitab;
- /* #define INNERFLOPS INNERFLOPS+2 */
- /* #endif */
-
- /* ## ELECTROSTATIC INTERACTIONS */
- /* #if 'electrostatics' in INTERACTION_FLAGS[I][J] */
-
- /* #if KERNEL_ELEC=='Coulomb' */
-
- /* COULOMB ELECTROSTATICS */
- velec = qq{I}{J}*rinv{I}{J};
- /* #define INNERFLOPS INNERFLOPS+1 */
- /* #if 'Force' in KERNEL_VF */
- felec = velec*rinvsq{I}{J};
- /* #define INNERFLOPS INNERFLOPS+2 */
- /* #endif */
-
- /* #elif KERNEL_ELEC=='ReactionField' */
-
- /* REACTION-FIELD ELECTROSTATICS */
- /* #if 'Potential' in KERNEL_VF */
- velec = qq{I}{J}*(rinv{I}{J}+krf*rsq{I}{J}-crf);
- /* #define INNERFLOPS INNERFLOPS+4 */
- /* #endif */
- /* #if 'Force' in KERNEL_VF */
- felec = qq{I}{J}*(rinv{I}{J}*rinvsq{I}{J}-krf2);
- /* #define INNERFLOPS INNERFLOPS+3 */
- /* #endif */
-
- /* #elif KERNEL_ELEC=='GeneralizedBorn' */
-
- /* GENERALIZED BORN AND COULOMB ELECTROSTATICS */
- isaprod = isai{I}*isaj{J};
- gbqqfactor = isaprod*(-qq{I}{J})*gbinvepsdiff;
- gbscale = isaprod*gbtabscale;
- dvdaj = dvda[jnr+{J}];
- /* #define INNERFLOPS INNERFLOPS+5 */
-
- /* Calculate generalized born table index - this is a separate table from the normal one,
- * but we use the same procedure by multiplying r with scale and truncating to integer.
- */
- rt = r{I}{J}*gbscale;
- gbitab = rt;
- gbeps = rt-gbitab;
- gbitab = 4*gbitab;
-
- Y = gbtab[gbitab];
- F = gbtab[gbitab+1];
- Geps = gbeps*gbtab[gbitab+2];
- Heps2 = gbeps*gbeps*gbtab[gbitab+3];
- Fp = F+Geps+Heps2;
- VV = Y+gbeps*Fp;
- vgb = gbqqfactor*VV;
- /* #define INNERFLOPS INNERFLOPS+10 */
-
- /* #if 'Force' in KERNEL_VF */
- FF = Fp+Geps+2.0*Heps2;
- fgb = gbqqfactor*FF*gbscale;
- dvdatmp = -0.5*(vgb+fgb*r{I}{J});
- dvdasum = dvdasum + dvdatmp;
- dvda[jnr] = dvdaj+dvdatmp*isaj{J}*isaj{J};
- /* #define INNERFLOPS INNERFLOPS+13 */
- /* #endif */
- velec = qq{I}{J}*rinv{I}{J};
- /* #define INNERFLOPS INNERFLOPS+1 */
- /* #if 'Force' in KERNEL_VF */
- felec = (velec*rinv{I}{J}-fgb)*rinv{I}{J};
- /* #define INNERFLOPS INNERFLOPS+3 */
- /* #endif */
-
- /* #elif KERNEL_ELEC=='Ewald' */
- /* EWALD ELECTROSTATICS */
-
- /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
- ewrt = r{I}{J}*ewtabscale;
- ewitab = ewrt;
- eweps = ewrt-ewitab;
- /* #define INNERFLOPS INNERFLOPS+2 */
- /* #if 'Potential' in KERNEL_VF or KERNEL_MOD_ELEC=='PotentialSwitch' */
- ewitab = 4*ewitab;
- felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
- /* #define INNERFLOPS INNERFLOPS+4 */
- /* #if KERNEL_MOD_ELEC=='PotentialShift' */
- velec = qq{I}{J}*((rinv{I}{J}-sh_ewald)-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
- /* #define INNERFLOPS INNERFLOPS+7 */
- /* #else */
- velec = qq{I}{J}*(rinv{I}{J}-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
- /* #define INNERFLOPS INNERFLOPS+6 */
- /* #endif */
- /* #if 'Force' in KERNEL_VF */
- felec = qq{I}{J}*rinv{I}{J}*(rinvsq{I}{J}-felec);
- /* #define INNERFLOPS INNERFLOPS+3 */
- /* #endif */
- /* #elif KERNEL_VF=='Force' */
- felec = (1.0-eweps)*ewtab[ewitab]+eweps*ewtab[ewitab+1];
- felec = qq{I}{J}*rinv{I}{J}*(rinvsq{I}{J}-felec);
- /* #define INNERFLOPS INNERFLOPS+7 */
- /* #endif */
-
- /* #elif KERNEL_ELEC=='CubicSplineTable' */
-
- /* CUBIC SPLINE TABLE ELECTROSTATICS */
- /* #if 'Potential' in KERNEL_VF */
- Y = vftab[vfitab];
- /* #endif */
- F = vftab[vfitab+1];
- Geps = vfeps*vftab[vfitab+2];
- Heps2 = vfeps*vfeps*vftab[vfitab+3];
- Fp = F+Geps+Heps2;
- /* #define INNERFLOPS INNERFLOPS+5 */
- /* #if 'Potential' in KERNEL_VF */
- VV = Y+vfeps*Fp;
- velec = qq{I}{J}*VV;
- /* #define INNERFLOPS INNERFLOPS+3 */
- /* #endif */
- /* #if 'Force' in KERNEL_VF */
- FF = Fp+Geps+2.0*Heps2;
- felec = -qq{I}{J}*FF*vftabscale*rinv{I}{J};
- /* #define INNERFLOPS INNERFLOPS+7 */
- /* #endif */
- /* #endif */
- /* ## End of check for electrostatics interaction forms */
- /* #endif */
- /* ## END OF ELECTROSTATIC INTERACTION CHECK FOR PAIR I-J */
-
- /* #if 'vdw' in INTERACTION_FLAGS[I][J] */
-
- /* #if KERNEL_VDW=='LennardJones' */
-
- /* LENNARD-JONES DISPERSION/REPULSION */
-
- rinvsix = rinvsq{I}{J}*rinvsq{I}{J}*rinvsq{I}{J};
- /* #define INNERFLOPS INNERFLOPS+2 */
- /* #if 'Potential' in KERNEL_VF or KERNEL_MOD_VDW=='PotentialSwitch' */
- vvdw6 = c6_{I}{J}*rinvsix;
- vvdw12 = c12_{I}{J}*rinvsix*rinvsix;
- /* #define INNERFLOPS INNERFLOPS+3 */
- /* #if KERNEL_MOD_VDW=='PotentialShift' */
- vvdw = (vvdw12 - c12_{I}{J}*sh_vdw_invrcut6*sh_vdw_invrcut6)*(1.0/12.0) - (vvdw6 - c6_{I}{J}*sh_vdw_invrcut6)*(1.0/6.0);
- /* #define INNERFLOPS INNERFLOPS+8 */
- /* #else */
- vvdw = vvdw12*(1.0/12.0) - vvdw6*(1.0/6.0);
- /* #define INNERFLOPS INNERFLOPS+3 */
- /* #endif */
- /* ## Check for force inside potential check, i.e. this means we already did the potential part */
- /* #if 'Force' in KERNEL_VF */
- fvdw = (vvdw12-vvdw6)*rinvsq{I}{J};
- /* #define INNERFLOPS INNERFLOPS+2 */
- /* #endif */
- /* #elif KERNEL_VF=='Force' */
- /* ## Force-only LennardJones makes it possible to save 1 flop (they do add up...) */
- fvdw = (c12_{I}{J}*rinvsix-c6_{I}{J})*rinvsix*rinvsq{I}{J};
- /* #define INNERFLOPS INNERFLOPS+4 */
- /* #endif */
-
- /* #elif KERNEL_VDW=='Buckingham' */
-
- /* BUCKINGHAM DISPERSION/REPULSION */
- rinvsix = rinvsq{I}{J}*rinvsq{I}{J}*rinvsq{I}{J};
- vvdw6 = c6_{I}{J}*rinvsix;
- br = cexp2_{I}{J}*r{I}{J};
- vvdwexp = cexp1_{I}{J}*exp(-br);
- /* ## Estimate exp() to 25 flops */
- /* #define INNERFLOPS INNERFLOPS+31 */
- /* #if 'Potential' in KERNEL_VF or KERNEL_MOD_VDW=='PotentialSwitch' */
- /* #if KERNEL_MOD_VDW=='PotentialShift' */
- vvdw = (vvdwexp-cexp1_{I}{J}*exp(-cexp2_{I}{J}*rvdw)) - (vvdw6 - c6_{I}{J}*sh_vdw_invrcut6)*(1.0/6.0);
- /* #define INNERFLOPS INNERFLOPS+33 */
- /* #else */
- vvdw = vvdwexp - vvdw6*(1.0/6.0);
- /* #define INNERFLOPS INNERFLOPS+2 */
- /* #endif */
- /* #endif */
- /* #if 'Force' in KERNEL_VF */
- fvdw = (br*vvdwexp-vvdw6)*rinvsq{I}{J};
- /* #define INNERFLOPS INNERFLOPS+3 */
- /* #endif */
-
- /* #elif KERNEL_VDW=='CubicSplineTable' */
-
- /* CUBIC SPLINE TABLE DISPERSION */
- vfitab += {TABLE_VDW_OFFSET};
- /* #if 'Potential' in KERNEL_VF */
- Y = vftab[vfitab];
- /* #endif */
- F = vftab[vfitab+1];
- Geps = vfeps*vftab[vfitab+2];
- Heps2 = vfeps*vfeps*vftab[vfitab+3];
- Fp = F+Geps+Heps2;
- /* #define INNERFLOPS INNERFLOPS+5 */
- /* #if 'Potential' in KERNEL_VF */
- VV = Y+vfeps*Fp;
- vvdw6 = c6_{I}{J}*VV;
- /* #define INNERFLOPS INNERFLOPS+3 */
- /* #endif */
- /* #if 'Force' in KERNEL_VF */
- FF = Fp+Geps+2.0*Heps2;
- fvdw6 = c6_{I}{J}*FF;
- /* #define INNERFLOPS INNERFLOPS+4 */
- /* #endif */
-
- /* CUBIC SPLINE TABLE REPULSION */
- /* #if 'Potential' in KERNEL_VF */
- Y = vftab[vfitab+4];
- /* #endif */
- F = vftab[vfitab+5];
- Geps = vfeps*vftab[vfitab+6];
- Heps2 = vfeps*vfeps*vftab[vfitab+7];
- Fp = F+Geps+Heps2;
- /* #define INNERFLOPS INNERFLOPS+5 */
- /* #if 'Potential' in KERNEL_VF */
- VV = Y+vfeps*Fp;
- vvdw12 = c12_{I}{J}*VV;
- /* #define INNERFLOPS INNERFLOPS+3 */
- /* #endif */
- /* #if 'Force' in KERNEL_VF */
- FF = Fp+Geps+2.0*Heps2;
- fvdw12 = c12_{I}{J}*FF;
- /* #define INNERFLOPS INNERFLOPS+4 */
- /* #endif */
- /* #if 'Potential' in KERNEL_VF */
- vvdw = vvdw12+vvdw6;
- /* #define INNERFLOPS INNERFLOPS+1 */
- /* #endif */
- /* #if 'Force' in KERNEL_VF */
- fvdw = -(fvdw6+fvdw12)*vftabscale*rinv{I}{J};
- /* #define INNERFLOPS INNERFLOPS+4 */
- /* #endif */
- /* #endif */
- /* ## End of check for vdw interaction forms */
- /* #endif */
- /* ## END OF VDW INTERACTION CHECK FOR PAIR I-J */
-
- /* #if 'switch' in INTERACTION_FLAGS[I][J] */
- d = r{I}{J}-rswitch;
- d = (d>0.0) ? d : 0.0;
- d2 = d*d;
- sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
- /* #define INNERFLOPS INNERFLOPS+9 */
-
- /* #if 'Force' in KERNEL_VF */
- dsw = d2*(swF2+d*(swF3+d*swF4));
- /* #define INNERFLOPS INNERFLOPS+5 */
- /* #endif */
-
- /* Evaluate switch function */
- /* #if 'Force' in KERNEL_VF */
- /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
- /* #if 'electrostatics' in INTERACTION_FLAGS[I][J] and KERNEL_MOD_ELEC=='PotentialSwitch' */
- felec = felec*sw - rinv{I}{J}*velec*dsw;
- /* #define INNERFLOPS INNERFLOPS+3 */
- /* #endif */
- /* #if 'vdw' in INTERACTION_FLAGS[I][J] and KERNEL_MOD_VDW=='PotentialSwitch' */
- fvdw = fvdw*sw - rinv{I}{J}*vvdw*dsw;
- /* #define INNERFLOPS INNERFLOPS+3 */
- /* #endif */
- /* #endif */
- /* #if 'Potential' in KERNEL_VF */
- /* #if 'electrostatics' in INTERACTION_FLAGS[I][J] and KERNEL_MOD_ELEC=='PotentialSwitch' */
- velec *= sw;
- /* #define INNERFLOPS INNERFLOPS+1 */
- /* #endif */
- /* #if 'vdw' in INTERACTION_FLAGS[I][J] and KERNEL_MOD_VDW=='PotentialSwitch' */
- vvdw *= sw;
- /* #define INNERFLOPS INNERFLOPS+1 */
- /* #endif */
- /* #endif */
- /* #endif */
-
- /* #if 'Potential' in KERNEL_VF */
- /* Update potential sums from outer loop */
- /* #if 'electrostatics' in INTERACTION_FLAGS[I][J] */
- velecsum += velec;
- /* #define INNERFLOPS INNERFLOPS+1 */
- /* #if KERNEL_ELEC=='GeneralizedBorn' */
- vgbsum += vgb;
- /* #define INNERFLOPS INNERFLOPS+1 */
- /* #endif */
- /* #endif */
- /* #if 'vdw' in INTERACTION_FLAGS[I][J] */
- vvdwsum += vvdw;
- /* #define INNERFLOPS INNERFLOPS+1 */
- /* #endif */
- /* #endif */
-
- /* #if 'Force' in KERNEL_VF */
-
- /* #if 'electrostatics' in INTERACTION_FLAGS[I][J] and 'vdw' in INTERACTION_FLAGS[I][J] */
- fscal = felec+fvdw;
- /* #define INNERFLOPS INNERFLOPS+1 */
- /* #elif 'electrostatics' in INTERACTION_FLAGS[I][J] */
- fscal = felec;
- /* #elif 'vdw' in INTERACTION_FLAGS[I][J] */
- fscal = fvdw;
- /* #endif */
-
- /* #if KERNEL_RESOLUTION=='CG' */
- if(force_cap>0 && (fabs(fscal)> force_cap))
- {
- fscal=force_cap*fscal/fabs(fscal);
- }
- /* #endif */
- fscal *= hybscal;
-
- /* Calculate temporary vectorial force */
- tx = fscal*dx{I}{J};
- ty = fscal*dy{I}{J};
- tz = fscal*dz{I}{J};
-
- /* Update vectorial force */
- fix{I} += tx;
- fiy{I} += ty;
- fiz{I} += tz;
- f[j_coord_offset+DIM*{J}+XX] -= tx;
- f[j_coord_offset+DIM*{J}+YY] -= ty;
- f[j_coord_offset+DIM*{J}+ZZ] -= tz;
-
- /* #define INNERFLOPS INNERFLOPS+9 */
- /* #endif */
-
- /* ## Note special check for TIP4P-TIP4P. Since we are cutting of all hydrogen interactions we also cut the LJ-only O-O interaction */
- /* #if 'exactcutoff' in INTERACTION_FLAGS[I][J] or (GEOMETRY_I=='Water4' and GEOMETRY_J=='Water4' and 'exactcutoff' in INTERACTION_FLAGS[1][1]) */
- /* #if 0 ## This and next two lines is a hack to maintain indentation in template file */
- {
- /* #endif */
- }
- /* #endif */
- /* ## End of check for the interaction being outside the cutoff */
-
- /* #endfor */
- /* ## End of loop over i-j interaction pairs */
-
- /* Inner loop uses {INNERFLOPS} flops */
- }
- /* End of innermost loop */
-
- /* #if 'Force' in KERNEL_VF */
- tx = ty = tz = 0;
- /* #for I in PARTICLES_I */
- f[i_coord_offset+DIM*{I}+XX] += fix{I};
- f[i_coord_offset+DIM*{I}+YY] += fiy{I};
- f[i_coord_offset+DIM*{I}+ZZ] += fiz{I};
- tx += fix{I};
- ty += fiy{I};
- tz += fiz{I};
- /* #define OUTERFLOPS OUTERFLOPS+6 */
- /* #endfor */
- fshift[i_shift_offset+XX] += tx;
- fshift[i_shift_offset+YY] += ty;
- fshift[i_shift_offset+ZZ] += tz;
- /* #define OUTERFLOPS OUTERFLOPS+3 */
- /* #endif */
-
- /* #if 'Potential' in KERNEL_VF */
- ggid = gid[iidx];
- /* Update potential energies */
- /* #if KERNEL_ELEC != 'None' */
- kernel_data->energygrp_elec[ggid] += velecsum;
- /* #define OUTERFLOPS OUTERFLOPS+1 */
- /* #endif */
- /* #if 'GeneralizedBorn' in KERNEL_ELEC */
- kernel_data->energygrp_polarization[ggid] += vgbsum;
- /* #define OUTERFLOPS OUTERFLOPS+1 */
- /* #endif */
- /* #if KERNEL_VDW != 'None' */
- kernel_data->energygrp_vdw[ggid] += vvdwsum;
- /* #define OUTERFLOPS OUTERFLOPS+1 */
- /* #endif */
- /* #endif */
- /* #if 'GeneralizedBorn' in KERNEL_ELEC and 'Force' in KERNEL_VF */
- dvda[inr] = dvda[inr] + dvdasum*isai{I}*isai{I};
- /* #endif */
-
- /* Increment number of inner iterations */
- inneriter += j_index_end - j_index_start;
-
- /* Outer loop uses {OUTERFLOPS} flops */
- }
-
- /* Increment number of outer iterations */
- outeriter += nri;
-
- /* Update outer/inner flops */
- /* ## NB: This is not important, it just affects the flopcount. However, since our preprocessor is */
- /* ## primitive and replaces aggressively even in strings inside these directives, we need to */
- /* ## assemble the main part of the name (containing KERNEL/ELEC/VDW) directly in the source. */
- /* #if GEOMETRY_I == 'Water3' */
- /* #define ISUFFIX '_W3' */
- /* #elif GEOMETRY_I == 'Water4' */
- /* #define ISUFFIX '_W4' */
- /* #else */
- /* #define ISUFFIX '' */
- /* #endif */
- /* #if GEOMETRY_J == 'Water3' */
- /* #define JSUFFIX 'W3' */
- /* #elif GEOMETRY_J == 'Water4' */
- /* #define JSUFFIX 'W4' */
- /* #else */
- /* #define JSUFFIX '' */
- /* #endif */
- /* #if 'PotentialAndForce' in KERNEL_VF */
- /* #define VFSUFFIX '_VF' */
- /* #elif 'Potential' in KERNEL_VF */
- /* #define VFSUFFIX '_V' */
- /* #else */
- /* #define VFSUFFIX '_F' */
- /* #endif */
-
- /* #if KERNEL_ELEC != 'None' and KERNEL_VDW != 'None' */
- inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW{ISUFFIX}{JSUFFIX}{VFSUFFIX},outeriter*{OUTERFLOPS} + inneriter*{INNERFLOPS});
- /* #elif KERNEL_ELEC != 'None' */
- inc_nrnb(nrnb,eNR_NBKERNEL_ELEC{ISUFFIX}{JSUFFIX}{VFSUFFIX},outeriter*{OUTERFLOPS} + inneriter*{INNERFLOPS});
- /* #else */
- inc_nrnb(nrnb,eNR_NBKERNEL_VDW{ISUFFIX}{JSUFFIX}{VFSUFFIX},outeriter*{OUTERFLOPS} + inneriter*{INNERFLOPS});
- /* #endif */
-}
#include "gromacs/fileio/txtdump.h"
#include "gromacs/gmxlib/nonbonded/nb_free_energy.h"
#include "gromacs/gmxlib/nonbonded/nb_generic.h"
-#include "gromacs/gmxlib/nonbonded/nb_generic_adress.h"
#include "gromacs/gmxlib/nonbonded/nb_generic_cg.h"
#include "gromacs/gmxlib/nonbonded/nb_kernel.h"
#include "gromacs/legacyheaders/names.h"
vdw_mod = eintmod_names[nl->ivdwmod];
geom = gmx_nblist_geometry_names[nl->igeometry];
- if (nl->type == GMX_NBLIST_INTERACTION_ADRESS)
- {
- nl->kernelptr_vf = (void *) gmx_nb_generic_adress_kernel;
- nl->kernelptr_f = (void *) gmx_nb_generic_adress_kernel;
- nl->simd_padding_width = 1;
- return;
- }
-
if (nl->type == GMX_NBLIST_INTERACTION_FREE_ENERGY)
{
nl->kernelptr_vf = (void *) gmx_nb_free_energy_kernel;
{ "NB Generic kernel", 1 },
{ "NB Generic charge grp kernel", 1 },
- { "NB Generic AdResS kernel", 1 },
{ "NB Free energy kernel", 1 },
{ "NB All-vs-all", 1 },
{ "NB All-vs-all, GB", 1 },
"NB: United atoms have the same atom numbers as normal ones.\n\n");
}
- if (ir->bAdress)
- {
- if ((ir->adress->const_wf > 1) || (ir->adress->const_wf < 0))
- {
- warning_error(wi, "AdResS contant weighting function should be between 0 and 1\n\n");
- }
- /** TODO check size of ex+hy width against box size */
- }
-
/* Check for errors in the input now, since they might cause problems
* during processing further down.
*/
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2009 Christoph Junghans, Brad Lambeth.
- * Copyright (c) 2012,2014,2015, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include <stdlib.h>
-#include <string.h>
-
-#include "gromacs/gmxpreprocess/readir.h"
-#include "gromacs/legacyheaders/names.h"
-#include "gromacs/mdtypes/inputrec.h"
-#include "gromacs/utility/cstringutil.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/smalloc.h"
-
-#define MAXPTR 254
-
-static char adress_refs[STRLEN], adress_tf_grp_names[STRLEN], adress_cg_grp_names[STRLEN];
-
-void read_adressparams(int *ninp_p, t_inpfile **inp_p, t_adress *adress, warninp_t wi)
-{
- int nadress_refs, i;
- const char *tmp;
- char *ptr1[MAXPTR];
-
-
- int ninp;
- t_inpfile *inp;
-
- ninp = *ninp_p;
- inp = *inp_p;
-
- EETYPE("adress_type", adress->type, eAdresstype_names);
- RTYPE ("adress_const_wf", adress->const_wf, 1);
- RTYPE ("adress_ex_width", adress->ex_width, 0);
- RTYPE ("adress_hy_width", adress->hy_width, 0);
- RTYPE ("adress_ex_forcecap", adress->ex_forcecap, 0);
- EETYPE("adress_interface_correction", adress->icor, eAdressICtype_names);
- EETYPE("adress_site", adress->site, eAdressSITEtype_names);
- STYPE ("adress_reference_coords", adress_refs, NULL);
- STYPE ("adress_tf_grp_names", adress_tf_grp_names, NULL);
- STYPE ("adress_cg_grp_names", adress_cg_grp_names, NULL);
- EETYPE("adress_do_hybridpairs", adress->do_hybridpairs, yesno_names);
-
- nadress_refs = str_nelem(adress_refs, MAXPTR, ptr1);
-
- for (i = 0; (i < nadress_refs); i++) /*read vector components*/
- {
- adress->refs[i] = strtod(ptr1[i], NULL);
- }
- for (; (i < DIM); i++) /*remaining undefined components of the vector set to zero*/
- {
- adress->refs[i] = 0;
- }
-
- *ninp_p = ninp;
- *inp_p = inp;
-}
-
-void do_adress_index(t_adress *adress, gmx_groups_t *groups, char **gnames, t_grpopts *opts, warninp_t wi)
-{
- int nr, i, j, k;
- char *ptr1[MAXPTR];
- int nadress_cg_grp_names, nadress_tf_grp_names;
-
- /* AdResS coarse grained groups input */
-
- nr = groups->grps[egcENER].nr;
- snew(adress->group_explicit, nr);
- for (i = 0; i < nr; i++)
- {
- adress->group_explicit[i] = TRUE;
- }
- adress->n_energy_grps = nr;
-
- nadress_cg_grp_names = str_nelem(adress_cg_grp_names, MAXPTR, ptr1);
-
- if (nadress_cg_grp_names > 0)
- {
- for (i = 0; i < nadress_cg_grp_names; i++)
- {
- /* search for the group name mathching the tf group name */
- k = 0;
- while ((k < nr) &&
- gmx_strcasecmp(ptr1[i], (char*)(gnames[groups->grps[egcENER].nm_ind[k]])))
- {
- k++;
- }
- if (k == nr)
- {
- gmx_fatal(FARGS, "Adress cg energy group %s not found\n", ptr1[i]);
- }
- adress->group_explicit[k] = FALSE;
- printf ("AdResS: Energy group %s is treated as coarse-grained \n",
- (char*)(gnames[groups->grps[egcENER].nm_ind[k]]));
- }
- /* set energy group exclusions between all coarse-grained and explicit groups */
- for (j = 0; j < nr; j++)
- {
- for (k = 0; k < nr; k++)
- {
- if (!(adress->group_explicit[k] == adress->group_explicit[j]))
- {
- opts->egp_flags[nr * j + k] |= EGP_EXCL;
- if (debug)
- {
- fprintf(debug, "AdResS excl %s %s \n",
- (char*)(gnames[groups->grps[egcENER].nm_ind[j]]),
- (char*)(gnames[groups->grps[egcENER].nm_ind[k]]));
- }
- }
- }
- }
- }
- else
- {
- warning(wi, "For an AdResS simulation at least one coarse-grained energy group has to be specified in adress_cg_grp_names");
- }
-
-
- /* AdResS multiple tf tables input */
- nadress_tf_grp_names = str_nelem(adress_tf_grp_names, MAXPTR, ptr1);
- adress->n_tf_grps = nadress_tf_grp_names;
- snew(adress->tf_table_index, nadress_tf_grp_names);
-
- nr = groups->grps[egcENER].nr;
-
- if (nadress_tf_grp_names > 0)
- {
- for (i = 0; i < nadress_tf_grp_names; i++)
- {
- /* search for the group name mathching the tf group name */
- k = 0;
- while ((k < nr) &&
- gmx_strcasecmp(ptr1[i], (char*)(gnames[groups->grps[egcENER].nm_ind[k]])))
- {
- k++;
- }
- if (k == nr)
- {
- gmx_fatal(FARGS, "Adress tf energy group %s not found\n", ptr1[i]);
- }
-
- adress->tf_table_index[i] = k;
- if (debug)
- {
- fprintf(debug, "found tf group %s id %d \n", ptr1[i], k);
- }
- if (adress->group_explicit[k])
- {
- gmx_fatal(FARGS, "Thermodynamic force group %s is not a coarse-grained group in adress_cg_grp_names. The thermodynamic force has to act on the coarse-grained vsite of a molecule.\n", ptr1[i]);
- }
-
- }
- }
- /* end AdResS multiple tf tables input */
-
-
-}
{
if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
{
- sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s a\
-nd %s",
+ sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
evdw_names[ir->vdwtype],
eintmod_names[eintmodPOTSHIFT],
eintmod_names[eintmodNONE]);
if (ir->bAdress)
{
- if (ir->cutoff_scheme != ecutsGROUP)
- {
- warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
- }
- if (!EI_SD(ir->eI))
- {
- warning_error(wi, "AdresS simulation supports only stochastic dynamics");
- }
- if (ir->epc != epcNO)
- {
- warning_error(wi, "AdresS simulation does not support pressure coupling");
- }
- if (EEL_FULL(ir->coulombtype))
- {
- warning_error(wi, "AdresS simulation does not support long-range electrostatics");
- }
+ gmx_fatal(FARGS, "AdResS simulations are no longer supported");
}
}
REM_TYPE("nstdihreout");
REM_TYPE("nstcheckpoint");
REM_TYPE("optimize-fft");
+ REM_TYPE("adress_type");
+ REM_TYPE("adress_const_wf");
+ REM_TYPE("adress_ex_width");
+ REM_TYPE("adress_hy_width");
+ REM_TYPE("adress_ex_forcecap");
+ REM_TYPE("adress_interface_correction");
+ REM_TYPE("adress_site");
+ REM_TYPE("adress_reference_coords");
+ REM_TYPE("adress_tf_grp_names");
+ REM_TYPE("adress_cg_grp_names");
+ REM_TYPE("adress_do_hybridpairs");
/* replace the following commands with the clearer new versions*/
REPL_TYPE("unconstrained-start", "continuation");
RTYPE("threshold", ir->swap->threshold, 1.0);
}
- /* AdResS defined thingies */
- CCTYPE ("AdResS parameters");
- EETYPE("adress", ir->bAdress, yesno_names);
- if (ir->bAdress)
- {
- snew(ir->adress, 1);
- read_adressparams(&ninp, &inp, ir->adress, wi);
- }
+ /* AdResS is no longer supported, but we need mdrun to be able to refuse to run old AdResS .tpr files */
+ EETYPE("adress", ir->bAdress, yesno_names);
/* User defined thingies */
CCTYPE ("User defined thingies");
decode_cos(is->efield_z, &(ir->ex[ZZ]));
decode_cos(is->efield_zt, &(ir->et[ZZ]));
- if (ir->bAdress)
- {
- do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
- }
-
for (i = 0; (i < grps->nr); i++)
{
sfree(gnames[i]);
struct gmx_mtop_t;
struct gmx_output_env_t;
struct pull_params_t;
-struct t_adress;
struct t_grpopts;
struct t_inputrec;
struct t_rot;
int str_nelem(const char *str, int maxptr, char *ptr[]);
/* helper function from readir.c to convert strings */
-void read_adressparams(int *ninp_p, t_inpfile **inp_p, t_adress *adress, warninp_t wi);
-/* Reads in AdResS related parameters */
-
-void do_adress_index(t_adress *adress, gmx_groups_t *groups, char **gnames, t_grpopts *opts, warninp_t wi);
-/* Generate adress groups */
-
char **read_rotparams(int *ninp_p, t_inpfile **inp, t_rot *rot, warninp_t wi);
/* Reads enforced rotation parameters, returns a list of the rot group names */
extern const char *eQMbasis_names[eQMbasisNR+1];
extern const char *eQMMMscheme_names[eQMMMschemeNR+1];
extern const char *eMultentOpt_names[eMultentOptNR+1];
-extern const char *eAdresstype_names[eAdressNR+1];
-extern const char *eAdressICtype_names[eAdressICNR+1];
-extern const char *eAdressSITEtype_names[eAdressSITENR+1];
extern const char *gmx_nblist_geometry_names[GMX_NBLIST_GEOMETRY_NR+1];
extern const char *gmx_nblist_interaction_names[GMX_NBLIST_INTERACTION_NR+1];
extern const char *gmx_nbkernel_elec_names[GMX_NBKERNEL_ELEC_NR+1];
#define EQMBASIS(e) ENUM_NAME(e, eQMbasisNR, eQMbasis_names)
#define EQMMMSCHEME(e) ENUM_NAME(e, eQMMMschemeNR, eQMMMscheme_names)
#define EMULTENTOPT(e) ENUM_NAME(e, eMultentOptNR, eMultentOpt_names)
-#define EADRESSTYPE(e) ENUM_NAME(e, eAdressNR, eAdresstype_names)
-#define EADRESSICTYPE(e) ENUM_NAME(e, eAdressICNR, eAdressICtype_names)
-#define EADRESSSITETYPE(e) ENUM_NAME(e, eAdressSITENR, eAdressSITEtype_names)
#define ELJPMECOMBNAMES(e) ENUM_NAME(e, eljpmeNR, eljpme_names)
#ifdef __cplusplus
efbposresCYLINDERX, efbposresCYLINDERY, efbposresCYLINDERZ, efbposresNR
};
-enum {
- eAdressOff, eAdressConst, eAdressXSplit, eAdressSphere, eAdressNR
-};
-
-enum {
- eAdressICOff, eAdressICThermoForce, eAdressICNR
-};
-
-enum {
- eAdressSITEcom, eAdressSITEcog, eAdressSITEatom, eAdressSITEatomatom, eAdressSITENR
-};
-
/* The interactions contained in a (possibly merged) table
* for computing electrostatic, VDW repulsion and/or VDW dispersion
{
GMX_NBLIST_INTERACTION_STANDARD,
GMX_NBLIST_INTERACTION_FREE_ENERGY,
- GMX_NBLIST_INTERACTION_ADRESS,
GMX_NBLIST_INTERACTION_NR
};
double t_wait;
int timesteps;
- /* parameter needed for AdResS simulation */
- int adress_type;
- gmx_bool badress_tf_full_box;
- real adress_const_wf;
- real adress_ex_width;
- real adress_hy_width;
- int adress_icor;
- int adress_site;
- rvec adress_refs;
- int n_adress_tf_grps;
- int * adress_tf_table_index;
- int *adress_group_explicit;
- t_forcetable * atf_tabs;
- real adress_ex_forcecap;
- gmx_bool adress_do_hybridpairs;
-
/* User determined parameters, copied from the inputrec */
int userint1;
int userint2;
#endif
-#define NO_TF_TABLE 255
-#define DEFAULT_TF_TABLE 0
-
typedef struct t_mdatoms {
real tmassA, tmassB, tmass;
int nr;
int homenr;
/* The lambda value used to create the contents of the struct */
real lambda;
- /* The AdResS weighting function */
- real *wf;
- unsigned short *tf_table_index; /* The tf table that will be applied, if thermodyn, force enabled*/
} t_mdatoms;
#ifdef __cplusplus
eNR_NBKERNEL_GENERIC = eNR_NBKERNEL_NR, /* Reuse number; KERNEL_NR is not an entry itself */
eNR_NBKERNEL_GENERIC_CG,
- eNR_NBKERNEL_GENERIC_ADRESS,
eNR_NBKERNEL_FREE_ENERGY, /* Add other generic kernels _before_ the free energy one */
eNR_NBKERNEL_ALLVSALL,
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2009 Christoph Junghans, Brad Lambeth.
- * Copyright (c) 2011,2012,2013,2014,2015, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-
-#include "gmxpre.h"
-
-#include "adress.h"
-
-#include <cmath>
-
-#include "gromacs/legacyheaders/types/forcerec.h"
-#include "gromacs/legacyheaders/types/ifunc.h"
-#include "gromacs/legacyheaders/types/mdatom.h"
-#include "gromacs/math/utilities.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/pbcutil/pbc.h"
-#include "gromacs/topology/topology.h"
-#include "gromacs/utility/fatalerror.h"
-
-real
-adress_weight(rvec x,
- int adresstype,
- real adressr,
- real adressw,
- rvec * ref,
- t_pbc * pbc,
- t_forcerec * fr )
-{
- int i;
- real l2 = adressr+adressw;
- real sqr_dl, dl;
- real tmp;
- rvec dx;
-
- sqr_dl = 0.0;
-
- if (pbc)
- {
- pbc_dx(pbc, (*ref), x, dx);
- }
- else
- {
- rvec_sub((*ref), x, dx);
- }
-
- switch (adresstype)
- {
- case eAdressOff:
- /* default to explicit simulation */
- return 1;
- case eAdressConst:
- /* constant value for weighting function = adressw */
- return fr->adress_const_wf;
- case eAdressXSplit:
- /* plane through center of ref, varies in x direction */
- sqr_dl = dx[0]*dx[0];
- break;
- case eAdressSphere:
- /* point at center of ref, assuming cubic geometry */
- for (i = 0; i < 3; i++)
- {
- sqr_dl += dx[i]*dx[i];
- }
- break;
- default:
- /* default to explicit simulation */
- return 1;
- }
-
- dl = std::sqrt(sqr_dl);
-
- /* molecule is coarse grained */
- if (dl > l2)
- {
- return 0;
- }
- /* molecule is explicit */
- else if (dl < adressr)
- {
- return 1;
- }
- /* hybrid region */
- else
- {
- tmp = std::cos((dl-adressr)*M_PI/2/adressw);
- return tmp*tmp;
- }
-}
-
-void
-update_adress_weights_com(FILE gmx_unused * fplog,
- int cg0,
- int cg1,
- t_block * cgs,
- rvec x[],
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- t_pbc * pbc)
-{
- int icg, k, k0, k1, d;
- real nrcg, inv_ncg, mtot, inv_mtot;
- int * cgindex;
- rvec ix;
- int adresstype;
- real adressr, adressw;
- rvec * ref;
- real * massT;
- real * wf;
-
-
- int n_hyb, n_ex, n_cg;
-
- n_hyb = 0;
- n_cg = 0;
- n_ex = 0;
-
- adresstype = fr->adress_type;
- adressr = fr->adress_ex_width;
- adressw = fr->adress_hy_width;
- massT = mdatoms->massT;
- wf = mdatoms->wf;
- ref = &(fr->adress_refs);
-
-
- /* Since this is center of mass AdResS, the vsite is not guaranteed
- * to be on the same node as the constructing atoms. Therefore we
- * loop over the charge groups, calculate their center of mass,
- * then use this to calculate wf for each atom. This wastes vsite
- * construction, but it's the only way to assure that the explicit
- * atoms have the same wf as their vsite. */
-
-#ifdef DEBUG
- fprintf(fplog, "Calculating center of mass for charge groups %d to %d\n",
- cg0, cg1);
-#endif
- cgindex = cgs->index;
-
- /* Compute the center of mass for all charge groups */
- for (icg = cg0; (icg < cg1); icg++)
- {
- k0 = cgindex[icg];
- k1 = cgindex[icg+1];
- nrcg = k1-k0;
- if (nrcg == 1)
- {
- wf[k0] = adress_weight(x[k0], adresstype, adressr, adressw, ref, pbc, fr);
- if (wf[k0] == 0)
- {
- n_cg++;
- }
- else if (wf[k0] == 1)
- {
- n_ex++;
- }
- else
- {
- n_hyb++;
- }
- }
- else
- {
- mtot = 0.0;
- for (k = k0; (k < k1); k++)
- {
- mtot += massT[k];
- }
- if (mtot > 0.0)
- {
- inv_mtot = 1.0/mtot;
-
- clear_rvec(ix);
- for (k = k0; (k < k1); k++)
- {
- for (d = 0; (d < DIM); d++)
- {
- ix[d] += x[k][d]*massT[k];
- }
- }
- for (d = 0; (d < DIM); d++)
- {
- ix[d] *= inv_mtot;
- }
- }
- /* Calculate the center of gravity if the charge group mtot=0 (only vsites) */
- else
- {
- inv_ncg = 1.0/nrcg;
-
- clear_rvec(ix);
- for (k = k0; (k < k1); k++)
- {
- for (d = 0; (d < DIM); d++)
- {
- ix[d] += x[k][d];
- }
- }
- for (d = 0; (d < DIM); d++)
- {
- ix[d] *= inv_ncg;
- }
- }
-
- /* Set wf of all atoms in charge group equal to wf of com */
- wf[k0] = adress_weight(ix, adresstype, adressr, adressw, ref, pbc, fr);
-
- if (wf[k0] == 0)
- {
- n_cg++;
- }
- else if (wf[k0] == 1)
- {
- n_ex++;
- }
- else
- {
- n_hyb++;
- }
-
- for (k = (k0+1); (k < k1); k++)
- {
- wf[k] = wf[k0];
- }
- }
- }
-}
-
-void update_adress_weights_atom_per_atom(
- int cg0,
- int cg1,
- t_block * cgs,
- rvec x[],
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- t_pbc * pbc)
-{
- int icg, k, k0, k1;
- int * cgindex;
- int adresstype;
- real adressr, adressw;
- rvec * ref;
- real * wf;
-
-
- int n_hyb, n_ex, n_cg;
-
- n_hyb = 0;
- n_cg = 0;
- n_ex = 0;
-
- adresstype = fr->adress_type;
- adressr = fr->adress_ex_width;
- adressw = fr->adress_hy_width;
- wf = mdatoms->wf;
- ref = &(fr->adress_refs);
-
- cgindex = cgs->index;
-
- /* Weighting function is determined for each atom individually.
- * This is an approximation
- * as in the theory requires an interpolation based on the center of masses.
- * Should be used with caution */
-
- for (icg = cg0; (icg < cg1); icg++)
- {
- k0 = cgindex[icg];
- k1 = cgindex[icg + 1];
-
- for (k = (k0); (k < k1); k++)
- {
- wf[k] = adress_weight(x[k], adresstype, adressr, adressw, ref, pbc, fr);
- if (wf[k] == 0)
- {
- n_cg++;
- }
- else if (wf[k] == 1)
- {
- n_ex++;
- }
- else
- {
- n_hyb++;
- }
- }
-
- }
-}
-
-void
-update_adress_weights_cog(t_iparams ip[],
- t_ilist ilist[],
- rvec x[],
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- t_pbc * pbc)
-{
- int i, j, nr, nra, inc;
- int ftype, adresstype;
- t_iatom avsite, ai, aj, ak, al;
- t_iatom * ia;
- real adressr, adressw;
- rvec * ref;
- real * wf;
- int n_hyb, n_ex, n_cg;
-
- adresstype = fr->adress_type;
- adressr = fr->adress_ex_width;
- adressw = fr->adress_hy_width;
- wf = mdatoms->wf;
- ref = &(fr->adress_refs);
-
-
- n_hyb = 0;
- n_cg = 0;
- n_ex = 0;
-
-
- /* Since this is center of geometry AdResS, we know the vsite
- * is in the same charge group node as the constructing atoms.
- * Loop over vsite types, calculate the weight of the vsite,
- * then assign that weight to the constructing atoms. */
-
- for (ftype = 0; (ftype < F_NRE); ftype++)
- {
- if (interaction_function[ftype].flags & IF_VSITE)
- {
- nra = interaction_function[ftype].nratoms;
- nr = ilist[ftype].nr;
- ia = ilist[ftype].iatoms;
-
- for (i = 0; (i < nr); )
- {
- /* The vsite and first constructing atom */
- avsite = ia[1];
- ai = ia[2];
- wf[avsite] = adress_weight(x[avsite], adresstype, adressr, adressw, ref, pbc, fr);
- wf[ai] = wf[avsite];
-
- if (wf[ai] == 0)
- {
- n_cg++;
- }
- else if (wf[ai] == 1)
- {
- n_ex++;
- }
- else
- {
- n_hyb++;
- }
-
- /* Assign the vsite wf to rest of constructing atoms depending on type */
- inc = nra+1;
- switch (ftype)
- {
- case F_VSITE2:
- aj = ia[3];
- wf[aj] = wf[avsite];
- break;
- case F_VSITE3:
- aj = ia[3];
- wf[aj] = wf[avsite];
- ak = ia[4];
- wf[ak] = wf[avsite];
- break;
- case F_VSITE3FD:
- aj = ia[3];
- wf[aj] = wf[avsite];
- ak = ia[4];
- wf[ak] = wf[avsite];
- break;
- case F_VSITE3FAD:
- aj = ia[3];
- wf[aj] = wf[avsite];
- ak = ia[4];
- wf[ak] = wf[avsite];
- break;
- case F_VSITE3OUT:
- aj = ia[3];
- wf[aj] = wf[avsite];
- ak = ia[4];
- wf[ak] = wf[avsite];
- break;
- case F_VSITE4FD:
- aj = ia[3];
- wf[aj] = wf[avsite];
- ak = ia[4];
- wf[ak] = wf[avsite];
- al = ia[5];
- wf[al] = wf[avsite];
- break;
- case F_VSITE4FDN:
- aj = ia[3];
- wf[aj] = wf[avsite];
- ak = ia[4];
- wf[ak] = wf[avsite];
- al = ia[5];
- wf[al] = wf[avsite];
- break;
- case F_VSITEN:
- inc = 3*ip[ia[0]].vsiten.n;
- for (j = 3; j < inc; j += 3)
- {
- ai = ia[j+2];
- wf[ai] = wf[avsite];
- }
- break;
- default:
- gmx_fatal(FARGS, "No such vsite type %d in %s, line %d",
- ftype, __FILE__, __LINE__);
- }
-
- /* Increment loop variables */
- i += inc;
- ia += inc;
- }
- }
- }
-}
-
-void
-update_adress_weights_atom(int cg0,
- int cg1,
- t_block * cgs,
- rvec x[],
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- t_pbc * pbc)
-{
- int icg, k, k0, k1;
- int * cgindex;
- int adresstype;
- real adressr, adressw;
- rvec * ref;
- real * wf;
-
- adresstype = fr->adress_type;
- adressr = fr->adress_ex_width;
- adressw = fr->adress_hy_width;
- wf = mdatoms->wf;
- ref = &(fr->adress_refs);
- cgindex = cgs->index;
-
- /* Only use first atom in charge group.
- * We still can't be sure that the vsite and constructing
- * atoms are on the same processor, so we must calculate
- * in the same way as com adress. */
-
- for (icg = cg0; (icg < cg1); icg++)
- {
- k0 = cgindex[icg];
- k1 = cgindex[icg+1];
- wf[k0] = adress_weight(x[k0], adresstype, adressr, adressw, ref, pbc, fr);
-
- /* Set wf of all atoms in charge group equal to wf of first atom in charge group*/
- for (k = (k0+1); (k < k1); k++)
- {
- wf[k] = wf[k0];
- }
- }
-}
-
-void
-adress_thermo_force(int start,
- int homenr,
- rvec x[],
- rvec f[],
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- t_pbc * pbc)
-{
- int iatom, n0, nnn, i;
- int adresstype;
- unsigned short * ptype;
- rvec * ref;
- real tabscale;
- real * ATFtab;
- rvec dr;
- real wt, rinv, sqr_dl, dl;
- real eps, eps2, F, Geps, Heps2, Fp, fscal;
-
- adresstype = fr->adress_type;
- ptype = mdatoms->ptype;
- ref = &(fr->adress_refs);
-
- for (iatom = start; (iatom < start+homenr); iatom++)
- {
- if (egp_coarsegrained(fr, mdatoms->cENER[iatom]))
- {
- if (ptype[iatom] == eptVSite)
- {
- /* is it hybrid or apply the thermodynamics force everywhere?*/
- if (mdatoms->tf_table_index[iatom] != NO_TF_TABLE)
- {
- if (fr->n_adress_tf_grps > 0)
- {
- /* multi component tf is on, select the right table */
- ATFtab = fr->atf_tabs[mdatoms->tf_table_index[iatom]].data;
- tabscale = fr->atf_tabs[mdatoms->tf_table_index[iatom]].scale;
- }
- else
- {
- /* just on component*/
- ATFtab = fr->atf_tabs[DEFAULT_TF_TABLE].data;
- tabscale = fr->atf_tabs[DEFAULT_TF_TABLE].scale;
- }
-
- if (pbc)
- {
- pbc_dx(pbc, (*ref), x[iatom], dr);
- }
- else
- {
- rvec_sub((*ref), x[iatom], dr);
- }
-
-
-
-
- /* calculate distace to adress center again */
- sqr_dl = 0.0;
- switch (adresstype)
- {
- case eAdressXSplit:
- /* plane through center of ref, varies in x direction */
- sqr_dl = dr[0]*dr[0];
- rinv = gmx_invsqrt(dr[0]*dr[0]);
- break;
- case eAdressSphere:
- /* point at center of ref, assuming cubic geometry */
- for (i = 0; i < 3; i++)
- {
- sqr_dl += dr[i]*dr[i];
- }
- rinv = gmx_invsqrt(sqr_dl);
- break;
- default:
- /* This case should not happen */
- rinv = 0.0;
- }
-
- dl = std::sqrt(sqr_dl);
- /* table origin is adress center */
- wt = dl*tabscale;
- n0 = static_cast<int>(wt);
- eps = wt-n0;
- eps2 = eps*eps;
- nnn = 4*n0;
- F = ATFtab[nnn+1];
- Geps = eps*ATFtab[nnn+2];
- Heps2 = eps2*ATFtab[nnn+3];
- Fp = F+Geps+Heps2;
- F = (Fp+Geps+2.0*Heps2)*tabscale;
-
- fscal = F*rinv;
-
- f[iatom][0] += fscal*dr[0];
- if (adresstype != eAdressXSplit)
- {
- f[iatom][1] += fscal*dr[1];
- f[iatom][2] += fscal*dr[2];
- }
- }
- }
- }
- }
-}
-
-gmx_bool egp_explicit(t_forcerec * fr, int egp_nr)
-{
- return fr->adress_group_explicit[egp_nr];
-}
-
-gmx_bool egp_coarsegrained(t_forcerec * fr, int egp_nr)
-{
- return !fr->adress_group_explicit[egp_nr];
-}
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2009 Christoph Junghans, Brad Lambeth.
- * Copyright (c) 2011,2014,2015, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-/*! \internal \file
- * \brief Implementation of the AdResS method.
- */
-#ifndef _adress_h_
-#define _adress_h_
-
-#include <cstdio>
-
-#include "gromacs/math/vectypes.h"
-#include "gromacs/topology/idef.h"
-#include "gromacs/utility/basedefinitions.h"
-#include "gromacs/utility/real.h"
-
-struct t_block;
-struct t_forcerec;
-struct t_mdatoms;
-struct t_pbc;
-
-/** \brief calculates the AdResS weight of a particle
- *
- * \param[in] x position of the particle
- * \param[in] adresstype type of address weight function
- * eAdressOff - explicit simulation
- * eAdressConst - constant weight all over the box
- * eAdressXSplit - split in x direction with ref as center
- * eAdressSphere - spherical splitting with ref as center
- * else - weight = 1 - explicit simulation
- * \param[in] adressr radius/size of the explicit zone
- * \param[in] adressw size of the hybrid zone
- * \param[in] ref center of the explicit zone
- * for adresstype 1 - unused
- * for adresstype 2 - only ref[0] is used
- * \param[in] pbc pbc struct for calculating shortest distance
- * \param[in] fr the forcerec containing all the parameters
- *
- * \return weight of the particle
- *
- */
-real
-adress_weight(rvec x,
- int adresstype,
- real adressr,
- real adressw,
- rvec * ref,
- struct t_pbc * pbc,
- t_forcerec * fr);
-
-/** \brief update the weight of all coarse-grained particles in several charge groups for com vsites
- *
- * \param[in,out] fplog log file in case of debug
- * \param[in] cg0 first charge group to update
- * \param[in] cg1 last+1 charge group to update
- * \param[in] cgs block containing the cg index
- * \param[in] x array with all the particle positions
- * \param[in] fr the forcerec containing all the parameters
- * \param[in,out] mdatoms the struct containing all the atoms properties
- * \param[in] pbc for shortest distance in adress_weight
- */
-void
-update_adress_weights_com(FILE * fplog,
- int cg0,
- int cg1,
- t_block * cgs,
- rvec x[],
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- struct t_pbc * pbc);
-
-/** \brief update the weight of all coarse-grained particles for cog vsites
- *
- * \param[in] ip contains interaction parameters, in this case the number of constructing atoms n for vsitesn
- * \param[in] ilist list of interaction types, in this case the virtual site types are what's important
- * \param[in] x array with all the particle positions
- * \param[in] fr the forcerec containing all the parameters
- * \param[in,out] mdatoms the struct containing all the atoms properties
- * \param[in] pbc for shortest distance in adress_weight
- */
-void
-update_adress_weights_cog(t_iparams ip[],
- t_ilist ilist[],
- rvec x[],
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- struct t_pbc * pbc);
-
-/** \brief update the weight of all coarse-grained particles in several charge groups for atom vsites
- *
- * \param[in] cg0 first charge group to update
- * \param[in] cg1 last+1 charge group to update
- * \param[in] cgs block containing the cg index
- * \param[in] x array with all the particle positions
- * \param[in] fr the forcerec containing all the parameters
- * \param[in,out] mdatoms the struct containing all the atoms properties
- * \param[in] pbc for shortest distance in adress_weight
- */
-void
-update_adress_weights_atom(int cg0,
- int cg1,
- t_block * cgs,
- rvec x[],
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- struct t_pbc * pbc);
-
-/** \brief update the weight on per atom basis of all coarse-grained particles in several charge groups for atom vsites
- *
- * \param[in] cg0 first charge group to update
- * \param[in] cg1 last+1 charge group to update
- * \param[in] cgs block containing the cg index
- * \param[in] x array with all the particle positions
- * \param[in] fr the forcerec containing all the parameters
- * \param[in,out] mdatoms the struct containing all the atoms properties
- * \param[in] pbc for shortest distance in adress_weight
- */
-void
-update_adress_weights_atom_per_atom(int cg0,
- int cg1,
- t_block * cgs,
- rvec x[],
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- struct t_pbc * pbc);
-
-/** \brief add AdResS IC thermodynamic force to f_novirsum
- *
- * \param[in] cg0 first charge group to update
- * \param[in] cg1 last+1 charge group to update
- * \param[in] x array with all the particle positions
- * \param[in,out] f the force array pointing at f_novirsum from sim_util.c
- * \param[in] fr the forcerec containing all the parameters
- * \param[in] mdatoms the struct containing all the atoms properties
- * \param[in] pbc for shortest distance to fr->adress_refs
- */
-void
-adress_thermo_force(int cg0,
- int cg1,
- rvec x[],
- rvec f[],
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- struct t_pbc * pbc);
-
-
-/** \brief checks weather a cpu calculates only coarse-grained or explicit interactions
- *
- * \param[in] n_ex number of explicit particles
- * \param[in] n_hyb number of hybrid particles
- * \param[in] n_cg number of coarse-grained particles
- * \param[in,out] mdatoms the struct containing all the atoms properties
- */
-void adress_set_kernel_flags(int n_ex, int n_hyb, int n_cg, t_mdatoms * mdatoms);
-
-/** \brief looks up if a energy group is explicit
- * \param[in] fr the forcerec containing all the parameters
- * \param[in] egp_nr energy group number
- * \return boolean if explicit or not
- */
-gmx_bool egp_explicit(t_forcerec * fr, int egp_nr);
-
-/** \brief looks up if a energy group is coarse-grained
- * \param[in] fr the forcerec containing all the parameters
- * \param[in] egp_nr energy group number
- * \return boolean if coarse-grained or not
- */
-gmx_bool egp_coarsegrained(t_forcerec * fr, int egp_nr);
-
-#endif
}
}
-static void bc_adress(const t_commrec *cr, t_adress *adress)
-{
- block_bc(cr, *adress);
- if (adress->n_tf_grps > 0)
- {
- snew_bc(cr, adress->tf_table_index, adress->n_tf_grps);
- nblock_bc(cr, adress->n_tf_grps, adress->tf_table_index);
- }
- if (adress->n_energy_grps > 0)
- {
- snew_bc(cr, adress->group_explicit, adress->n_energy_grps);
- nblock_bc(cr, adress->n_energy_grps, adress->group_explicit);
- }
-}
-
static void bc_imd(const t_commrec *cr, t_IMD *imd)
{
block_bc(cr, *imd);
snew_bc(cr, inputrec->swap, 1);
bc_swapions(cr, inputrec->swap);
}
- if (inputrec->bAdress)
- {
- snew_bc(cr, inputrec->adress, 1);
- bc_adress(cr, inputrec->adress);
- }
}
static void bc_moltype(const t_commrec *cr, t_symtab *symtab,
return cutoff;
}
-static void make_adress_tf_tables(FILE *fp,
- t_forcerec *fr, const t_inputrec *ir,
- const char *tabfn, const gmx_mtop_t *mtop,
- matrix box)
-{
- char buf[STRLEN];
- int i, j;
-
- if (tabfn == NULL)
- {
- gmx_fatal(FARGS, "No thermoforce table file given. Use -tabletf to specify a file\n");
- return;
- }
-
- snew(fr->atf_tabs, ir->adress->n_tf_grps);
-
- sprintf(buf, "%s", tabfn);
- for (i = 0; i < ir->adress->n_tf_grps; i++)
- {
- j = ir->adress->tf_table_index[i]; /* get energy group index */
- sprintf(buf + strlen(tabfn) - strlen(ftp2ext(efXVG)) - 1, "tf_%s.%s",
- *(mtop->groups.grpname[mtop->groups.grps[egcENER].nm_ind[j]]), ftp2ext(efXVG));
- if (fp)
- {
- fprintf(fp, "loading tf table for energygrp index %d from %s\n", ir->adress->tf_table_index[i], buf);
- }
- fr->atf_tabs[i] = make_atf_table(fp, fr, buf, box);
- }
-
-}
-
gmx_bool can_use_allvsall(const t_inputrec *ir, gmx_bool bPrintNote, t_commrec *cr, FILE *fp)
{
gmx_bool bAllvsAll;
const t_commrec *cr,
matrix box,
const char *tabfn,
- const char *tabafn,
const char *tabpfn,
const char *tabbfn,
const char *nbpu_opt,
fr->n_tpi = 0;
}
- /* Copy AdResS parameters */
if (ir->bAdress)
{
- fr->adress_type = ir->adress->type;
- fr->adress_const_wf = ir->adress->const_wf;
- fr->adress_ex_width = ir->adress->ex_width;
- fr->adress_hy_width = ir->adress->hy_width;
- fr->adress_icor = ir->adress->icor;
- fr->adress_site = ir->adress->site;
- fr->adress_ex_forcecap = ir->adress->ex_forcecap;
- fr->adress_do_hybridpairs = ir->adress->do_hybridpairs;
-
-
- snew(fr->adress_group_explicit, ir->adress->n_energy_grps);
- for (i = 0; i < ir->adress->n_energy_grps; i++)
- {
- fr->adress_group_explicit[i] = ir->adress->group_explicit[i];
- }
-
- fr->n_adress_tf_grps = ir->adress->n_tf_grps;
- snew(fr->adress_tf_table_index, fr->n_adress_tf_grps);
- for (i = 0; i < fr->n_adress_tf_grps; i++)
- {
- fr->adress_tf_table_index[i] = ir->adress->tf_table_index[i];
- }
- copy_rvec(ir->adress->refs, fr->adress_refs);
- }
- else
- {
- fr->adress_type = eAdressOff;
- fr->adress_do_hybridpairs = FALSE;
+ gmx_fatal(FARGS, "AdResS simulations are no longer supported");
}
/* Copy the user determined parameters */
fr->bF_NoVirSum = (EEL_FULL(fr->eeltype) || EVDW_PME(fr->vdwtype) ||
gmx_mtop_ftype_count(mtop, F_POSRES) > 0 ||
gmx_mtop_ftype_count(mtop, F_FBPOSRES) > 0 ||
- IR_ELEC_FIELD(*ir) ||
- (fr->adress_icor != eAdressICOff)
+ IR_ELEC_FIELD(*ir)
);
if (fr->cutoff_scheme == ecutsGROUP &&
}
}
- if (ir->adress)
- {
- fr->nnblists *= 2;
- }
-
snew(fr->nblists, fr->nnblists);
/* This code automatically gives table length tabext without cut-off's,
if (bSomeNormalNbListsAreInUse)
{
make_nbf_tables(fp, fr, rtab, tabfn, NULL, NULL, &fr->nblists[0]);
- if (ir->adress)
- {
- make_nbf_tables(fp, fr, rtab, tabfn, NULL, NULL, &fr->nblists[fr->nnblists/2]);
- }
if (!bMakeSeparate14Table)
{
fr->tab14 = fr->nblists[0].table_elec_vdw;
*mtop->groups.grpname[nm_ind[egi]],
*mtop->groups.grpname[nm_ind[egj]],
&fr->nblists[m]);
- if (ir->adress)
- {
- make_nbf_tables(fp, fr, rtab, tabfn,
- *mtop->groups.grpname[nm_ind[egi]],
- *mtop->groups.grpname[nm_ind[egj]],
- &fr->nblists[fr->nnblists/2+m]);
- }
m++;
}
else if (fr->nnblists > 1)
GMX_MAKETABLES_14ONLY);
}
- /* Read AdResS Thermo Force table if needed */
- if (fr->adress_icor == eAdressICThermoForce)
- {
- /* old todo replace */
-
- if (ir->adress->n_tf_grps > 0)
- {
- make_adress_tf_tables(fp, fr, ir, tabfn, mtop, box);
-
- }
- else
- {
- /* load the default table */
- snew(fr->atf_tabs, 1);
- fr->atf_tabs[DEFAULT_TF_TABLE] = make_atf_table(fp, fr, tabafn, box);
- }
- }
-
/* Wall stuff */
fr->nwall = ir->nwall;
if (ir->nwall && ir->wall_type == ewtTABLE)
if (cr->duty & DUTY_PP)
{
gmx_nonbonded_setup(fr, bGenericKernelOnly);
- /*
- if (ir->bAdress)
- {
- gmx_setup_adress_kernels(fp,bGenericKernelOnly);
- }
- */
}
/* Initialize the thread working data for bonded interactions */
* \param[in] mtop Molecular topology
* \param[in] cr Communication structures
* \param[in] box Simulation box
- * \param[in] tabfn Table potential file
- * \param[in] tabafn Table potential file for angles
- * \param[in] tabpfn Table potential file for proper dihedrals
- * \param[in] tabbfn Table potential file for bonds
+ * \param[in] tabfn Table potential file for non-bonded interactions
+ * \param[in] tabpfn Table potential file for pair interactions
+ * \param[in] tabbfn Table potential file for bonded interactions
* \param[in] nbpu_opt Nonbonded Processing Unit (GPU/CPU etc.)
* \param[in] bNoSolvOpt Do not use solvent optimization
* \param[in] print_force Print forces for atoms with force >= print_force
const t_commrec *cr,
matrix box,
const char *tabfn,
- const char *tabafn,
const char *tabpfn,
const char *tabbfn,
const char *nbpu_opt,
{
srenew(md->bQM, md->nalloc);
}
- if (ir->bAdress)
- {
- srenew(md->wf, md->nalloc);
- srenew(md->tf_table_index, md->nalloc);
- }
}
alook = gmx_mtop_atomlookup_init(mtop);
md->bQM[i] = FALSE;
}
}
- /* Initialize AdResS weighting functions to adressw */
- if (ir->bAdress)
- {
- md->wf[i] = 1.0;
- /* if no tf table groups specified, use default table */
- md->tf_table_index[i] = DEFAULT_TF_TABLE;
- if (ir->adress->n_tf_grps > 0)
- {
- /* if tf table groups specified, tf is only applied to thoose energy groups*/
- md->tf_table_index[i] = NO_TF_TABLE;
- /* check wether atom is in one of the relevant energy groups and assign a table index */
- for (g = 0; g < ir->adress->n_tf_grps; g++)
- {
- if (md->cENER[i] == ir->adress->tf_table_index[g])
- {
- md->tf_table_index[i] = g;
- }
- }
- }
- }
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
}
snew(md, 1);
- md->bVir = TRUE;
- md->bPress = TRUE;
- md->bSurft = TRUE;
- md->bMu = TRUE;
-
if (EI_DYNAMICS(ir->eI))
{
md->delta_t = ir->delta_t;
}
}
- /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
- if (ir->bAdress && !debug)
- {
- for (i = 0; i < F_NRE; i++)
- {
- md->bEner[i] = FALSE;
- if (i == F_EKIN)
- {
- md->bEner[i] = TRUE;
- }
- if (i == F_TEMP)
- {
- md->bEner[i] = TRUE;
- }
- }
- md->bVir = FALSE;
- md->bPress = FALSE;
- md->bSurft = FALSE;
- md->bMu = FALSE;
- }
-
md->f_nre = 0;
for (i = 0; i < F_NRE; i++)
{
md->isvir = get_ebin_space(md->ebin, asize(sv_nm), sv_nm, unit_energy);
md->ifvir = get_ebin_space(md->ebin, asize(fv_nm), fv_nm, unit_energy);
}
- if (md->bVir)
- {
- md->ivir = get_ebin_space(md->ebin, asize(vir_nm), vir_nm, unit_energy);
- }
- if (md->bPress)
- {
- md->ipres = get_ebin_space(md->ebin, asize(pres_nm), pres_nm, unit_pres_bar);
- }
- if (md->bSurft)
- {
- md->isurft = get_ebin_space(md->ebin, asize(surft_nm), surft_nm,
- unit_surft_bar);
- }
+ md->ivir = get_ebin_space(md->ebin, asize(vir_nm), vir_nm, unit_energy);
+ md->ipres = get_ebin_space(md->ebin, asize(pres_nm), pres_nm, unit_pres_bar);
+ md->isurft = get_ebin_space(md->ebin, asize(surft_nm), surft_nm,
+ unit_surft_bar);
if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
{
md->ipc = get_ebin_space(md->ebin, md->bTricl ? 6 : 3,
}
}
- n = groups->grps[egcENER].nr;
- /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
- if (!ir->bAdress)
- {
- /*standard simulation*/
- md->nEg = n;
- md->nE = (n*(n+1))/2;
- }
- else if (!debug)
- {
- /*AdResS simulation*/
- md->nU = 0;
- md->nEg = 0;
- md->nE = 0;
- md->nEc = 0;
- md->isvir = FALSE;
- }
+ n = groups->grps[egcENER].nr;
+ md->nEg = n;
+ md->nE = (n*(n+1))/2;
+
snew(md->igrp, md->nE);
if (md->nE > 1)
{
add_ebin(md->ebin, md->isvir, 9, svir[0], bSum);
add_ebin(md->ebin, md->ifvir, 9, fvir[0], bSum);
}
- if (md->bVir)
- {
- add_ebin(md->ebin, md->ivir, 9, vir[0], bSum);
- }
- if (md->bPress)
- {
- add_ebin(md->ebin, md->ipres, 9, pres[0], bSum);
- }
- if (md->bSurft)
- {
- tmp = (pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ];
- add_ebin(md->ebin, md->isurft, 1, &tmp, bSum);
- }
+ add_ebin(md->ebin, md->ivir, 9, vir[0], bSum);
+ add_ebin(md->ebin, md->ipres, 9, pres[0], bSum);
+ tmp = (pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ];
+ add_ebin(md->ebin, md->isurft, 1, &tmp, bSum);
if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
{
tmp6[0] = state->boxv[XX][XX];
pr_ebin(log, md->ebin, md->ifvir, 9, 3, mode, FALSE);
fprintf(log, "\n");
}
- if (md->bVir)
- {
- fprintf(log, " Total Virial (%s)\n", unit_energy);
- pr_ebin(log, md->ebin, md->ivir, 9, 3, mode, FALSE);
- fprintf(log, "\n");
- }
- if (md->bPress)
- {
- fprintf(log, " Pressure (%s)\n", unit_pres_bar);
- pr_ebin(log, md->ebin, md->ipres, 9, 3, mode, FALSE);
- fprintf(log, "\n");
- }
+ fprintf(log, " Total Virial (%s)\n", unit_energy);
+ pr_ebin(log, md->ebin, md->ivir, 9, 3, mode, FALSE);
+ fprintf(log, "\n");
+ fprintf(log, " Pressure (%s)\n", unit_pres_bar);
+ pr_ebin(log, md->ebin, md->ipres, 9, 3, mode, FALSE);
+ fprintf(log, "\n");
if (md->bMu)
{
fprintf(log, " Total Dipole (%s)\n", unit_dipole_D);
gmx_bool bMTTK;
gmx_bool bMu; /* true if dipole is calculated */
gmx_bool bDiagPres;
- gmx_bool bVir;
- gmx_bool bPress;
- gmx_bool bSurft;
int f_nre;
int epc;
real ref_p;
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
-#include "adress.h"
-
/*
* E X C L U S I O N H A N D L I N G
*/
{
nbl = &(fr->nblists[i]);
- if ((fr->adress_type != eAdressOff) && (i >= fr->nnblists/2))
- {
- type = GMX_NBLIST_INTERACTION_ADRESS;
- }
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],
}
}
-static void
-put_in_list_adress(gmx_bool bHaveVdW[],
- int ngid,
- t_mdatoms * md,
- int icg,
- int jgid,
- int nj,
- int jjcg[],
- int index[],
- t_excl bExcl[],
- int shift,
- t_forcerec * fr,
- gmx_bool bLR,
- gmx_bool bDoVdW,
- gmx_bool bDoCoul,
- int gmx_unused solvent_opt)
-{
- /* The a[] index has been removed,
- * to put it back in i_atom should be a[i0] and jj should be a[jj].
- */
- t_nblist * vdwc;
- t_nblist * vdw;
- t_nblist * coul;
- t_nblist * vdwc_adress = NULL;
- t_nblist * vdw_adress = NULL;
- t_nblist * coul_adress = NULL;
-
- int i, j, jcg, igid, gid, nbl_ind, nbl_ind_adress;
- int jj, jj0, jj1, i_atom;
- int i0, nicg;
-
- int *cginfo;
- int *type;
- real *charge;
- real *wf;
- real qi;
- gmx_bool bNotEx;
- gmx_bool bDoVdW_i, bDoCoul_i;
- gmx_bool b_hybrid;
- t_nblist *nlist, *nlist_adress;
- gmx_bool bEnergyGroupCG;
-
- /* Copy some pointers */
- cginfo = fr->cginfo;
- charge = md->chargeA;
- type = md->typeA;
- wf = md->wf;
-
- /* Get atom range */
- i0 = index[icg];
- nicg = index[icg+1]-i0;
-
- /* Get the i charge group info */
- igid = GET_CGINFO_GID(cginfo[icg]);
-
- if (md->nPerturbed)
- {
- gmx_fatal(FARGS, "AdResS does not support free energy pertubation\n");
- }
-
- /* Unpack pointers to neighbourlist structs */
- if (fr->nnblists == 2)
- {
- nbl_ind = 0;
- nbl_ind_adress = 1;
- }
- else
- {
- nbl_ind = fr->gid2nblists[GID(igid, jgid, ngid)];
- nbl_ind_adress = nbl_ind+fr->nnblists/2;
- }
- if (bLR)
- {
- nlist = fr->nblists[nbl_ind].nlist_lr;
- nlist_adress = fr->nblists[nbl_ind_adress].nlist_lr;
- }
- else
- {
- nlist = fr->nblists[nbl_ind].nlist_sr;
- nlist_adress = fr->nblists[nbl_ind_adress].nlist_sr;
- }
-
-
- vdwc = &nlist[eNL_VDWQQ];
- vdw = &nlist[eNL_VDW];
- coul = &nlist[eNL_QQ];
-
- vdwc_adress = &nlist_adress[eNL_VDWQQ];
- vdw_adress = &nlist_adress[eNL_VDW];
- coul_adress = &nlist_adress[eNL_QQ];
-
- /* We do not support solvent optimization with AdResS for now.
- For this we would need hybrid solvent-other kernels */
-
- /* no solvent as i charge group */
- /* Loop over the atoms in the i charge group */
- for (i = 0; i < nicg; i++)
- {
- i_atom = i0+i;
- gid = GID(igid, jgid, ngid);
- qi = charge[i_atom];
-
- /* Create new i_atom for each energy group */
- if (bDoVdW && bDoCoul)
- {
- new_i_nblist(vdwc, i_atom, shift, gid);
- new_i_nblist(vdwc_adress, i_atom, shift, gid);
-
- }
- if (bDoVdW)
- {
- new_i_nblist(vdw, i_atom, shift, gid);
- new_i_nblist(vdw_adress, i_atom, shift, gid);
-
- }
- if (bDoCoul)
- {
- new_i_nblist(coul, i_atom, shift, gid);
- new_i_nblist(coul_adress, i_atom, shift, gid);
- }
- bDoVdW_i = (bDoVdW && bHaveVdW[type[i_atom]]);
- bDoCoul_i = (bDoCoul && qi != 0);
-
- /* Here we find out whether the energy groups interaction belong to a
- * coarse-grained (vsite) or atomistic interaction. Note that, beacuse
- * interactions between coarse-grained and other (atomistic) energygroups
- * are excluded automatically by grompp, it is sufficient to check for
- * the group id of atom i (igid) */
- bEnergyGroupCG = !egp_explicit(fr, igid);
-
- if (bDoVdW_i || bDoCoul_i)
- {
- /* Loop over the j charge groups */
- for (j = 0; (j < nj); j++)
- {
- jcg = jjcg[j];
-
- /* Check for large charge groups */
- if (jcg == icg)
- {
- jj0 = i0 + i + 1;
- }
- else
- {
- jj0 = index[jcg];
- }
-
- jj1 = index[jcg+1];
- /* Finally loop over the atoms in the j-charge group */
- for (jj = jj0; jj < jj1; jj++)
- {
- bNotEx = NOTEXCL(bExcl, i, jj);
-
- /* Now we have to exclude interactions which will be zero
- * anyway due to the AdResS weights (in previous implementations
- * this was done in the force kernel). This is necessary as
- * pure interactions (those with b_hybrid=false, i.e. w_i*w_j==1 or 0)
- * are put into neighbour lists which will be passed to the
- * standard (optimized) kernels for speed. The interactions with
- * b_hybrid=true are placed into the _adress neighbour lists and
- * processed by the generic AdResS kernel.
- */
- if ( (bEnergyGroupCG &&
- wf[i_atom] >= 1-GMX_REAL_EPS && wf[jj] >= 1-GMX_REAL_EPS ) ||
- ( !bEnergyGroupCG && wf[jj] <= GMX_REAL_EPS ) )
- {
- continue;
- }
-
- b_hybrid = !((wf[i_atom] >= 1-GMX_REAL_EPS && wf[jj] >= 1-GMX_REAL_EPS) ||
- (wf[i_atom] <= GMX_REAL_EPS && wf[jj] <= GMX_REAL_EPS));
-
- if (bNotEx)
- {
- if (!bDoVdW_i)
- {
- if (charge[jj] != 0)
- {
- if (!b_hybrid)
- {
- add_j_to_nblist(coul, jj, bLR);
- }
- else
- {
- add_j_to_nblist(coul_adress, jj, bLR);
- }
- }
- }
- else if (!bDoCoul_i)
- {
- if (bHaveVdW[type[jj]])
- {
- if (!b_hybrid)
- {
- add_j_to_nblist(vdw, jj, bLR);
- }
- else
- {
- add_j_to_nblist(vdw_adress, jj, bLR);
- }
- }
- }
- else
- {
- if (bHaveVdW[type[jj]])
- {
- if (charge[jj] != 0)
- {
- if (!b_hybrid)
- {
- add_j_to_nblist(vdwc, jj, bLR);
- }
- else
- {
- add_j_to_nblist(vdwc_adress, jj, bLR);
- }
- }
- else
- {
- if (!b_hybrid)
- {
- add_j_to_nblist(vdw, jj, bLR);
- }
- else
- {
- add_j_to_nblist(vdw_adress, jj, bLR);
- }
-
- }
- }
- else if (charge[jj] != 0)
- {
- if (!b_hybrid)
- {
- add_j_to_nblist(coul, jj, bLR);
- }
- else
- {
- add_j_to_nblist(coul_adress, jj, bLR);
- }
-
- }
- }
- }
- }
- }
-
- close_i_nblist(vdw);
- close_i_nblist(coul);
- close_i_nblist(vdwc);
- close_i_nblist(vdw_adress);
- close_i_nblist(coul_adress);
- close_i_nblist(vdwc_adress);
- }
- }
-}
-
static void
put_in_list_qmmm(gmx_bool gmx_unused bHaveVdW[],
int ngid,
{
continue;
}
- /* Adress: an explicit cg that has a weigthing function of 0 is excluded
- * from the neigbour list as it will not interact */
- if (fr->adress_type != eAdressOff)
- {
- if (md->wf[cgs->index[icg]] <= GMX_REAL_EPS && egp_explicit(fr, igid))
- {
- continue;
- }
- }
/* Get shift vector */
shift = XYZ2IS(tx, ty, tz);
#ifdef NS5DB
fill_grid(NULL, ns->grid, fr->hcg, fr->hcg-1, fr->hcg, fr->cg_cm);
}
- if (fr->adress_type == eAdressOff)
+ if (!fr->ns.bCGlist)
{
- if (!fr->ns.bCGlist)
- {
- put_in_list = put_in_list_at;
- }
- else
- {
- put_in_list = put_in_list_cg;
- }
+ put_in_list = put_in_list_at;
}
else
{
- put_in_list = put_in_list_adress;
+ put_in_list = put_in_list_cg;
}
/* Do the core! */
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/sysinfo.h"
-#include "adress.h"
#include "nbnxn_gpu.h"
void print_time(FILE *out,
double mu[2*DIM];
gmx_bool bStateChanged, bNS, bFillGrid, bCalcCGCM;
gmx_bool bDoLongRangeNS, bDoForces, bSepLRF;
- gmx_bool bDoAdressWF;
- t_pbc pbc;
float cycles_pme, cycles_force;
start = 0;
bSepLRF = ((inputrec->nstcalclr > 1) && bDoForces &&
(flags & GMX_FORCE_SEPLRF) && (flags & GMX_FORCE_DO_LR));
- /* should probably move this to the forcerec since it doesn't change */
- bDoAdressWF = ((fr->adress_type != eAdressOff));
-
if (bStateChanged)
{
update_forcerec(fr, box);
wallcycle_stop(wcycle, ewcMOVEX);
}
- /* update adress weight beforehand */
- if (bStateChanged && bDoAdressWF)
- {
- /* need pbc for adress weight calculation with pbc_dx */
- set_pbc(&pbc, inputrec->ePBC, box);
- if (fr->adress_site == eAdressSITEcog)
- {
- update_adress_weights_cog(top->idef.iparams, top->idef.il, x, fr, mdatoms,
- inputrec->ePBC == epbcNONE ? NULL : &pbc);
- }
- else if (fr->adress_site == eAdressSITEcom)
- {
- update_adress_weights_com(fplog, cg0, cg1, &(top->cgs), x, fr, mdatoms,
- inputrec->ePBC == epbcNONE ? NULL : &pbc);
- }
- else if (fr->adress_site == eAdressSITEatomatom)
- {
- update_adress_weights_atom_per_atom(cg0, cg1, &(top->cgs), x, fr, mdatoms,
- inputrec->ePBC == epbcNONE ? NULL : &pbc);
- }
- else
- {
- update_adress_weights_atom(cg0, cg1, &(top->cgs), x, fr, mdatoms,
- inputrec->ePBC == epbcNONE ? NULL : &pbc);
- }
- }
-
if (NEED_MUTOT(*inputrec))
{
inputrec->ex, inputrec->et, t);
}
- if (bDoAdressWF && fr->adress_icor == eAdressICThermoForce)
- {
- /* Compute thermodynamic force in hybrid AdResS region */
- adress_thermo_force(start, homenr, x, fr->f_novirsum, fr, mdatoms,
- inputrec->ePBC == epbcNONE ? NULL : &pbc);
- }
-
/* Communicate the forces */
if (DOMAINDECOMP(cr))
{
mtop, ir, mdoutf_get_fp_dhdl(*outf));
}
- if (ir->bAdress)
- {
- please_cite(fplog, "Fritsch12");
- please_cite(fplog, "Junghans10");
- }
/* Initiate variables */
clear_mat(force_vir);
clear_mat(shake_vir);
int QMconstraints; /* constraints on QM bonds */
int QMMMscheme; /* Scheme: ONIOM or normal */
real scalefactor; /* factor for scaling the MM charges in QM calc.*/
- /* parameter needed for AdResS simulation */
- gmx_bool bAdress; /* Is AdResS enabled ? */
- t_adress *adress; /* The data for adress simulations */
+
+ /* Fields for removed features go here (better caching) */
+ gmx_bool bAdress;
} t_inputrec;
#define DEFORM(ir) ((ir).deform[XX][XX] != 0 || (ir).deform[YY][YY] != 0 || (ir).deform[ZZ][ZZ] != 0 || (ir).deform[YY][XX] != 0 || (ir).deform[ZZ][XX] != 0 || (ir).deform[ZZ][YY] != 0)
}
-t_forcetable make_atf_table(FILE *out,
- const t_forcerec *fr,
- const char *fn,
- matrix box)
-{
- t_tabledata *td;
- real rtab;
- real rx, ry, rz, box_r;
-
- t_forcetable table;
-
-
- /* Set the table dimensions for ATF, not really necessary to
- * use etiNR (since we only have one table, but ...)
- */
- snew(td, 1);
-
- if (fr->adress_type == eAdressSphere)
- {
- /* take half box diagonal direction as tab range */
- rx = 0.5*box[0][0]+0.5*box[1][0]+0.5*box[2][0];
- ry = 0.5*box[0][1]+0.5*box[1][1]+0.5*box[2][1];
- rz = 0.5*box[0][2]+0.5*box[1][2]+0.5*box[2][2];
- box_r = sqrt(rx*rx+ry*ry+rz*rz);
-
- }
- else
- {
- /* xsplit: take half box x direction as tab range */
- box_r = box[0][0]/2;
- }
- table.r = box_r;
- table.scale = 0;
- table.n = 0;
-
- read_tables(out, fn, 1, 0, td);
- rtab = td[0].x[td[0].nx-1];
-
- if (fr->adress_type == eAdressXSplit && (rtab < box[0][0]/2))
- {
- gmx_fatal(FARGS, "AdResS full box therm force table in file %s extends to %f:\n"
- "\tshould extend to at least half the length of the box in x-direction"
- "%f\n", fn, rtab, box[0][0]/2);
- }
- if (rtab < box_r)
- {
- gmx_fatal(FARGS, "AdResS full box therm force table in file %s extends to %f:\n"
- "\tshould extend to at least for spherical adress"
- "%f (=distance from center to furthermost point in box \n", fn, rtab, box_r);
- }
-
-
- table.n = td[0].nx;
- table.scale = td[0].tabscale;
-
- table.interaction = GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP;
- table.format = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH;
- table.formatsize = 4;
- table.ninteractions = 1;
- table.stride = table.formatsize*table.ninteractions;
-
- /* Each table type (e.g. coul,lj6,lj12) requires four numbers per
- * datapoint. For performance reasons we want the table data to be
- * aligned to a 32-byte boundary. This new pointer must not be
- * used in a free() call, but thankfully we're sloppy enough not
- * to do this :-)
- */
-
- snew_aligned(table.data, table.stride*table.n, 32);
-
- copy2table(table.n, 0, table.stride, td[0].x, td[0].v, td[0].f, 1.0, table.data);
-
- done_tabledata(&(td[0]));
- sfree(td);
-
- return table;
-}
-
bondedtable_t make_bonded_table(FILE *fplog, char *fn, int angle)
{
t_tabledata td;
/* Return a table for GB calculations */
t_forcetable make_gb_table(const t_forcerec *fr);
-/* Read a table for AdResS Thermo Force calculations */
-t_forcetable make_atf_table(FILE *out,
- const t_forcerec *fr,
- const char *fn,
- matrix box);
-
#endif /* GMX_TABLES_FORCETABLE_H */
}
}
}
-static void cmp_adress(FILE *fp, t_adress *ad1, t_adress *ad2,
- real ftol, real abstol)
-{
- cmp_int(fp, "ir->adress->type", -1, ad1->type, ad2->type);
- cmp_real(fp, "ir->adress->const_wf", -1, ad1->const_wf, ad2->const_wf, ftol, abstol);
- cmp_real(fp, "ir->adress->ex_width", -1, ad1->ex_width, ad2->ex_width, ftol, abstol);
- cmp_real(fp, "ir->adress->hy_width", -1, ad1->hy_width, ad2->hy_width, ftol, abstol);
- cmp_int(fp, "ir->adress->icor", -1, ad1->icor, ad2->icor);
- cmp_int(fp, "ir->adress->site", -1, ad1->site, ad2->site);
- cmp_rvec(fp, "ir->adress->refs", -1, ad1->refs, ad2->refs, ftol, abstol);
- cmp_real(fp, "ir->adress->ex_forcecap", -1, ad1->ex_forcecap, ad2->ex_forcecap, ftol, abstol);
-}
-
static void cmp_pull(FILE *fp)
{
fprintf(fp, "WARNING: Both files use COM pulling, but comparing of the pull struct is not implemented (yet). The pull parameters could be the same or different.\n");
cmp_rvec(fp, "inputrec->deform(c)", -1, ir1->deform[ZZ], ir2->deform[ZZ], ftol, abstol);
- cmp_bool(fp, "ir->bAdress->type", -1, ir1->bAdress, ir2->bAdress);
- if (ir1->bAdress && ir2->bAdress)
- {
- cmp_adress(fp, ir1->adress, ir2->adress, ftol, abstol);
- }
-
cmp_int(fp, "inputrec->userint1", -1, ir1->userint1, ir2->userint1);
cmp_int(fp, "inputrec->userint2", -1, ir1->userint2, ir2->userint2);
cmp_int(fp, "inputrec->userint3", -1, ir1->userint3, ir2->userint3);
{ efXVG, "-dhdl", "dhdl", ffOPTWR },
{ efXVG, "-field", "field", ffOPTWR },
{ efXVG, "-table", "table", ffOPTRD },
- { efXVG, "-tabletf", "tabletf", ffOPTRD },
{ efXVG, "-tablep", "tablep", ffOPTRD },
{ efXVG, "-tableb", "table", ffOPTRD },
{ efTRX, "-rerun", "rerun", ffOPTRD },
fr->gpu_opt = &hw_opt->gpu_opt;
init_forcerec(fplog, fr, fcd, inputrec, mtop, cr, box,
opt2fn("-table", nfile, fnm),
- opt2fn("-tabletf", nfile, fnm),
opt2fn("-tablep", nfile, fnm),
opt2fn("-tableb", nfile, fnm),
nbpu_opt,
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff