if (CMAKE_CXX_COMPILER_ID MATCHES "XL" AND CMAKE_CXX_COMPILER_LOADED)
GMX_TEST_CXXFLAG(XLC_BLUEGENEQ_CXXFLAG "-qarch=qp -qtune=qp" ACCELERATION_CXX_FLAGS)
try_compile(TEST_QPX ${CMAKE_BINARY_DIR}
- "cmake/TestQPX.c"
- COMPILE_DEFINITIONS "${ACCELERATION_CXX_FLAGS")
+ "${CMAKE_SOURCE_DIR}/cmake/TestQPX.c"
+ COMPILE_DEFINITIONS "${ACCELERATION_CXX_FLAGS}")
if(NOT TEST_QPX)
message(FATAL_ERROR "Cannot compile the requested IBM QPX intrinsics.")
endif()
set(COREWRAP_INCLUDE_DIR "${CMAKE_SOURCE_DIR}/../corewrap" CACHE STRING
"Path to swindirect.h")
include_directories(${COREWRAP_INCLUDE_DIR})
- set_property(CACHE GMX_COOL_QUOTES VALUE OFF)
+ set_property(CACHE GMX_COOL_QUOTES PROPERTY VALUE OFF)
endif(GMX_FAHCORE)
# # # # # # # # # # NO MORE TESTS AFTER THIS LINE! # # # # # # # # # # #
endif (CMAKE_USE_PTHREADS_INIT)
# Turns on thread_mpi core threading functions.
-# options are:
MACRO(TMPI_ENABLE_CORE INCDIR)
TMPI_TEST_ATOMICS(${INCDIR})
ENDMACRO(TMPI_ENABLE_CORE)
+++ /dev/null
-/*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * If you want to redistribute modifications, 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 www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
- *
- * And Hey:
- * Green Red Orange Magenta Azure Cyan Skyblue
- */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stddef.h>
-
-#include "types/commrec.h"
-#include "vec.h"
-#include "calcpot.h"
-#include "nrnb.h"
-#include "mdebin.h"
-#include "mshift.h"
-#include "smalloc.h"
-#include "force.h"
-#include "main.h"
-#include "filenm.h"
-#include "gmx_fatal.h"
-#include "mdrun.h"
-#include "ns.h"
-#include "txtdump.h"
-#include "mdatoms.h"
-#include "main.h"
-#include "mtop_util.h"
-#include "chargegroup.h"
-
-static void c_tabpot(real tabscale, real VFtab[],
- int nri, int iinr[],
- int shift[],
- int jindex[], int jjnr[],
- real pos[],
- real facel, real charge[],
- real pot[], real shiftvec[])
-{
- /* Local variables */
- const real nul = 0.000000;
-
- /* Table stuff */
- const real one = 1.000000;
- const real two = 2.000000;
- real r1, r1t, fijC, eps, eps2, Y, F, Fp, Geps, Heps2, VV, FF;
- int n0, n1, nnn;
-
- /* General and coulomb stuff */
- int ii, k, n, jnr, ii3, nj0, nj1, is3, j3, ggid;
- real fxJ, fyJ, fzJ, fxO, fyO, fzO;
- real ixO, iyO, izO, dxO, dyO, dzO;
- real txO, tyO, tzO, vcO, fsO, qO, rsqO, rinv1O, rinv2O;
- real qqO, qj;
- real jx, jy, jz, shX, shY, shZ, poti;
-
- /* Outer loop (over i particles) starts here */
- for (n = 0; (n < nri); n++)
- {
-
- /* Unpack shift vector */
- is3 = 3*shift[n];
- shX = shiftvec[is3];
- shY = shiftvec[is3+1];
- shZ = shiftvec[is3+2];
-
- /* Unpack I particle */
- ii = iinr[n];
- ii3 = 3*ii;
-
- /* Charge of i particle(s) divided by 4 pi eps0 */
- qO = facel*charge[ii];
-
- /* Bounds for the innerloop */
- nj0 = jindex[n];
- nj1 = jindex[n+1];
-
- /* Compute shifted i position */
- ixO = shX + pos[ii3];
- iyO = shY + pos[ii3+1];
- izO = shZ + pos[ii3+2];
- poti = nul;
-
- /* Inner loop (over j-particles) starts right here */
- for (k = nj0; (k < nj1); k++)
- {
-
- /* Unpack neighbourlist */
- jnr = jjnr[k];
- j3 = 3*jnr;
- qj = facel*charge[jnr];
- jx = pos[j3];
- jy = pos[j3+1];
- jz = pos[j3+2];
-
- /* First one is for oxygen, with LJ */
- dxO = ixO - jx;
- dyO = iyO - jy;
- dzO = izO - jz;
- rsqO = dxO*dxO + dyO*dyO + dzO*dzO;
-
- /* Doing fast gmx_invsqrt */
- rinv1O = gmx_invsqrt(rsqO);
-
- /* O block */
- r1 = one/rinv1O;
- r1t = r1*tabscale;
- n0 = r1t;
- n1 = 12*n0;
- eps = r1t-n0;
- eps2 = eps*eps;
- nnn = n1;
- 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;
-
- pot[jnr] += VV*qO;
- poti += VV*qj;
-
- }
- pot[ii] += poti;
- }
-}
-
-static void low_calc_pot(FILE *log, int nl_type, t_forcerec *fr,
- rvec x[], t_mdatoms *mdatoms, real pot[])
-{
- t_nblist *nlist;
-
- nlist = &fr->nblists[0].nlist_sr[nl_type];
-
- c_tabpot(fr->nblists[0].table_elec_vdw.scale, fr->nblists[0].table_elec_vdw.data,
- nlist->nri, nlist->iinr,
- nlist->shift, nlist->jindex, nlist->jjnr,
- x[0], fr->epsfac, mdatoms->chargeA, pot, fr->shift_vec[0]);
-
- fprintf(log, "There were %d interactions\n", nlist->nrj);
-}
-
-void calc_pot(FILE *logf, t_commrec *cr,
- gmx_mtop_t *mtop,
- t_inputrec *inputrec, gmx_localtop_t *top, rvec x[],
- t_forcerec *fr, gmx_enerdata_t *enerd,
- t_mdatoms *mdatoms, real pot[], matrix box, t_graph *graph)
-{
- static t_nrnb nrnb;
- real lam[efptNR], dum[efptNR];
- rvec box_size;
- int i, m;
-
- /* Calc the force */
- fprintf(stderr, "Doing single force calculation...\n");
-
- if (inputrec->ePBC != epbcNONE)
- {
- calc_shifts(box, fr->shift_vec);
- }
-
- put_charge_groups_in_box(logf, 0, top->cgs.nr, fr->ePBC, box, &(top->cgs),
- x, fr->cg_cm);
- if (graph)
- {
- mk_mshift(logf, graph, fr->ePBC, box, x);
- }
- /* Do the actual neighbour searching and if twin range electrostatics
- * also do the calculation of long range forces and energies.
- */
- ns(logf, fr, x, box, &mtop->groups, &(inputrec->opts), top, mdatoms, cr,
- &nrnb, &lam[0], &dum[0], &enerd->grpp, TRUE, FALSE);
- for (m = 0; (m < DIM); m++)
- {
- box_size[m] = box[m][m];
- }
- for (i = 0; (i < mdatoms->nr); i++)
- {
- pot[i] = 0;
- }
- if (debug)
- {
- pr_rvecs(debug, 0, "x", x, mdatoms->nr);
- pr_rvecs(debug, 0, "cgcm", fr->cg_cm, top->cgs.nr);
- }
- /* electrostatics from any atom to atoms without LJ */
- low_calc_pot(logf, eNL_QQ, fr, x, mdatoms, pot);
- /* electrostatics from any atom to atoms without charge */
- low_calc_pot(logf, eNL_VDW, fr, x, mdatoms, pot);
- /* electrostatics from any atom to atoms with LJ */
- low_calc_pot(logf, eNL_VDWQQ, fr, x, mdatoms, pot);
-}
-
-FILE *init_calcpot(const char *log, const char *tpx, const char *table,
- gmx_mtop_t *mtop, gmx_localtop_t *top,
- t_inputrec *inputrec, t_commrec **cr,
- t_graph **graph, t_mdatoms **mdatoms,
- t_forcerec **fr,
- gmx_enerdata_t *enerd,
- real **pot,
- matrix box, rvec **x, const output_env_t oenv)
-{
- gmx_localtop_t *ltop;
- double t, t0;
- real lam[efptNR];
- int fep_state;
- gmx_bool bNEMD, bSA;
- int traj = 0, xtc_traj = 0;
- t_state *state;
- rvec mutot;
- t_nrnb nrnb;
- int m;
- rvec box_size;
- tensor force_vir, shake_vir;
- FILE *fplog;
-
- /* Initiate */
- *cr = init_par(NULL, NULL);
- gmx_log_open(log, *cr, FALSE, 0, &fplog);
-
- init_nrnb(&nrnb);
- snew(state, 1);
- read_tpx_state(tpx, inputrec, state, NULL, mtop);
- set_state_entries(state, inputrec, 1);
- if (fplog)
- {
- pr_inputrec(fplog, 0, "Input Parameters", inputrec, FALSE);
- }
-
-
-
- if (inputrec->efep)
- {
- fprintf(stderr, "WARNING: turning of free energy, will use lambda=0\n");
- inputrec->efep = 0;
- }
-
- clear_rvec(mutot);
- init_md(fplog, *cr, inputrec, oenv, &t, &t0, lam, &fep_state, NULL,
- &nrnb, mtop, NULL, -1, NULL, NULL, NULL,
- force_vir, shake_vir, mutot, &bSA, NULL, NULL, 0);
-
- init_enerdata(mtop->groups.grps[egcENER].nr, 0, enerd);
-
- ltop = gmx_mtop_generate_local_top(mtop, inputrec);
- *top = *ltop;
- sfree(ltop);
-
- *mdatoms = init_mdatoms(fplog, mtop, FALSE);
- atoms2md(mtop, inputrec, 0, NULL, 0, mtop->natoms, *mdatoms);
-
- if (inputrec->ePBC == epbcXYZ)
- {
- /* Calculate intramolecular shift vectors to make molecules whole again */
- *graph = mk_graph(fplog, &(top->idef), 0, mtop->natoms, FALSE, FALSE);
- mk_mshift(fplog, *graph, inputrec->ePBC, state->box, state->x);
- }
- else
- {
- *graph = NULL;
- }
-
- /* Turn off twin range if appropriate */
- inputrec->rvdw = inputrec->rcoulomb;
- inputrec->rlist = inputrec->rcoulomb;
- fprintf(stderr, "Using a coulomb cut-off of %g nm\n", inputrec->rcoulomb);
-
- /* Turn off free energy computation */
- inputrec->efep = 0;
-
- /* Set vanderwaals to shift, to force tables */
- inputrec->vdwtype = evdwSHIFT;
- inputrec->rvdw_switch = 0.0;
- inputrec->eDispCorr = edispcNO;
-
- /* Initiate forcerecord */
- *fr = mk_forcerec();
- init_forcerec(fplog, oenv, *fr, NULL, inputrec, mtop, *cr,
- state->box, FALSE, table, NULL, table, NULL, NULL, TRUE, -1);
-
- /* Remove periodicity */
- for (m = 0; (m < DIM); m++)
- {
- box_size[m] = state->box[m][m];
- }
- if (inputrec->ePBC != epbcNONE)
- {
- do_pbc_first(fplog, state->box, *fr, *graph, state->x);
- }
-
- copy_mat(state->box, box);
- *x = state->x;
- state->x = NULL;
- done_state(state);
- sfree(state);
-
- snew(*pot, mtop->natoms);
-
- return fplog;
-}
+++ /dev/null
-/*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * If you want to redistribute modifications, 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 www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
- *
- * And Hey:
- * Green Red Orange Magenta Azure Cyan Skyblue
- */
-
-extern FILE *init_calcpot(const char *log, const char *tpx, const char *table,
- gmx_mtop_t *mtop,
- gmx_localtop_t *top, t_inputrec *inputrec,
- t_commrec **cr,
- t_graph **graph, t_mdatoms **mdatoms,
- t_forcerec **fr,
- gmx_enerdata_t *enerdata,
- real **coulomb,
- matrix box, rvec **x, const output_env_t oenv);
-
-extern void calc_pot(FILE *logf, t_commrec *cr,
- gmx_mtop_t *mtop,
- t_inputrec *inputrec, gmx_localtop_t *top, rvec x[],
- t_forcerec *fr, gmx_enerdata_t *enerd,
- t_mdatoms *mdatoms, real coulomb[], matrix box, t_graph *graph);
-
-extern void write_pdb_coul();
-
-extern void delete_atom(gmx_localtop_t *top, int inr);
-/* Delete an atom from a topology */
-
-extern void replace_atom(gmx_localtop_t *top, int inr, char *anm, char *resnm,
- real q, real m, int type);
-/* Replace an atom in a topology by someting else */
}
else if (!scprev && !sc)
{
- gmx_fatal(FARGS, "There is no path from lambda=%g -> %g that is covered by foreign lambdas:\ncannot proceed with BAR.\nUse thermodynamic integration of dH/dl by calculating the averages of dH/dl\nwith g_analyze and integrating them.\nAlternatively, use the -extp option if (and only if) the Hamiltonian\ndepends linearly on lambda, which is NOT normally the case.\n", bl->prev->lambda, bl->lambda);
+ gmx_fatal(FARGS, "There is no path from lambda=%f -> %f that is covered by foreign lambdas:\ncannot proceed with BAR.\nUse thermodynamic integration of dH/dl by calculating the averages of dH/dl\nwith g_analyze and integrating them.\nAlternatively, use the -extp option if (and only if) the Hamiltonian\ndepends linearly on lambda, which is NOT normally the case.\n", bl->prev->lambda, bl->lambda);
}
/* normal delta H */
if (!scprev)
{
- gmx_fatal(FARGS, "Could not find a set for foreign lambda = %g\nin the files for lambda = %g", bl->lambda, bl->prev->lambda);
+ gmx_fatal(FARGS, "Could not find a set for foreign lambda = %f\nin the files for lambda = %f", bl->lambda, bl->prev->lambda);
}
if (!sc)
{
- gmx_fatal(FARGS, "Could not find a set for foreign lambda = %g\nin the files for lambda = %g", bl->prev->lambda, bl->lambda);
+ gmx_fatal(FARGS, "Could not find a set for foreign lambda = %f\nin the files for lambda = %f", bl->prev->lambda, bl->lambda);
}
br->a = scprev;
br->b = sc;
{
DG1 = 0.5*(DG0 + DG2);
- /*printf("Wfac1=%g, Wfac2=%g, beta=%g, DG1=%g\n",Wfac1,Wfac2,beta,
- DG1);*/
-
/* calculate the BAR averages */
dDG1 = 0.;
{
if (ba->temp <= 0)
{
- gmx_fatal(FARGS, "Found temperature of %g in file '%s'",
+ gmx_fatal(FARGS, "Found temperature of %f in file '%s'",
ba->temp, fn);
}
}
(blk->sub[1].nr < 1) )
{
gmx_fatal(FARGS,
- "Unexpected/corrupted block data in file %s around time %g.",
+ "Unexpected/corrupted block data in file %s around time %f.",
filename, start_time);
}
lambda_vec_print(foreign_lambda, buf, FALSE);
lambda_vec_print(s->foreign_lambda, buf2, FALSE);
fprintf(stderr, "Got foreign lambda=%s, expected: %s\n", buf, buf2);
- gmx_fatal(FARGS, "Corrupted data in file %s around t=%g.",
+ gmx_fatal(FARGS, "Corrupted data in file %s around t=%f.",
filename, start_time);
}
(blk->sub[1].nr < 2) )
{
gmx_fatal(FARGS,
- "Unexpected/corrupted block data in file %s around time %g",
+ "Unexpected/corrupted block data in file %s around time %f",
filename, start_time);
}
if (nhist > 2)
{
gmx_fatal(FARGS,
- "Unexpected/corrupted block data in file %s around time %g",
+ "Unexpected/corrupted block data in file %s around time %f",
filename, start_time);
}
if (fr->block[i].nsub < 2)
{
gmx_fatal(FARGS,
- "No lambda vector, but start_lambda=%g\n",
+ "No lambda vector, but start_lambda=%f\n",
old_start_lambda);
}
n_lambda_vec = fr->block[i].sub[1].ival[1];
/* check the native lambda */
if (!lambda_vec_same(&start_lambda, native_lambda) )
{
- gmx_fatal(FARGS, "Native lambda not constant in file %s: started at %g, and becomes %g at time %g",
+ gmx_fatal(FARGS, "Native lambda not constant in file %s: started at %f, and becomes %f at time %f",
fn, native_lambda, start_lambda, start_time);
}
/* check the number of samples against the previous number */
"[TT]*[tt] lam_B: the [GRK]lambda[grk] values for point B.[BR]",
"[TT]*[tt] DG: the free energy estimate.[BR]",
"[TT]*[tt] s_A: an estimate of the relative entropy of B in A.[BR]",
- "[TT]*[tt] s_A: an estimate of the relative entropy of A in B.[BR]",
+ "[TT]*[tt] s_B: an estimate of the relative entropy of A in B.[BR]",
"[TT]*[tt] stdev: an estimate expected per-sample standard deviation.[PAR]",
"The relative entropy of both states in each other's ensemble can be ",
#include "futil.h"
#include "maths.h"
#include "macros.h"
-#include "physics.h"
#include "vec.h"
#include "tpxio.h"
#include "mdrun.h"
-#include "calcpot.h"
#include "main.h"
#include "random.h"
#include "index.h"
static void insert_ion(int nsa, int *nwater,
gmx_bool bSet[], int repl[], atom_id index[],
- real pot[], rvec x[], t_pbc *pbc,
+ rvec x[], t_pbc *pbc,
int sign, int q, const char *ionname,
- t_mdatoms *mdatoms,
- real rmin, gmx_bool bRandom, int *seed)
+ t_atoms *atoms,
+ real rmin, int *seed)
{
- int i, ii, ei, owater, wlast, m, nw;
- real extr_e, poti, rmin2;
- rvec xei, dx;
- gmx_bool bSub = FALSE;
+ int i, ei,nw;
+ real rmin2;
+ rvec dx;
gmx_large_int_t maxrand;
ei = -1;
nw = *nwater;
maxrand = nw;
maxrand *= 1000;
- if (bRandom)
+
+ do
{
- do
- {
- ei = nw*rando(seed);
- maxrand--;
- }
- while (bSet[ei] && (maxrand > 0));
- if (bSet[ei])
- {
- gmx_fatal(FARGS, "No more replaceable solvent!");
- }
+ ei = nw*rando(seed);
+ maxrand--;
}
- else
+ while (bSet[ei] && (maxrand > 0));
+ if (bSet[ei])
{
- extr_e = 0;
- for (i = 0; (i < nw); i++)
- {
- if (!bSet[i])
- {
- ii = index[nsa*i];
- poti = pot[ii];
- if (q > 0)
- {
- if ((poti <= extr_e) || !bSub)
- {
- extr_e = poti;
- ei = i;
- bSub = TRUE;
- }
- }
- else
- {
- if ((poti >= extr_e) || !bSub)
- {
- extr_e = poti;
- ei = i;
- bSub = TRUE;
- }
- }
- }
- }
- if (ei == -1)
- {
- gmx_fatal(FARGS, "No more replaceable solvent!");
- }
+ gmx_fatal(FARGS, "No more replaceable solvent!");
}
+
fprintf(stderr, "Replacing solvent molecule %d (atom %d) with %s\n",
ei, index[nsa*ei], ionname);
/* Replace solvent molecule charges with ion charge */
bSet[ei] = TRUE;
repl[ei] = sign;
- mdatoms->chargeA[index[nsa*ei]] = q;
+
+ atoms->atom[index[nsa*ei]].q = q;
for (i = 1; i < nsa; i++)
{
- mdatoms->chargeA[index[nsa*ei+i]] = 0;
+ atoms->atom[index[nsa*ei+i]].q = 0;
}
/* Mark all solvent molecules within rmin as unavailable for substitution */
}
}
+
static char *aname(const char *mname)
{
char *str;
int gmx_genion(int argc, char *argv[])
{
const char *desc[] = {
- "[TT]genion[tt] replaces solvent molecules by monoatomic ions at",
- "the position of the first atoms with the most favorable electrostatic",
- "potential or at random. The potential is calculated on all atoms, using",
- "normal GROMACS particle-based methods (in contrast to other methods",
- "based on solving the Poisson-Boltzmann equation).",
- "The potential is recalculated after every ion insertion.",
- "If specified in the run input file, a reaction field, shift function",
- "or user function can be used. For the user function a table file",
- "can be specified with the option [TT]-table[tt].",
+ "[TT]genion[tt] randomly replaces solvent molecules with monoatomic ions.",
"The group of solvent molecules should be continuous and all molecules",
"should have the same number of atoms.",
"The user should add the ion molecules to the topology file or use",
"either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
"[PAR]Ions which can have multiple charge states get the multiplicity",
"added, without sign, for the uncommon states only.[PAR]",
- "With the option [TT]-pot[tt] the potential can be written as B-factors",
- "in a [TT].pdb[tt] file (for visualisation using e.g. Rasmol).",
- "The unit of the potential is 1000 kJ/(mol e), the scaling be changed",
- "with the [TT]-scale[tt] option.[PAR]",
"For larger ions, e.g. sulfate we recommended using [TT]genbox[tt]."
};
const char *bugs[] = {
- "Calculation of the potential is not reliable, therefore the [TT]-random[tt] option is now turned on by default.",
- "If you specify a salt concentration existing ions are not taken into account. In effect you therefore specify the amount of salt to be added."
+ "If you specify a salt concentration existing ions are not taken into "
+ "account. In effect you therefore specify the amount of salt to be added.",
};
static int p_num = 0, n_num = 0, p_q = 1, n_q = -1;
static const char *p_name = "NA", *n_name = "CL";
- static real rmin = 0.6, scale = 0.001, conc = 0;
+ static real rmin = 0.6, conc = 0;
static int seed = 1993;
- static gmx_bool bRandom = TRUE, bNeutral = FALSE;
+ static gmx_bool bNeutral = FALSE;
static t_pargs pa[] = {
{ "-np", FALSE, etINT, {&p_num}, "Number of positive ions" },
{ "-pname", FALSE, etSTR, {&p_name}, "Name of the positive ion" },
{ "-nname", FALSE, etSTR, {&n_name}, "Name of the negative ion" },
{ "-nq", FALSE, etINT, {&n_q}, "Charge of the negative ion" },
{ "-rmin", FALSE, etREAL, {&rmin}, "Minimum distance between ions" },
- { "-random", FALSE, etBOOL, {&bRandom}, "Use random placement of ions instead of based on potential. The rmin option should still work" },
{ "-seed", FALSE, etINT, {&seed}, "Seed for random number generator" },
- { "-scale", FALSE, etREAL, {&scale}, "Scaling factor for the potential for [TT]-pot[tt]" },
{ "-conc", FALSE, etREAL, {&conc},
"Specify salt concentration (mol/liter). This will add sufficient ions to reach up to the specified concentration as computed from the volume of the cell in the input [TT].tpr[tt] file. Overrides the [TT]-np[tt] and [TT]-nn[tt] options." },
{ "-neutral", FALSE, etBOOL, {&bNeutral}, "This option will add enough ions to neutralize the system. These ions are added on top of those specified with [TT]-np[tt]/[TT]-nn[tt] or [TT]-conc[tt]. "}
};
- gmx_mtop_t *mtop;
- gmx_localtop_t *top;
- t_inputrec inputrec;
- t_commrec *cr;
- t_mdatoms *mdatoms;
- gmx_enerdata_t enerd;
- t_graph *graph;
- t_forcerec *fr;
+ t_topology top;
rvec *x, *v;
- real *pot, vol, qtot;
+ real vol, qtot;
matrix box;
t_atoms atoms;
t_pbc pbc;
- int *repl;
+ int *repl, ePBC;
atom_id *index;
- char *grpname;
- gmx_bool *bSet, bPDB;
+ char *grpname, title[STRLEN];
+ gmx_bool *bSet;
int i, nw, nwa, nsa, nsalt, iqtot;
- FILE *fplog;
output_env_t oenv;
t_filenm fnm[] = {
{ efTPX, NULL, NULL, ffREAD },
- { efXVG, "-table", "table", ffOPTRD },
{ efNDX, NULL, NULL, ffOPTRD },
{ efSTO, "-o", NULL, ffWRITE },
- { efLOG, "-g", "genion", ffWRITE },
- { efPDB, "-pot", "pot", ffOPTWR },
{ efTOP, "-p", "topol", ffOPTRW }
};
#define NFILE asize(fnm)
parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
asize(desc), desc, asize(bugs), bugs, &oenv);
- bPDB = ftp2bSet(efPDB, NFILE, fnm);
- if (bRandom && bPDB)
- {
- fprintf(stderr, "Not computing potential with random option!\n");
- bPDB = FALSE;
- }
/* Check input for something sensible */
if ((p_num < 0) || (n_num < 0))
gmx_fatal(FARGS, "Negative number of ions to add?");
}
- snew(mtop, 1);
- snew(top, 1);
- fplog = init_calcpot(ftp2fn(efLOG, NFILE, fnm), ftp2fn(efTPX, NFILE, fnm),
- opt2fn("-table", NFILE, fnm), mtop, top, &inputrec, &cr,
- &graph, &mdatoms, &fr, &enerd, &pot, box, &x, oenv);
+ if (conc > 0 && (p_num > 0 || n_num > 0))
+ {
+ fprintf(stderr, "WARNING: -conc specified, overriding -nn and -np.\n");
+ }
- atoms = gmx_mtop_global_atoms(mtop);
+ /* Read atom positions and charges */
+ read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, &x, &v, box, FALSE);
+ atoms = top.atoms;
+ /* Compute total charge */
qtot = 0;
for (i = 0; (i < atoms.nr); i++)
{
if ((p_num == 0) && (n_num == 0))
{
- if (!bPDB)
- {
- fprintf(stderr, "No ions to add and no potential to calculate.\n");
- exit(0);
- }
- nw = 0;
- nsa = 0; /* to keep gcc happy */
+ fprintf(stderr, "No ions to add.\n");
+ exit(0);
}
else
{
snew(v, atoms.nr);
snew(atoms.pdbinfo, atoms.nr);
- set_pbc(&pbc, inputrec.ePBC, box);
+ set_pbc(&pbc, ePBC, box);
/* Now loop over the ions that have to be placed */
- do
+ while (p_num-- > 0)
{
- if (!bRandom)
- {
- calc_pot(fplog, cr, mtop, &inputrec, top, x, fr, &enerd, mdatoms, pot, box, graph);
- if (bPDB || debug)
- {
- char buf[STRLEN];
-
- if (debug)
- {
- sprintf(buf, "%d_%s", p_num+n_num, ftp2fn(efPDB, NFILE, fnm));
- }
- else
- {
- strcpy(buf, ftp2fn(efPDB, NFILE, fnm));
- }
- for (i = 0; (i < atoms.nr); i++)
- {
- atoms.pdbinfo[i].bfac = pot[i]*scale;
- }
- write_sto_conf(buf, "Potential calculated by genion",
- &atoms, x, v, inputrec.ePBC, box);
- bPDB = FALSE;
- }
- }
- if ((p_num > 0) && (p_num >= n_num))
- {
- insert_ion(nsa, &nw, bSet, repl, index, pot, x, &pbc,
- 1, p_q, p_name, mdatoms, rmin, bRandom, &seed);
- p_num--;
- }
- else if (n_num > 0)
- {
- insert_ion(nsa, &nw, bSet, repl, index, pot, x, &pbc,
- -1, n_q, n_name, mdatoms, rmin, bRandom, &seed);
- n_num--;
- }
+ insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
+ 1, p_q, p_name, &atoms, rmin, &seed);
+ }
+ while (n_num-- > 0)
+ {
+ insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
+ -1, n_q, n_name, &atoms, rmin, &seed);
}
- while (p_num+n_num > 0);
fprintf(stderr, "\n");
if (nw)
sfree(atoms.pdbinfo);
atoms.pdbinfo = NULL;
- write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *mtop->name, &atoms, x, NULL,
- inputrec.ePBC, box);
+ write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, NULL, ePBC,
+ box);
thanx(stderr);
- gmx_log_close(fplog);
-
return 0;
}
/* Dont change szProgram arbitrarily - it must be argv[0], i.e. the
* name of a file. Otherwise, we won't be able to find the library dir.
*/
+
#define NCR (int)asize(CopyrightText)
/* TODO: Is this exception still needed? */
#ifdef GMX_FAHCORE
ster_print(out, GromacsVersion());
fprintf(out, "\n");
+ if (getenv("GMX_NO_CREDITS"))
+ {
+ return;
+ }
+
/* fprintf(out,"\n");*/
/* sp_print(out,"PLEASE NOTE: THIS IS A BETA VERSION\n");
* reserved for that part.
* This routine does internal renumbering so we get continuous indices, and also
* decodes the actual number of packages,cores-per-package and hwthreads-per-core.
+ * Returns: 0 on success, non-zero on failure.
*/
-static void
+static int
cpuid_x86_decode_apic_id(gmx_cpuid_t cpuid, int *apic_id, int core_bits, int hwthread_bits)
{
int i, idx;
cpuid->ncores_per_package = cpuid_renumber_elements(cpuid->core_id, cpuid->nproc);
cpuid->nhwthreads_per_core = cpuid_renumber_elements(cpuid->hwthread_id, cpuid->nproc);
+ /* now check for consistency */
+ if ( (cpuid->npackages * cpuid->ncores_per_package *
+ cpuid->nhwthreads_per_core) != cpuid->nproc )
+ {
+ /* the packages/cores-per-package/hwthreads-per-core counts are
+ inconsistent. */
+ return -1;
+ }
+
/* Create a locality order array, i.e. first all resources in package0, which in turn
* are sorted so we first have all resources in core0, where threads are sorted in order, etc.
*/
+
for (i = 0; i < cpuid->nproc; i++)
{
idx = (cpuid->package_id[i]*cpuid->ncores_per_package + cpuid->core_id[i])*cpuid->nhwthreads_per_core + cpuid->hwthread_id[i];
cpuid->locality_order[idx] = i;
}
+ return 0;
}
static int
cpuid_check_amd_x86(gmx_cpuid_t cpuid)
{
- int max_stdfn, max_extfn;
+ int max_stdfn, max_extfn, ret;
unsigned int eax, ebx, ecx, edx;
int hwthread_bits, core_bits;
int * apic_id;
;
}
}
- cpuid_x86_decode_apic_id(cpuid, apic_id, core_bits, hwthread_bits);
- cpuid->have_cpu_topology = 1;
+ ret = cpuid_x86_decode_apic_id(cpuid, apic_id, core_bits,
+ hwthread_bits);
+ cpuid->have_cpu_topology = (ret == 0);
#endif
}
return 0;
static int
cpuid_check_intel_x86(gmx_cpuid_t cpuid)
{
- unsigned int max_stdfn, max_extfn;
+ unsigned int max_stdfn, max_extfn, ret;
unsigned int eax, ebx, ecx, edx;
unsigned int max_logical_cores, max_physical_cores;
int hwthread_bits, core_bits;
hwthread_bits = eax & 0x1F;
execute_x86cpuid(0xB, 1, &eax, &ebx, &ecx, &edx);
core_bits = (eax & 0x1F) - hwthread_bits;
- cpuid_x86_decode_apic_id(cpuid, apic_id, core_bits, hwthread_bits);
- cpuid->have_cpu_topology = 1;
+ ret = cpuid_x86_decode_apic_id(cpuid, apic_id, core_bits,
+ hwthread_bits);
+ cpuid->have_cpu_topology = (ret == 0);
#endif
}
return 0;
return;
}
-#ifndef __APPLE__
/* If the tMPI thread affinity setting is not supported encourage the user
* to report it as it's either a bug or an exotic platform which we might
* want to support. */
if (tMPI_Thread_setaffinity_support() != TMPI_SETAFFINITY_SUPPORT_YES)
{
+ /* we know Mac OS doesn't support setting thread affinity, so there's
+ no point in warning the user in that case. In any other case
+ the user might be able to do something about it. */
+#ifndef __APPLE__
md_print_warn(NULL, fplog,
"Can not set thread affinities on the current platform. On NUMA systems this\n"
"can cause performance degradation. If you think your platform should support\n"
"setting affinities, contact the GROMACS developers.");
+#endif /* __APPLE__ */
return;
}
-#endif /* __APPLE__ */
/* threads on this MPI process or TMPI thread */
if (cr->duty & DUTY_PP)
#endif
const __m128d signmask = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF) );
- const __m128i signbit_epi32 = _mm_set1_epi32(0x80000000);
const __m128d tabscale = _mm_set1_pd(32.0/M_PI);
const __m128d invtabscale0 = _mm_set1_pd(9.81747508049011230469e-02);
__m128d sinpoint, cospoint;
__m128d xsign, ssign, csign;
__m128i imask, sswapsign, cswapsign;
- __m128d minusone;
xsign = _mm_andnot_pd(signmask, x);
xabs = _mm_and_pd(x, signmask);
const __m128d limit1 = _mm_set1_pd(0.625);
const __m128d limit2 = _mm_set1_pd(1e-8);
const __m128d one = _mm_set1_pd(1.0);
- const __m128d halfpi = _mm_set1_pd(M_PI/2.0);
const __m128d quarterpi = _mm_set1_pd(M_PI/4.0);
const __m128d morebits = _mm_set1_pd(6.123233995736765886130e-17);
__m128d sign;
__m128d mask;
__m128d xabs;
- __m128d zz, ww, z, q, w, y, zz2, ww2;
+ __m128d zz, ww, z, q, w, zz2, ww2;
__m128d PA, PB;
__m128d QA, QB;
__m128d RA, RB;
static __m128d
gmx_mm_acos_pd(__m128d x)
{
- const __m128d signmask = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF) );
const __m128d one = _mm_set1_pd(1.0);
const __m128d half = _mm_set1_pd(0.5);
- const __m128d pi = _mm_set1_pd(M_PI);
const __m128d quarterpi0 = _mm_set1_pd(7.85398163397448309616e-1);
const __m128d quarterpi1 = _mm_set1_pd(6.123233995736765886130e-17);
const __m128 CE_0 = _mm_set1_ps(1.3429983063133937f);
const __m128 CE_1 = _mm_set1_ps(1.807420826584643f);
- __m128 fexp, fexp1;
+ __m128 fexp;
__m128i iexp;
__m128 mask;
- __m128 x1, x2;
+ __m128 x2;
__m128 y;
__m128 pA, pB, pC, pD, pE, tB, tC, tD, tE;
{
zones->cg_range[i] = dd->ncg_home;
}
+ /* zone_ncg1[0] should always be equal to ncg_home */
+ dd->comm->zone_ncg1[0] = dd->ncg_home;
}
static void rebuild_cgindex(gmx_domdec_t *dd,
t_block *cgs_gl;
gmx_large_int_t step_pcoupl;
rvec cell_ns_x0, cell_ns_x1;
- int i, j, n, cg0 = 0, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
+ int i, j, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
gmx_bool bRedist, bSortCG, bResortAll;
ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
{
/* Clear the old state */
clear_dd_indices(dd, 0, 0);
+ ncgindex_set = 0;
set_ddbox(dd, bMasterState, cr, ir, state_global->box,
TRUE, cgs_gl, state_global->x, &ddbox);
inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
-
- cg0 = 0;
}
else if (state_local->ddp_count != dd->ddp_count)
{
/* Build the new indices */
rebuild_cgindex(dd, cgs_gl->index, state_local);
make_dd_indices(dd, cgs_gl->index, 0);
+ ncgindex_set = dd->ncg_home;
if (fr->cutoff_scheme == ecutsGROUP)
{
/* Clear the non-home indices */
clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
+ ncgindex_set = 0;
/* Avoid global communication for dim's without pbc and -gcom */
if (!bNStGlobalComm)
dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
state_local, f, fr, mdatoms,
- !bSortCG, nrnb, &cg0, &ncg_moved);
+ !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
}
dd_sort_state(dd, ir->ePBC, fr->cg_cm, fr, state_local,
bResortAll ? -1 : ncg_home_old);
/* Rebuild all the indices */
- cg0 = 0;
ga2la_clear(dd->ga2la);
+ ncgindex_set = 0;
wallcycle_sub_stop(wcycle, ewcsDD_GRID);
}
setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
/* Set the indices */
- make_dd_indices(dd, cgs_gl->index, cg0);
+ make_dd_indices(dd, cgs_gl->index, ncgindex_set);
/* Set the charge group boundaries for neighbor searching */
set_cg_boundaries(&comm->zones);
}
/* < 0 means stop at next step, > 0 means stop at next NS step */
- if ( (gs.set[eglsSTOPCOND] < 0 ) ||
- ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist == 0)) )
+ if ( (gs.set[eglsSTOPCOND] < 0) ||
+ ( (gs.set[eglsSTOPCOND] > 0) && (bNStList || ir->nstlist == 0) ) )
{
bLastStep = TRUE;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff