tmpi_get_source_list(THREAD_MPI_SOURCES ${CMAKE_SOURCE_DIR}/src/external/thread_mpi/src)
list(APPEND LIBGROMACS_SOURCES ${THREAD_MPI_SOURCES})
-list(APPEND LIBGROMACS_SOURCES ${TNG_SOURCES})
-tng_set_source_properties(WITH_ZLIB ${HAVE_ZLIB})
+if(GMX_USE_TNG)
+ list(APPEND LIBGROMACS_SOURCES ${TNG_SOURCES})
+ tng_set_source_properties(WITH_ZLIB ${HAVE_ZLIB})
+endif()
get_lmfit_properties(LMFIT_SOURCES LMFIT_LIBRARIES_TO_LINK LMFIT_INCLUDE_DIRECTORY LMFIT_INCLUDE_DIR_ORDER)
include_directories(${LMFIT_INCLUDE_DIR_ORDER} SYSTEM "${LMFIT_INCLUDE_DIRECTORY}")
if (dd->reverse_top->bExclRequired)
{
dd->nbonded_local += nexcl;
-
- forcerec_set_excl_load(fr, ltop);
}
ltop->atomtypes = mtop->atomtypes;
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016, 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.
#include "gromacs/mdtypes/forcerec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
/* There's nothing special to do here if just masses are perturbed,
* but if either charge or type is perturbed then the implementation
* perturbations. The parameter vectors for LJ-PME are likewise
* undefined when LJ-PME is not active. This works because
* bHaveChargeOrTypePerturbed handles the control flow. */
-void ewald_LRcorrection(int start, int end,
- t_commrec *cr, int thread, t_forcerec *fr,
+void ewald_LRcorrection(int numAtomsLocal,
+ t_commrec *cr,
+ int numThreads, int thread,
+ t_forcerec *fr,
real *chargeA, real *chargeB,
real *C6A, real *C6B,
real *sigmaA, real *sigmaB,
real lambda_q, real lambda_lj,
real *dvdlambda_q, real *dvdlambda_lj)
{
+ int numAtomsToBeCorrected;
+ if (calc_excl_corr)
+ {
+ /* We need to correct all exclusion pairs (cutoff-scheme = group) */
+ numAtomsToBeCorrected = excl->nr;
+
+ GMX_RELEASE_ASSERT(numAtomsToBeCorrected >= numAtomsLocal, "We might need to do self-corrections");
+ }
+ else
+ {
+ /* We need to correct only self interactions */
+ numAtomsToBeCorrected = numAtomsLocal;
+ }
+ int start = (numAtomsToBeCorrected* thread )/numThreads;
+ int end = (numAtomsToBeCorrected*(thread + 1))/numThreads;
+
int i, i1, i2, j, k, m, iv, jv, q;
int *AA;
double Vexcl_q, dvdl_excl_q, dvdl_excl_lj; /* Necessary for precision */
}
}
/* Dipole correction on force */
- if (dipole_coeff != 0)
+ if (dipole_coeff != 0 && i < numAtomsLocal)
{
for (j = 0; (j < DIM); j++)
{
}
}
/* Dipole correction on force */
- if (dipole_coeff != 0)
+ if (dipole_coeff != 0 && i < numAtomsLocal)
{
for (j = 0; (j < DIM); j++)
{
if (debug)
{
fprintf(debug, "Long Range corrections for Ewald interactions:\n");
- fprintf(debug, "start=%d,natoms=%d\n", start, end-start);
fprintf(debug, "q2sum = %g, Vself_q=%g c6sum = %g, Vself_lj=%g\n",
L1_q*fr->q2sum[0]+lambda_q*fr->q2sum[1], L1_q*Vself_q[0]+lambda_q*Vself_q[1], L1_lj*fr->c6sum[0]+lambda_lj*fr->c6sum[1], L1_lj*Vself_lj[0]+lambda_lj*Vself_lj[1]);
fprintf(debug, "Electrostatic Long Range correction: Vexcl=%g\n", Vexcl_q);
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016, 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.
* For both cutoff schemes, but only for Coulomb interactions,
* calculates correction for surface dipole terms. */
void
-ewald_LRcorrection(int start, int end,
- t_commrec *cr, int thread, t_forcerec *fr,
+ewald_LRcorrection(int numAtomsLocal,
+ t_commrec *cr,
+ int numThreads, int thread,
+ t_forcerec *fr,
real *chargeA, real *chargeB,
real *C6A, real *C6B,
real *sigmaA, real *sigmaB,
{
case efGRO:
out = gmx_fio_fopen(outfile, "w");
- write_hconf_indexed_p(out, title, atoms, nindex, index, 3, x, v, box);
+ write_hconf_indexed_p(out, title, atoms, nindex, index, x, v, box);
gmx_fio_fclose(out);
break;
case efG96:
switch (ftp)
{
case efGRO:
- write_conf_p(outfile, title, atoms, 3, x, v, box);
+ write_conf_p(outfile, title, atoms, x, v, box);
break;
case efG96:
clear_trxframe(&fr, TRUE);
{
case efGRO:
out = gmx_fio_fopen(outfile, "w");
- write_hconf_mtop(out, title, mtop, 3, x, v, box);
+ write_hconf_mtop(out, title, mtop, x, v, box);
gmx_fio_fclose(out);
break;
default:
gmx_fio_fclose(in);
}
+/* Note that the .gro reading routine still support variable precision
+ * for backward compatibility with old .gro files.
+ * We have removed writing of variable precision to avoid compatibility
+ * issues with other software packages.
+ */
static gmx_bool get_w_conf(FILE *in, const char *infile, char *title,
t_symtab *symtab, t_atoms *atoms, int *ndec,
rvec x[], rvec *v, matrix box)
return fr->natoms;
}
-static void make_hconf_format(int pr, gmx_bool bVel, char format[])
+static const char *get_hconf_format(bool haveVelocities)
{
- int l, vpr;
-
- /* build format string for printing,
- something like "%8.3f" for x and "%8.4f" for v */
- if (pr < 0)
- {
- pr = 0;
- }
- if (pr > 30)
- {
- pr = 30;
- }
- l = pr+5;
- vpr = pr+1;
- if (bVel)
+ if (haveVelocities)
{
- sprintf(format, "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
- l, pr, l, pr, l, pr, l, vpr, l, vpr, l, vpr);
+ return "%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n";
}
else
{
- sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
+ return "%8.3f%8.3f%8.3f\n";
}
}
-static void write_hconf_box(FILE *out, int pr, const matrix box)
+static void write_hconf_box(FILE *out, const matrix box)
{
- char format[100];
- int l;
-
- if (pr < 5)
- {
- pr = 5;
- }
- l = pr+5;
-
if (box[XX][YY] || box[XX][ZZ] || box[YY][XX] || box[YY][ZZ] ||
box[ZZ][XX] || box[ZZ][YY])
{
- sprintf(format, "%%%d.%df%%%d.%df%%%d.%df"
- "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
- l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr);
- fprintf(out, format,
+ fprintf(out, "%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n",
box[XX][XX], box[YY][YY], box[ZZ][ZZ],
box[XX][YY], box[XX][ZZ], box[YY][XX],
box[YY][ZZ], box[ZZ][XX], box[ZZ][YY]);
}
else
{
- sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
- fprintf(out, format,
+ fprintf(out, "%10.5f%10.5f%10.5f\n",
box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
}
}
void write_hconf_indexed_p(FILE *out, const char *title, const t_atoms *atoms,
- int nx, const int index[], int pr,
+ int nx, const int index[],
const rvec *x, const rvec *v, const matrix box)
{
- char resnm[6], nm[6], format[100];
+ char resnm[6], nm[6];
int ai, i, resind, resnr;
fprintf(out, "%s\n", (title && title[0]) ? title : gmx::bromacs().c_str());
fprintf(out, "%5d\n", nx);
- make_hconf_format(pr, v != NULL, format);
+ const char *format = get_hconf_format(v != NULL);
for (i = 0; (i < nx); i++)
{
}
}
- write_hconf_box(out, pr, box);
+ write_hconf_box(out, box);
fflush(out);
}
-void write_hconf_mtop(FILE *out, const char *title, gmx_mtop_t *mtop, int pr,
+void write_hconf_mtop(FILE *out, const char *title, gmx_mtop_t *mtop,
const rvec *x, const rvec *v, const matrix box)
{
- char format[100];
int i, resnr;
gmx_mtop_atomloop_all_t aloop;
t_atom *atom;
fprintf(out, "%s\n", (title && title[0]) ? title : gmx::bromacs().c_str());
fprintf(out, "%5d\n", mtop->natoms);
- make_hconf_format(pr, v != NULL, format);
+ const char *format = get_hconf_format(v != NULL);
aloop = gmx_mtop_atomloop_all_init(mtop);
while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
}
}
- write_hconf_box(out, pr, box);
+ write_hconf_box(out, box);
fflush(out);
}
-void write_hconf_p(FILE *out, const char *title, const t_atoms *atoms, int pr,
+void write_hconf_p(FILE *out, const char *title, const t_atoms *atoms,
const rvec *x, const rvec *v, const matrix box)
{
int *aa;
{
aa[i] = i;
}
- write_hconf_indexed_p(out, title, atoms, atoms->nr, aa, pr, x, v, box);
+ write_hconf_indexed_p(out, title, atoms, atoms->nr, aa, x, v, box);
sfree(aa);
}
void write_conf_p(const char *outfile, const char *title,
- const t_atoms *atoms, int pr,
+ const t_atoms *atoms,
const rvec *x, const rvec *v, const matrix box)
{
FILE *out;
out = gmx_fio_fopen(outfile, "w");
- write_hconf_p(out, title, atoms, pr, x, v, box);
+ write_hconf_p(out, title, atoms, x, v, box);
gmx_fio_fclose(out);
}
/* read first/next x and/or v frame from gro file */
void write_hconf_indexed_p(FILE *out, const char *title, const t_atoms *atoms,
- int nx, const int index[], int ndec,
+ int nx, const int index[],
const rvec *x, const rvec *v, const matrix box);
-void write_hconf_mtop(FILE *out, const char *title, struct gmx_mtop_t *mtop, int pr,
+void write_hconf_mtop(FILE *out, const char *title, struct gmx_mtop_t *mtop,
const rvec *x, const rvec *v, const matrix box);
-void write_hconf_p(FILE *out, const char *title, const t_atoms *atoms, int ndec,
+void write_hconf_p(FILE *out, const char *title, const t_atoms *atoms,
const rvec *x, const rvec *v, const matrix box);
/* Write a Gromos file with precision ndec: number of decimal places in x,
* v has one place more. */
void write_conf_p(const char *outfile, const char *title,
- const t_atoms *atoms, int pr,
+ const t_atoms *atoms,
const rvec *x, const rvec *v, const matrix box);
#ifdef __cplusplus
if (ftp == efGRO)
{
write_hconf_indexed_p(gmx_fio_getfp(status->fio), title, fr->atoms, nind, ind,
- prec2ndec(prec),
fr->x, fr->bV ? fr->v : NULL, fr->box);
}
else
if (gmx_fio_getftp(status->fio) == efGRO)
{
write_hconf_p(gmx_fio_getfp(status->fio), title, fr->atoms,
- prec2ndec(prec), fr->x, fr->bV ? fr->v : NULL, fr->box);
+ fr->x, fr->bV ? fr->v : NULL, fr->box);
}
else
{
}
}
+// Note it is the caller's responsibility to ensure that exiting is correct behaviour
+static gmx_noreturn void check_warning_error_impl(warninp_t wi, int f_errno, const char *file, int line)
+{
+ print_warn_count("note", wi->nwarn_note);
+ print_warn_count("warning", wi->nwarn_warn);
+
+ gmx_fatal(f_errno, file, line, "There %s %d error%s in input file(s)",
+ (wi->nwarn_error == 1) ? "was" : "were", wi->nwarn_error,
+ (wi->nwarn_error == 1) ? "" : "s");
+}
+
void check_warning_error(warninp_t wi, int f_errno, const char *file, int line)
{
if (wi->nwarn_error > 0)
{
- print_warn_count("note", wi->nwarn_note);
- print_warn_count("warning", wi->nwarn_warn);
-
- gmx_fatal(f_errno, file, line, "There %s %d error%s in input file(s)",
- (wi->nwarn_error == 1) ? "was" : "were", wi->nwarn_error,
- (wi->nwarn_error == 1) ? "" : "s");
+ check_warning_error_impl(wi, f_errno, file, line);
}
}
+void warning_error_and_exit(warninp_t wi, const char *s, int f_errno, const char *file, int line)
+{
+ warning_error(wi, s);
+ check_warning_error_impl(wi, f_errno, file, line);
+}
+
gmx_bool warning_errors_exist(warninp_t wi)
{
return (wi->nwarn_error > 0);
void done_warning(warninp_t wi, int f_errno, const char *file, int line)
{
+ // If we've had an error, then this will report the number of
+ // notes and warnings, and then exit.
+ check_warning_error(wi, f_errno, file, line);
+
+ // Otherwise, we report the number of notes and warnings.
print_warn_count("note", wi->nwarn_note);
print_warn_count("warning", wi->nwarn_warn);
- check_warning_error(wi, f_errno, file, line);
-
if (wi->maxwarn >= 0 && wi->nwarn_warn > wi->maxwarn)
{
gmx_fatal(f_errno, file, line,
* are printed, nwarn_error (local) is incremented.
*/
+/*! \brief Issue an error with warning_error() and prevent further
+ * processing by calling check_warning_error().
+ *
+ * This is intended for use where there is no way to produce a data
+ * structure that would prevent execution from segfaulting. */
+gmx_noreturn void warning_error_and_exit(warninp_t wi, const char *s, int f_errno, const char *file, int line);
+
gmx_bool warning_errors_exist(warninp_t wi);
/* Return whether any error-level warnings were issued to wi. */
fp = gmx_ffopen(outfile, "w");
if (!bOne)
{
- write_hconf_p(fp, *top1->name, atoms1, 3, x1, v1, box1);
+ write_hconf_p(fp, *top1->name, atoms1, x1, v1, box1);
}
- write_hconf_p(fp, *top2->name, atoms2, 3, x2, v2, box2);
+ write_hconf_p(fp, *top2->name, atoms2, x2, v2, box2);
gmx_ffclose(fp);
break;
default:
switch (ftp)
{
case efGRO:
- write_hconf_p(out, title, &useatoms, prec2ndec(frout.prec),
+ write_hconf_p(out, title, &useatoms,
frout.x, frout.bV ? frout.v : NULL, frout.box);
break;
case efPDB:
rpos[XX][mol], rpos[YY][mol], rpos[ZZ][mol]);
++mol;
++failed;
+ firstTrial = trial;
+ continue;
}
// Insert at positions taken from option -ip file.
offset_x[XX] = rpos[XX][mol] + deltaR[XX]*(2 * dist(rng)-1);
{
init_block2(&(block2[nmol-1]), mi0->atoms.nr);
}
- push_excl(pline, &(block2[nmol-1]));
+ push_excl(pline, &(block2[nmol-1]), wi);
break;
case d_system:
trim(pline);
mi0->plist,
&nnb,
&(mi0->excls));
- merge_excl(&(mi0->excls), &(block2[whichmol]));
+ merge_excl(&(mi0->excls), &(block2[whichmol]), wi);
done_block2(&(block2[whichmol]));
make_shake(mi0->plist, &mi0->atoms, opts->nshake);
convert_moltype_couple(mi0, dcatt, *fudgeQQ,
opts->couple_lam0, opts->couple_lam1,
opts->bCoupleIntra,
- nb_funct, &(plist[nb_funct]));
+ nb_funct, &(plist[nb_funct]), wi);
}
stupid_fill_block(&mi0->mols, mi0->atoms.nr, TRUE);
mi0->bProcessed = TRUE;
static void generate_qmexcl_moltype(gmx_moltype_t *molt, unsigned char *grpnr,
- t_inputrec *ir)
+ t_inputrec *ir, warninp_t wi)
{
/* This routine expects molt->ilist to be of size F_NRE and ordered. */
t_block2 qmexcl2;
init_block2(&qmexcl2, molt->atoms.nr);
b_to_b2(&qmexcl, &qmexcl2);
- merge_excl(&(molt->excls), &qmexcl2);
+ merge_excl(&(molt->excls), &qmexcl2, wi);
done_block2(&qmexcl2);
/* Finally, we also need to get rid of the pair interactions of the
/* Set the molecule type for the QMMM molblock */
molb->type = sys->nmoltype - 1;
}
- generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir);
+ generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir, wi);
}
if (grpnr)
{
int i, j, k = -1, nf;
int nr, nrfp;
real c, bi, bj, ci, cj, ci0, ci1, ci2, cj0, cj1, cj2;
- char errbuf[256];
+ char errbuf[STRLEN];
/* Lean mean shortcuts */
nr = get_atomtype_ntypes(atype);
break;
default:
- gmx_fatal(FARGS, "No such combination rule %d", comb);
+ sprintf(errbuf, "No such combination rule %d", comb);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
if (plist->nr != k)
{
double c[MAXFORCEPARAM];
double radius, vol, surftens, gb_radius, S_hct;
char tmpfield[12][100]; /* Max 12 fields of width 100 */
- char errbuf[256];
+ char errbuf[STRLEN];
t_atom *atom;
t_param *param;
int atomnr;
break;
default:
- gmx_fatal(FARGS, "Invalid function type %d in push_at %s %d", nb_funct,
- __FILE__, __LINE__);
+ sprintf(errbuf, "Invalid function type %d in push_at", nb_funct);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
for (j = nfp0; (j < MAXFORCEPARAM); j++)
{
if (strlen(type) == 1 && isdigit(type[0]))
{
- gmx_fatal(FARGS, "Atom type names can't be single digits.");
+ warning_error_and_exit(wi, "Atom type names can't be single digits.", FARGS);
}
if (strlen(btype) == 1 && isdigit(btype[0]))
{
- gmx_fatal(FARGS, "Bond atom type names can't be single digits.");
+ warning_error_and_exit(wi, "Bond atom type names can't be single digits.", FARGS);
}
/* Hack to read old topologies */
}
if (j == eptNR)
{
- gmx_fatal(FARGS, "Invalid particle type %s on line %s",
- ptype, line);
+ sprintf(errbuf, "Invalid particle type %s", ptype);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
pt = xl[j].ptype;
if (debug)
if ((nr = set_atomtype(nr, at, symtab, atom, type, param, batype_nr,
radius, vol, surftens, atomnr, gb_radius, S_hct)) == NOTSET)
{
- gmx_fatal(FARGS, "Replacing atomtype %s failed", type);
+ sprintf(errbuf, "Replacing atomtype %s failed", type);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
}
else if ((add_atomtype(at, symtab, atom, type, param,
batype_nr, radius, vol,
surftens, atomnr, gb_radius, S_hct)) == NOTSET)
{
- gmx_fatal(FARGS, "Adding atomtype %s failed", type);
+ sprintf(errbuf, "Adding atomtype %s failed", type);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
else
{
gmx_bool bTest, bFound, bCont, bId;
int nr = bt->nr;
int nrfp = NRFP(ftype);
- char errbuf[256];
+ char errbuf[STRLEN];
/* If bAllowRepeat is TRUE, we allow multiple entries as long as they
are on directly _adjacent_ lines.
/* One force parameter more, so we can check if we read too many */
double c[MAXFORCEPARAM+1];
t_param p;
- char errbuf[256];
+ char errbuf[STRLEN];
if ((bat && at) || (!bat && !at))
{
{
if (at && ((p.a[i] = get_atomtype_type(alc[i], at)) == NOTSET))
{
- gmx_fatal(FARGS, "Unknown atomtype %s\n", alc[i]);
+ sprintf(errbuf, "Unknown atomtype %s\n", alc[i]);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
{
- gmx_fatal(FARGS, "Unknown bond_atomtype %s\n", alc[i]);
+ sprintf(errbuf, "Unknown bond_atomtype %s\n", alc[i]);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
}
for (i = 0; (i < nrfp); i++)
double c[MAXFORCEPARAM];
t_param p;
gmx_bool bAllowRepeat;
- char errbuf[256];
+ char errbuf[STRLEN];
/* This routine accepts dihedraltypes defined from either 2 or 4 atoms.
*
{
if ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET)
{
- gmx_fatal(FARGS, "Unknown bond_atomtype %s", alc[i]);
+ sprintf(errbuf, "Unknown bond_atomtype %s", alc[i]);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
}
}
int ai, aj;
t_nbparam *nbp;
gmx_bool bId;
- char errbuf[256];
+ char errbuf[STRLEN];
if (sscanf (pline, "%s%s%d", a0, a1, &f) != 3)
{
}
else
{
- gmx_fatal(FARGS, "Number of force parameters for nonbonded interactions is %d"
- " in file %s, line %d", nrfp, __FILE__, __LINE__);
+ sprintf(errbuf, "Number of force parameters for nonbonded interactions is %d", nrfp);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
for (i = 0; (i < nrfp); i++)
{
/* Put the parameters in the matrix */
if ((ai = get_atomtype_type (a0, atype)) == NOTSET)
{
- gmx_fatal(FARGS, "Atomtype %s not found", a0);
+ sprintf(errbuf, "Atomtype %s not found", a0);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
if ((aj = get_atomtype_type (a1, atype)) == NOTSET)
{
- gmx_fatal(FARGS, "Atomtype %s not found", a1);
+ sprintf(errbuf, "Atomtype %s not found", a1);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
nbp = &(nbt[std::max(ai, aj)][std::min(ai, aj)]);
int nxcmap, nycmap, ncmap, read_cmap, sl, nct;
char s[20], alc[MAXATOMLIST+2][20];
t_param p;
- char errbuf[256];
+ char errbuf[STRLEN];
/* Keep the compiler happy */
read_cmap = 0;
/* Check for equal grid spacing in x and y dims */
if (nxcmap != nycmap)
{
- gmx_fatal(FARGS, "Not the same grid spacing in x and y for cmap grid: x=%d, y=%d", nxcmap, nycmap);
+ sprintf(errbuf, "Not the same grid spacing in x and y for cmap grid: x=%d, y=%d", nxcmap, nycmap);
+ warning_error(wi, errbuf);
}
ncmap = nxcmap*nycmap;
}
else
{
- gmx_fatal(FARGS, "Error in reading cmap parameter for angle %s %s %s %s %s", alc[0], alc[1], alc[2], alc[3], alc[4]);
+ sprintf(errbuf, "Error in reading cmap parameter for angle %s %s %s %s %s", alc[0], alc[1], alc[2], alc[3], alc[4]);
+ warning_error(wi, errbuf);
}
}
{
if (at && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
{
- gmx_fatal(FARGS, "Unknown atomtype %s\n", alc[i]);
+ sprintf(errbuf, "Unknown atomtype %s\n", alc[i]);
+ warning_error(wi, errbuf);
}
else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
{
- gmx_fatal(FARGS, "Unknown bond_atomtype %s\n", alc[i]);
+ sprintf(errbuf, "Unknown bond_atomtype %s\n", alc[i]);
+ warning_error(wi, errbuf);
}
/* Assign a grid number to each cmap_type */
/* Check for the correct number of atoms (again) */
if (bt[F_CMAP].nct != nct)
{
- gmx_fatal(FARGS, "Incorrect number of atom types (%d) in cmap type %d\n", nct, bt[F_CMAP].nc);
+ sprintf(errbuf, "Incorrect number of atom types (%d) in cmap type %d\n", nct, bt[F_CMAP].nc);
+ warning_error(wi, errbuf);
}
/* Is this correct?? */
int type, char *ctype, int ptype,
char *resnumberic,
char *resname, char *name, real m0, real q0,
- int typeB, char *ctypeB, real mB, real qB)
+ int typeB, char *ctypeB, real mB, real qB,
+ warninp_t wi)
{
int j, resind = 0, resnr;
unsigned char ric;
int nr = at->nr;
+ char errbuf[STRLEN];
if (((nr == 0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1)))
{
- gmx_fatal(FARGS, "Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)", atomnr, at->nr);
+ sprintf(errbuf, "Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)", atomnr, at->nr);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
j = strlen(resnumberic) - 1;
ric = resnumberic[j];
if (j == 0 || !isdigit(resnumberic[j-1]))
{
- gmx_fatal(FARGS, "Invalid residue number '%s' for atom %d",
- resnumberic, atomnr);
+ sprintf(errbuf, "Invalid residue number '%s' for atom %d",
+ resnumberic, atomnr);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
}
resnr = strtol(resnumberic, NULL, 10);
resnumberic[STRLEN], resname[STRLEN], name[STRLEN], check[STRLEN];
double m, q, mb, qb;
real m0, q0, mB, qB;
+ char errbuf[STRLEN];
/* Make a shortcut for writing in this molecule */
nr = at->nr;
sscanf(id, "%d", &atomnr);
if ((type = get_atomtype_type(ctype, atype)) == NOTSET)
{
- gmx_fatal(FARGS, "Atomtype %s not found", ctype);
+ sprintf(errbuf, "Atomtype %s not found", ctype);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
ptype = get_atomtype_ptype(type, atype);
{
if ((typeB = get_atomtype_type(ctypeB, atype)) == NOTSET)
{
- gmx_fatal(FARGS, "Atomtype %s not found", ctypeB);
+ sprintf(errbuf, "Atomtype %s not found", ctypeB);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
qB = get_atomtype_qA(typeB, atype);
mB = get_atomtype_massA(typeB, atype);
push_atom_now(symtab, at, atomnr, get_atomtype_atomnumber(type, atype),
type, ctype, ptype, resnumberic,
resname, name, m0, q0, typeB,
- typeB == type ? ctype : ctypeB, mB, qB);
+ typeB == type ? ctype : ctypeB, mB, qB, wi);
}
void push_molt(t_symtab *symtab, int *nmol, t_molinfo **mol, char *line,
char type[STRLEN];
int nrexcl, i;
t_molinfo *newmol;
+ char errbuf[STRLEN];
if ((sscanf(line, "%s%d", type, &nrexcl)) != 2)
{
{
if (strcmp(*((*mol)[i].name), type) == 0)
{
- gmx_fatal(FARGS, "moleculetype %s is redefined", type);
+ sprintf(errbuf, "moleculetype %s is redefined", type);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
i++;
}
static gmx_bool default_cmap_params(t_params bondtype[],
t_atoms *at, gpp_atomtype_t atype,
t_param *p, gmx_bool bB,
- int *cmap_type, int *nparam_def)
+ int *cmap_type, int *nparam_def,
+ warninp_t wi)
{
- int i, nparam_found;
- int ct;
- gmx_bool bFound = FALSE;
+ int i, nparam_found;
+ int ct;
+ gmx_bool bFound = FALSE;
+ char errbuf[STRLEN];
nparam_found = 0;
ct = 0;
/* If we did not find a matching type for this cmap torsion */
if (!bFound)
{
- gmx_fatal(FARGS, "Unknown cmap torsion between atoms %d %d %d %d %d\n",
- p->a[0]+1, p->a[1]+1, p->a[2]+1, p->a[3]+1, p->a[4]+1);
+ sprintf(errbuf, "Unknown cmap torsion between atoms %d %d %d %d %d",
+ p->a[0]+1, p->a[1]+1, p->a[2]+1, p->a[3]+1, p->a[4]+1);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
*nparam_def = nparam_found;
t_param param, *param_defA, *param_defB;
gmx_bool bFoundA = FALSE, bFoundB = FALSE, bDef, bPert, bSwapParity = FALSE;
int nparam_defA, nparam_defB;
- char errbuf[256];
+ char errbuf[STRLEN];
nparam_defA = nparam_defB = 0;
case F_VSITE3OUT:
break;
default:
- gmx_fatal(FARGS, "Negative function types only allowed for %s and %s",
- interaction_function[F_VSITE3FAD].longname,
- interaction_function[F_VSITE3OUT].longname);
+ sprintf(errbuf, "Negative function types only allowed for %s and %s",
+ interaction_function[F_VSITE3FAD].longname,
+ interaction_function[F_VSITE3OUT].longname);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
}
}
{
if (aa[i] < 1 || aa[i] > at->nr)
{
- gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
- "Atom index (%d) in %s out of bounds (1-%d).\n"
- "This probably means that you have inserted topology section \"%s\"\n"
- "in a part belonging to a different molecule than you intended to.\n"
- "In that case move the \"%s\" section to the right molecule.",
- get_warning_file(wi), get_warning_line(wi),
- aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
+ sprintf(errbuf, "Atom index (%d) in %s out of bounds (1-%d).\n"
+ "This probably means that you have inserted topology section \"%s\"\n"
+ "in a part belonging to a different molecule than you intended to.\n"
+ "In that case move the \"%s\" section to the right molecule.",
+ aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
+ warning_error_and_exit(wi, errbuf, FARGS);
}
for (j = i+1; (j < nral); j++)
{
if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
!(ftype == F_LJC14_Q && nread == 1))
{
- gmx_fatal(FARGS, "Incorrect number of parameters - found %d, expected %d or %d for %s.",
- nread, NRFPA(ftype), NRFP(ftype),
- interaction_function[ftype].longname);
+ sprintf(errbuf, "Incorrect number of parameters - found %d, expected %d or %d for %s.",
+ nread, NRFPA(ftype), NRFP(ftype), interaction_function[ftype].longname);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
for (j = 0; (j < nread); j++)
if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ)
&& param.c[5] != param.c[2])
{
- gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
- " %s multiplicity can not be perturbed %f!=%f",
- get_warning_file(wi), get_warning_line(wi),
- interaction_function[ftype].longname,
- param.c[2], param.c[5]);
+ sprintf(errbuf, "%s multiplicity can not be perturbed %f!=%f",
+ interaction_function[ftype].longname,
+ param.c[2], param.c[5]);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
if (IS_TABULATED(ftype) && param.c[2] != param.c[0])
{
- gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
- " %s table number can not be perturbed %d!=%d",
- get_warning_file(wi), get_warning_line(wi),
- interaction_function[ftype].longname,
- (int)(param.c[0]+0.5), (int)(param.c[2]+0.5));
+ sprintf(errbuf, "%s table number can not be perturbed %d!=%d",
+ interaction_function[ftype].longname,
+ (int)(param.c[0]+0.5), (int)(param.c[2]+0.5));
+ warning_error_and_exit(wi, errbuf, FARGS);
}
/* Dont add R-B dihedrals where all parameters are zero (no interaction) */
int i, j, ftype, nral, nread, ncmap_params;
int cmap_type;
int aa[MAXATOMLIST+1];
- char errbuf[256];
+ char errbuf[STRLEN];
gmx_bool bFound;
t_param param;
{
if (aa[i] < 1 || aa[i] > at->nr)
{
- gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
- "Atom index (%d) in %s out of bounds (1-%d).\n"
- "This probably means that you have inserted topology section \"%s\"\n"
- "in a part belonging to a different molecule than you intended to.\n"
- "In that case move the \"%s\" section to the right molecule.",
- get_warning_file(wi), get_warning_line(wi),
- aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
+ sprintf(errbuf, "Atom index (%d) in %s out of bounds (1-%d).\n"
+ "This probably means that you have inserted topology section \"%s\"\n"
+ "in a part belonging to a different molecule than you intended to.\n"
+ "In that case move the \"%s\" section to the right molecule.",
+ aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
+ warning_error_and_exit(wi, errbuf, FARGS);
}
for (j = i+1; (j < nral); j++)
}
/* Get the cmap type for this cmap angle */
- bFound = default_cmap_params(bondtype, at, atype, ¶m, FALSE, &cmap_type, &ncmap_params);
+ bFound = default_cmap_params(bondtype, at, atype, ¶m, FALSE, &cmap_type, &ncmap_params, wi);
/* We want exactly one parameter (the cmap type in state A (currently no state B) back */
if (bFound && ncmap_params == 1)
else
{
/* This is essentially the same check as in default_cmap_params() done one more time */
- gmx_fatal(FARGS, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
- param.a[0]+1, param.a[1]+1, param.a[2]+1, param.a[3]+1, param.a[4]+1);
+ sprintf(errbuf, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
+ param.a[0]+1, param.a[1]+1, param.a[2]+1, param.a[3]+1, param.a[4]+1);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
}
int *atc = NULL;
double *weight = NULL, weight_tot;
t_param param;
+ char errbuf[STRLEN];
/* default force parameters */
for (j = 0; (j < MAXATOMLIST); j++)
ptr += n;
if (ret == 0)
{
- gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
- " Expected an atom index in section \"%s\"",
- get_warning_file(wi), get_warning_line(wi),
- dir2str(d));
+ sprintf(errbuf, "Expected an atom index in section \"%s\"",
+ dir2str(d));
+ warning_error_and_exit(wi, errbuf, FARGS);
}
param.a[0] = a - 1;
ptr += n;
if (weight[nj] < 0)
{
- gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
- " No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
- get_warning_file(wi), get_warning_line(wi),
- nj+1, atc[nj]+1);
+ sprintf(errbuf, "No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
+ nj+1, atc[nj]+1);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
break;
default:
- gmx_fatal(FARGS, "Unknown vsiten type %d", type);
+ sprintf(errbuf, "Unknown vsiten type %d", type);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
weight_tot += weight[nj];
nj++;
if (nj == 0)
{
- gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
- " Expected more than one atom index in section \"%s\"",
- get_warning_file(wi), get_warning_line(wi),
- dir2str(d));
+ sprintf(errbuf, "Expected more than one atom index in section \"%s\"",
+ dir2str(d));
+ warning_error_and_exit(wi, errbuf, FARGS);
}
if (weight_tot == 0)
{
- gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
- " The total mass of the construting atoms is zero",
- get_warning_file(wi), get_warning_line(wi));
+ warning_error_and_exit(wi, "The total mass of the construting atoms is zero", FARGS);
}
for (j = 0; j < nj; j++)
warninp_t wi)
{
char type[STRLEN];
+ char errbuf[STRLEN];
if (sscanf(pline, "%s%d", type, nrcopies) != 2)
{
// avoid matching case-insensitive when we have multiple matches
if (nrci > 1)
{
- gmx_fatal(FARGS, "For moleculetype '%s' in [ system ] %d case insensitive matches, but %d case sensitive matches were found. Check the case of the characters in the moleculetypes.", type, nrci, nrcs);
+ sprintf(errbuf, "For moleculetype '%s' in [ system ] %d case insensitive matches, but %d case sensitive matches were found. Check the case of the characters in the moleculetypes.", type, nrci, nrcs);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
if (nrci == 1)
{
}
else
{
- gmx_fatal(FARGS, "No such moleculetype %s", type);
+ sprintf(errbuf, "No such moleculetype %s", type);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
}
}
}
}
-void push_excl(char *line, t_block2 *b2)
+void push_excl(char *line, t_block2 *b2, warninp_t wi)
{
int i, j;
int n;
char base[STRLEN], format[STRLEN];
+ char errbuf[STRLEN];
if (sscanf(line, "%d", &i) == 0)
{
}
else
{
- gmx_fatal(FARGS, "Invalid Atomnr j: %d, b2->nr: %d\n", j, b2->nr);
+ sprintf(errbuf, "Invalid Atomnr j: %d, b2->nr: %d\n", j, b2->nr);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
}
}
return (*((int *) v1))-(*((int *) v2));
}
-void merge_excl(t_blocka *excl, t_block2 *b2)
+void merge_excl(t_blocka *excl, t_block2 *b2, warninp_t wi)
{
int i, k;
int j;
int nra;
+ char errbuf[STRLEN];
if (!b2->nr)
{
}
else if (b2->nr != excl->nr)
{
- gmx_fatal(FARGS, "DEATH HORROR: b2->nr = %d, while excl->nr = %d",
- b2->nr, excl->nr);
+ sprintf(errbuf, "DEATH HORROR: b2->nr = %d, while excl->nr = %d",
+ b2->nr, excl->nr);
+ warning_error_and_exit(wi, errbuf, FARGS);
}
else if (debug)
{
plist[F_LJ14].param = NULL;
}
-static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp)
+static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp, warninp_t wi)
{
int n, ntype, i, j, k;
t_atom *atom;
t_blocka *excl;
gmx_bool bExcl;
t_param param;
+ char errbuf[STRLEN];
n = mol->atoms.nr;
atom = mol->atoms.atom;
{
if (nb_funct != F_LJ)
{
- gmx_fatal(FARGS, "Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
+ sprintf(errbuf, "Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
+ warning_error_and_exit(wi, errbuf, FARGS);
}
param.a[0] = i;
param.a[1] = j;
static void decouple_atoms(t_atoms *atoms, int atomtype_decouple,
int couple_lam0, int couple_lam1,
- const char *mol_name)
+ const char *mol_name, warninp_t wi)
{
- int i;
+ int i;
+ char errbuf[STRLEN];
for (i = 0; i < atoms->nr; i++)
{
if (atom->qB != atom->q || atom->typeB != atom->type)
{
- gmx_fatal(FARGS, "Atom %d in molecule type '%s' has different A and B state charges and/or atom types set in the topology file as well as through the mdp option '%s'. You can not use both these methods simultaneously.",
- i + 1, mol_name, "couple-moltype");
+ sprintf(errbuf, "Atom %d in molecule type '%s' has different A and B state charges and/or atom types set in the topology file as well as through the mdp option '%s'. You can not use both these methods simultaneously.",
+ i + 1, mol_name, "couple-moltype");
+ warning_error_and_exit(wi, errbuf, FARGS);
}
if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW)
void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, real fudgeQQ,
int couple_lam0, int couple_lam1,
- gmx_bool bCoupleIntra, int nb_funct, t_params *nbp)
+ gmx_bool bCoupleIntra, int nb_funct, t_params *nbp,
+ warninp_t wi)
{
convert_pairs_to_pairsQ(mol->plist, fudgeQQ, &mol->atoms);
if (!bCoupleIntra)
{
- generate_LJCpairsNB(mol, nb_funct, nbp);
+ generate_LJCpairsNB(mol, nb_funct, nbp, wi);
set_excl_all(&mol->excls);
}
decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1,
- *mol->name);
+ *mol->name, wi);
}
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016, 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.
void done_block2(t_block2 *b2);
-void push_excl(char *line, t_block2 *b2);
+void push_excl(char *line, t_block2 *b2, warninp_t wi);
-void merge_excl(t_blocka *excl, t_block2 *b2);
+void merge_excl(t_blocka *excl, t_block2 *b2, warninp_t wi);
void b_to_b2(t_blocka *b, t_block2 *b2);
real fudgeQQ,
int couple_lam0, int couple_lam1,
gmx_bool bCoupleIntra,
- int nb_funct, t_params *nbp);
+ int nb_funct, t_params *nbp,
+ warninp_t wi);
/* Setup mol such that the B-state has no interaction with the rest
* of the system, but full interaction with itself.
*/
/* This file is completely threadsafe - keep it that way! */
#include "gmxpre.h"
+#include "config.h" /* for GMX_MAX_OPENMP_THREADS */
+
#include <algorithm>
+#include "gromacs/math/vec.h"
#include "gromacs/math/vectypes.h"
#include "gromacs/mdlib/force.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
#include "gromacs/pbcutil/ishift.h"
#include "gromacs/pbcutil/mshift.h"
#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/utility/gmxassert.h"
#define XXXX 0
#define XXYY 1
vir[ZZ] -= 0.5*dvz;
}
-void calc_vir(int nxf, rvec x[], rvec f[], tensor vir,
- gmx_bool bScrewPBC, matrix box)
+static void calc_x_times_f(int nxf, const rvec x[], const rvec f[],
+ gmx_bool bScrewPBC, const matrix box,
+ matrix x_times_f)
{
- int i;
- double dvxx = 0, dvxy = 0, dvxz = 0, dvyx = 0, dvyy = 0, dvyz = 0, dvzx = 0, dvzy = 0, dvzz = 0;
+ clear_mat(x_times_f);
-#pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntDefault)) \
- schedule(static) \
- reduction(+: dvxx, dvxy, dvxz, dvyx, dvyy, dvyz, dvzx, dvzy, dvzz)
- for (i = 0; i < nxf; i++)
+ for (int i = 0; i < nxf; i++)
{
- dvxx += x[i][XX]*f[i][XX];
- dvxy += x[i][XX]*f[i][YY];
- dvxz += x[i][XX]*f[i][ZZ];
-
- dvyx += x[i][YY]*f[i][XX];
- dvyy += x[i][YY]*f[i][YY];
- dvyz += x[i][YY]*f[i][ZZ];
-
- dvzx += x[i][ZZ]*f[i][XX];
- dvzy += x[i][ZZ]*f[i][YY];
- dvzz += x[i][ZZ]*f[i][ZZ];
+ for (int d = 0; d < DIM; d++)
+ {
+ for (int n = 0; n < DIM; n++)
+ {
+ x_times_f[d][n] += x[i][d]*f[i][n];
+ }
+ }
if (bScrewPBC)
{
/* We should correct all odd x-shifts, but the range of isx is -2 to 2 */
if (isx == 1 || isx == -1)
{
- dvyy += box[YY][YY]*f[i][YY];
- dvyz += box[YY][YY]*f[i][ZZ];
-
- dvzy += box[ZZ][ZZ]*f[i][YY];
- dvzz += box[ZZ][ZZ]*f[i][ZZ];
+ for (int d = 0; d < DIM; d++)
+ {
+ for (int n = 0; n < DIM; n++)
+ {
+ x_times_f[d][n] += box[d][d]*f[i][n];
+ }
+ }
}
}
}
+}
- upd_vir(vir[XX], dvxx, dvxy, dvxz);
- upd_vir(vir[YY], dvyx, dvyy, dvyz);
- upd_vir(vir[ZZ], dvzx, dvzy, dvzz);
+void calc_vir(int nxf, rvec x[], rvec f[], tensor vir,
+ gmx_bool bScrewPBC, matrix box)
+{
+ matrix x_times_f;
+
+ int nthreads = gmx_omp_nthreads_get_simple_rvec_task(emntDefault, nxf*9);
+
+ GMX_ASSERT(nthreads >= 1, "Avoids uninitialized x_times_f (warning)");
+
+ if (nthreads == 1)
+ {
+ calc_x_times_f(nxf, x, f, bScrewPBC, box, x_times_f);
+ }
+ else
+ {
+ /* Use a buffer on the stack for storing thread-local results.
+ * We use 2 extra elements (=18 reals) per thread to separate thread
+ * local data by at least a cache line. Element 0 is not used.
+ */
+ matrix xf_buf[GMX_OPENMP_MAX_THREADS*3];
+
+#pragma omp parallel for num_threads(nthreads) schedule(static)
+ for (int thread = 0; thread < nthreads; thread++)
+ {
+ int start = (nxf*thread)/nthreads;
+ int end = std::min(nxf*(thread + 1)/nthreads, nxf);
+
+ calc_x_times_f(end - start, x + start, f + start, bScrewPBC, box,
+ // cppcheck-suppress uninitvar
+ thread == 0 ? x_times_f : xf_buf[thread*3]);
+ }
+
+ for (int thread = 1; thread < nthreads; thread++)
+ {
+ m_add(x_times_f, xf_buf[thread*3], x_times_f);
+ }
+ }
+
+ for (int d = 0; d < DIM; d++)
+ {
+ upd_vir(vir[d], x_times_f[d][XX], x_times_f[d][YY], x_times_f[d][ZZ]);
+ }
}
* exclusion forces) are calculated, so we can store
* the forces in the normal, single fr->f_novirsum array.
*/
- ewald_LRcorrection(fr->excl_load[t], fr->excl_load[t+1],
- cr, t, fr,
+ ewald_LRcorrection(md->homenr, cr, nthreads, t, fr,
md->chargeA, md->chargeB,
md->sqrt_c6A, md->sqrt_c6B,
md->sigmaA, md->sigmaB,
fr->nthread_ewc = gmx_omp_nthreads_get(emntBonded);
snew(fr->ewc_t, fr->nthread_ewc);
- snew(fr->excl_load, fr->nthread_ewc + 1);
/* fr->ic is used both by verlet and group kernels (to some extent) now */
init_interaction_const(fp, &fr->ic, fr);
fflush(fp);
}
-void forcerec_set_excl_load(t_forcerec *fr,
- const gmx_localtop_t *top)
-{
- const int *ind, *a;
- int t, i, j, ntot, n, ntarget;
-
- ind = top->excls.index;
- a = top->excls.a;
-
- ntot = 0;
- for (i = 0; i < top->excls.nr; i++)
- {
- for (j = ind[i]; j < ind[i+1]; j++)
- {
- if (a[j] > i)
- {
- ntot++;
- }
- }
- }
-
- fr->excl_load[0] = 0;
- n = 0;
- i = 0;
- for (t = 1; t <= fr->nthread_ewc; t++)
- {
- ntarget = (ntot*t)/fr->nthread_ewc;
- while (i < top->excls.nr && n < ntarget)
- {
- for (j = ind[i]; j < ind[i+1]; j++)
- {
- if (a[j] > i)
- {
- n++;
- }
- }
- i++;
- }
- fr->excl_load[t] = i;
- }
-}
-
/* Frees GPU memory and destroys the GPU context.
*
* Note that this function needs to be called even if GPUs are not used
*top = gmx_mtop_generate_local_top(top_global, ir->efep != efepNO);
- forcerec_set_excl_load(fr, *top);
-
setup_bonded_threading(fr, &(*top)->idef);
if (ir->ePBC != epbcNONE && !fr->bMolPBC)
#endif /* LJ_EWALD_COMB_GEOM */
#endif /* LJ_EWALD */
+#ifdef LJ_POT_SWITCH
+#ifdef CALC_ENERGIES
+ calculate_potential_switch_F_E(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
+#else
+ calculate_potential_switch_F(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
+#endif /* CALC_ENERGIES */
+#endif /* LJ_POT_SWITCH */
+
#ifdef VDW_CUTOFF_CHECK
/* Separate VDW cut-off check to enable twin-range cut-offs
* (rvdw < rcoulomb <= rlist)
#endif
#endif /* VDW_CUTOFF_CHECK */
-#ifdef LJ_POT_SWITCH
-#ifdef CALC_ENERGIES
- calculate_potential_switch_F_E(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
-#else
- calculate_potential_switch_F(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
-#endif /* CALC_ENERGIES */
-#endif /* LJ_POT_SWITCH */
-
#ifdef CALC_ENERGIES
E_lj += E_lj_p;
#endif
#endif /* LJ_EWALD_COMB_GEOM */
#endif /* LJ_EWALD */
+#ifdef LJ_POT_SWITCH
+#ifdef CALC_ENERGIES
+ calculate_potential_switch_F_E(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
+#else
+ calculate_potential_switch_F(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
+#endif /* CALC_ENERGIES */
+#endif /* LJ_POT_SWITCH */
+
#ifdef VDW_CUTOFF_CHECK
/* Separate VDW cut-off check to enable twin-range cut-offs
* (rvdw < rcoulomb <= rlist)
#endif
#endif /* VDW_CUTOFF_CHECK */
-#ifdef LJ_POT_SWITCH
-#ifdef CALC_ENERGIES
- calculate_potential_switch_F_E(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
-#else
- calculate_potential_switch_F(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
-#endif /* CALC_ENERGIES */
-#endif /* LJ_POT_SWITCH */
-
#ifdef CALC_ENERGIES
E_lj += E_lj_p;
#endif
#endif /* LJ_EWALD_COMB_GEOM */
#endif /* LJ_EWALD */
+#ifdef LJ_POT_SWITCH
+#ifdef CALC_ENERGIES
+ calculate_potential_switch_F_E(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
+#else
+ calculate_potential_switch_F(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
+#endif /* CALC_ENERGIES */
+#endif /* LJ_POT_SWITCH */
+
#ifdef VDW_CUTOFF_CHECK
/* Separate VDW cut-off check to enable twin-range cut-offs
* (rvdw < rcoulomb <= rlist)
#endif
#endif /* VDW_CUTOFF_CHECK */
-#ifdef LJ_POT_SWITCH
-#ifdef CALC_ENERGIES
- calculate_potential_switch_F_E(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
-#else
- calculate_potential_switch_F(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
-#endif /* CALC_ENERGIES */
-#endif /* LJ_POT_SWITCH */
-
#ifdef CALC_ENERGIES
E_lj += E_lj_p;
#endif /* LJ_EWALD_COMB_GEOM */
#endif /* LJ_EWALD */
+#ifdef LJ_POT_SWITCH
+#ifdef CALC_ENERGIES
+ calculate_potential_switch_F_E(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
+#else
+ calculate_potential_switch_F(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
+#endif /* CALC_ENERGIES */
+#endif /* LJ_POT_SWITCH */
+
#ifdef VDW_CUTOFF_CHECK
/* Separate VDW cut-off check to enable twin-range cut-offs
* (rvdw < rcoulomb <= rlist)
#endif
#endif /* VDW_CUTOFF_CHECK */
-#ifdef LJ_POT_SWITCH
-#ifdef CALC_ENERGIES
- calculate_potential_switch_F_E(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
-#else
- calculate_potential_switch_F(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
-#endif /* CALC_ENERGIES */
-#endif /* LJ_POT_SWITCH */
-
#ifdef CALC_ENERGIES
E_lj += E_lj_p;
#endif /* LJ_EWALD_COMB_GEOM */
#endif /* LJ_EWALD */
+#ifdef LJ_POT_SWITCH
+#ifdef CALC_ENERGIES
+ calculate_potential_switch_F_E(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
+#else
+ calculate_potential_switch_F(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
+#endif /* CALC_ENERGIES */
+#endif /* LJ_POT_SWITCH */
+
#ifdef VDW_CUTOFF_CHECK
/* Separate VDW cut-off check to enable twin-range cut-offs
* (rvdw < rcoulomb <= rlist)
#endif
#endif /* VDW_CUTOFF_CHECK */
-#ifdef LJ_POT_SWITCH
-#ifdef CALC_ENERGIES
- calculate_potential_switch_F_E(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
-#else
- calculate_potential_switch_F(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
-#endif /* CALC_ENERGIES */
-#endif /* LJ_POT_SWITCH */
-
#ifdef CALC_ENERGIES
E_lj += E_lj_p;
};
-static void do_update_md(int start, int nrend, double dt,
+static void do_update_md(int start, int nrend,
+ double dt, int nstpcouple,
t_grp_tcstat *tcstat,
double nh_vxi[],
gmx_bool bNEMD, t_grp_acc *gstat, rvec accel[],
if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
{
vnrel = (lg*vrel[d] + dt*(imass*f[n][d] - 0.5*vxi*vrel[d]
- - iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
+ - nstpcouple*iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
/* do not scale the mean velocities u */
vn = gstat[ga].u[d] + accel[ga][d]*dt + vnrel;
v[n][d] = vn;
}
} /* do_update_vv_pos */
-static void do_update_visc(int start, int nrend, double dt,
+static void do_update_visc(int start, int nrend,
+ double dt, int nstpcouple,
t_grp_tcstat *tcstat,
double nh_vxi[],
real invmass[],
if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
{
vn = (lg*vrel[d] + dt*(imass*f[n][d] - 0.5*vxi*vrel[d]
- - iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
+ - nstpcouple*iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
if (d == XX)
{
vn += vc + dt*cosz*cos_accel;
case (eiMD):
if (ekind->cosacc.cos_accel == 0)
{
- do_update_md(start_th, end_th, dt,
+ do_update_md(start_th, end_th,
+ dt, inputrec->nstpcouple,
ekind->tcstat, state->nosehoover_vxi,
ekind->bNEMD, ekind->grpstat, inputrec->opts.acc,
inputrec->opts.nFreeze,
}
else
{
- do_update_visc(start_th, end_th, dt,
+ do_update_visc(start_th, end_th,
+ dt, inputrec->nstpcouple,
ekind->tcstat, state->nosehoover_vxi,
md->invmass, md->ptype,
md->cTC, state->x, upd->xp, state->v, f, M,
/* If the energy group pair is excluded, we don't need a table */
if (!(fr->egp_flags[egp*ir->opts.ngener+negp_pp+w] & EGP_EXCL))
{
- fr->wall_tab[w][egp] = make_tables(fplog, fr, buf, 0,
- GMX_MAKETABLES_FORCEUSER);
sprintf(buf, "%s", tabfn);
sprintf(buf + strlen(tabfn) - strlen(ftp2ext(efXVG)) - 1, "_%s_%s.%s",
*groups->grpname[nm_ind[egp]],
*groups->grpname[nm_ind[negp_pp+w]],
ftp2ext(efXVG));
+ fr->wall_tab[w][egp] = make_tables(fplog, fr, buf, 0,
+ GMX_MAKETABLES_FORCEUSER);
/* Since wall have no charge, we can compress the table */
for (int i = 0; i <= fr->wall_tab[w][egp]->n; i++)
/* Ewald correction thread local virial and energy data */
int nthread_ewc;
ewald_corr_thread_t *ewc_t;
- /* Ewald charge correction load distribution over the threads */
- int *excl_load;
} t_forcerec;
/* Important: Starting with Gromacs-4.6, the values of c6 and c12 in the nbfp array have
{
top = gmx_mtop_generate_local_top(top_global, ir->efep != efepNO);
- forcerec_set_excl_load(fr, top);
-
state = serial_init_local_state(state_global);
atoms2md(top_global, ir, 0, NULL, top_global->natoms, mdatoms);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff