From: Mark Abraham Date: Fri, 26 Oct 2018 13:47:06 +0000 (+0200) Subject: Merge branch release-2018 into release-2019 X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=fb755982116311bee2d1f9770678d221112e1ddb;p=alexxy%2Fgromacs.git Merge branch release-2018 into release-2019 Restored some text about gmx_enerdata_t fields that got list when master branch moved some content from forcerec.h into the new enerdata.h. Change-Id: I45b956e02d6713190fca456ad0a5948db5001737 --- fb755982116311bee2d1f9770678d221112e1ddb diff --cc docs/CMakeLists.txt index f7af04cd47,22e8abeaec..6b72f9d299 --- a/docs/CMakeLists.txt +++ b/docs/CMakeLists.txt @@@ -353,27 -150,28 +353,28 @@@ if (SPHINX_FOUND release-notes/2016/major/removed-features.rst release-notes/2016/major/miscellaneous.rst release-notes/older/index.rst - user-guide/index.rst + # the entry for user-guide/index.rst should not appear here, + # as it will be included conditionally further down depending on + # if the documentation will be build with the full reference + # manual or without. + user-guide/cmdline.rst user-guide/cutoff-schemes.rst - user-guide/getting-started.rst - user-guide/force-fields.rst + user-guide/deprecation-policy.rst + user-guide/environment-variables.rst user-guide/faq.rst - user-guide/flow.rst user-guide/floating-point.rst - user-guide/system-preparation.rst + user-guide/flow.rst + user-guide/force-fields.rst + user-guide/getting-started.rst + user-guide/index.rst user-guide/managing-simulations.rst + user-guide/mdp-options.rst user-guide/mdrun-features.rst user-guide/mdrun-performance.rst - user-guide/mdp-options.rst user-guide/run-time-errors.rst - user-guide/file-formats.rst - user-guide/cmdline.rst - user-guide/environment-variables.rst + user-guide/security.rst + user-guide/system-preparation.rst user-guide/terminology.rst - user-guide/plotje.gif - user-guide/xvgr.gif - conf.py - links.dat ) include(SphinxMacros.cmake) diff --cc docs/user-guide/index.rst index 9330890fb1,5a352a5242..ad509c9eca --- a/docs/user-guide/index.rst +++ b/docs/user-guide/index.rst @@@ -41,4 -37,4 +41,5 @@@ For background on algorithms and implem terminology environment-variables floating-point + security + deprecation-policy diff --cc src/gromacs/ewald/pme.cpp index 91675f18e9,8b861b3163..717eadbfad --- a/src/gromacs/ewald/pme.cpp +++ b/src/gromacs/ewald/pme.cpp @@@ -1125,10 -1013,20 +1125,20 @@@ int gmx_pme_do(struct gmx_pme_t *pme gmx_bool bFirst, bDoSplines; int fep_state; int fep_states_lj = pme->bFEP_lj ? 2 : 1; - const gmx_bool bCalcEnerVir = flags & GMX_PME_CALC_ENER_VIR; - const gmx_bool bBackFFT = flags & (GMX_PME_CALC_F | GMX_PME_CALC_POT); - const gmx_bool bCalcF = flags & GMX_PME_CALC_F; + const gmx_bool bCalcEnerVir = (flags & GMX_PME_CALC_ENER_VIR) != 0; + const gmx_bool bBackFFT = (flags & (GMX_PME_CALC_F | GMX_PME_CALC_POT)) != 0; + const gmx_bool bCalcF = (flags & GMX_PME_CALC_F) != 0; + /* We could be passing lambda!=1 while no q or LJ is actually perturbed */ + if (!pme->bFEP_q) + { + lambda_q = 1; + } + if (!pme->bFEP_lj) + { + lambda_lj = 1; + } + assert(pme->nnodes > 0); assert(pme->nnodes == 1 || pme->ndecompdim > 0); diff --cc src/gromacs/mdlib/force.cpp index 93e623f128,c939c4e8a7..5e089e79f6 --- a/src/gromacs/mdlib/force.cpp +++ b/src/gromacs/mdlib/force.cpp @@@ -772,11 -791,30 +763,30 @@@ void sum_dhdl(gmx_enerdata_t *enerd, gm current lambda, because the contributions to the current lambda are automatically zeroed */ + double &enerpart_lambda = enerd->enerpart_lambda[i + 1]; + - for (size_t j = 0; j < lambda.size(); j++) + for (gmx::index j = 0; j < lambda.size(); j++) { /* Note that this loop is over all dhdl components, not just the separated ones */ - dlam = (fepvals->all_lambda[j][i] - lambda[j]); - enerd->enerpart_lambda[i+1] += dlam*enerd->dvdl_lin[j]; + const double dlam = fepvals->all_lambda[j][i] - lambda[j]; + + enerpart_lambda += dlam*enerd->dvdl_lin[j]; + + /* Constraints can not be evaluated at foreign lambdas, so we add + * a linear extrapolation. This is an approximation, but usually + * quite accurate since constraints change little between lambdas. + */ + if ((j == efptBONDED && fepvals->separate_dvdl[efptBONDED]) || + (j == efptFEP && !fepvals->separate_dvdl[efptBONDED])) + { + enerpart_lambda += dlam*enerd->term[F_DVDL_CONSTR]; + } + + if (j == efptMASS) + { + enerpart_lambda += dlam*enerd->term[F_DKDL]; + } + if (debug) { fprintf(debug, "enerdiff lam %g: (%15s), non-linear %f linear %f*%f\n", diff --cc src/gromacs/mdlib/membed.cpp index ee1b0c3b71,0000000000..c564756966 mode 100644,000000..100644 --- a/src/gromacs/mdlib/membed.cpp +++ b/src/gromacs/mdlib/membed.cpp @@@ -1,1321 -1,0 +1,1321 @@@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2010,2011,2012,2013,2014,2015,2016,2017,2018, 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 "membed.h" + +#include +#include + +#include "gromacs/commandline/filenm.h" +#include "gromacs/fileio/readinp.h" +#include "gromacs/fileio/warninp.h" +#include "gromacs/gmxlib/network.h" +#include "gromacs/math/vec.h" +#include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/inputrec.h" +#include "gromacs/mdtypes/md_enums.h" +#include "gromacs/mdtypes/state.h" +#include "gromacs/pbcutil/pbc.h" +#include "gromacs/topology/index.h" +#include "gromacs/topology/mtop_lookup.h" +#include "gromacs/topology/mtop_util.h" +#include "gromacs/topology/topology.h" +#include "gromacs/utility/cstringutil.h" +#include "gromacs/utility/exceptions.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/filestream.h" +#include "gromacs/utility/futil.h" +#include "gromacs/utility/smalloc.h" + +/* information about scaling center */ +typedef struct { + rvec xmin; /* smallest coordinates of all embedded molecules */ + rvec xmax; /* largest coordinates of all embedded molecules */ + rvec *geom_cent; /* scaling center of each independent molecule to embed */ + int pieces; /* number of molecules to embed independently */ + int *nidx; /* n atoms for every independent embedded molecule (index in subindex) */ + int **subindex; /* atomids for independent molecule * + * atoms of piece i run from subindex[i][0] to subindex[i][nidx[i]] */ +} pos_ins_t; + +/* variables needed in do_md */ +struct gmx_membed_t { + int it_xy; /* number of iterations (steps) used to grow something in the xy-plane */ + int it_z; /* same, but for z */ + real xy_step; /* stepsize used during resize in xy-plane */ + real z_step; /* same, but in z */ + rvec fac; /* initial scaling of the molecule to grow into the membrane */ + rvec *r_ins; /* final coordinates of the molecule to grow */ + pos_ins_t *pos_ins; /* scaling center for each piece to embed */ +}; + +/* membrane related variables */ +typedef struct { + char *name; /* name of index group to embed molecule into (usually membrane) */ + t_block mem_at; /* list all atoms in membrane */ + int nmol; /* number of membrane molecules overlapping with the molecule to embed */ + int *mol_id; /* list of molecules in membrane that overlap with the molecule to embed */ + real lip_area; /* average area per lipid in membrane (only correct for homogeneous bilayers)*/ + real zmin; /* minimum z coordinate of membrane */ + real zmax; /* maximum z coordinate of membrane */ + real zmed; /* median z coordinate of membrane */ +} mem_t; + +/* Lists all molecules in the membrane that overlap with the molecule to be embedded. * + * These will then be removed from the system */ +typedef struct { + int nr; /* number of molecules to remove */ + int *mol; /* list of molecule ids to remove */ + int *block; /* id of the molblock that the molecule to remove is part of */ +} rm_t; + +/* Get the global molecule id, and the corresponding molecule type and id of the * + * molblock from the global atom nr. */ +static int get_mol_id(int at, gmx_mtop_t *mtop, int *type, int *block) +{ + int mol_id = 0; + int i; + int atnr_mol; + + *block = 0; + mtopGetMolblockIndex(mtop, at, block, &mol_id, &atnr_mol); + for (i = 0; i < *block; i++) + { + mol_id += mtop->molblock[i].nmol; + } + *type = mtop->molblock[*block].type; + + return mol_id; +} + +/* Get the id of the molblock from a global molecule id */ +static int get_molblock(int mol_id, const std::vector &mblock) +{ + int nmol = 0; + + for (size_t i = 0; i < mblock.size(); i++) + { + nmol += mblock[i].nmol; + if (mol_id < nmol) + { + return i; + } + } + + gmx_fatal(FARGS, "mol_id %d larger than total number of molecules %d.\n", mol_id, nmol); +} + +/* Get a list of all the molecule types that are present in a group of atoms. * + * Because all interaction within the group to embed are removed on the topology * + * level, if the same molecule type is found in another part of the system, these * + * would also be affected. Therefore we have to check if the embedded and rest group * + * share common molecule types. If so, membed will stop with an error. */ +static int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist) +{ + int i, j, nr; + int type = 0, block = 0; + gmx_bool bNEW; + + nr = 0; + snew(tlist->index, at->nr); + for (i = 0; i < at->nr; i++) + { + bNEW = TRUE; + get_mol_id(at->index[i], mtop, &type, &block); + for (j = 0; j < nr; j++) + { + if (tlist->index[j] == type) + { + bNEW = FALSE; + } + } + + if (bNEW) + { + tlist->index[nr] = type; + nr++; + } + } + srenew(tlist->index, nr); + return nr; +} + +/* Do the actual check of the molecule types between embedded and rest group */ +static void check_types(t_block *ins_at, t_block *rest_at, gmx_mtop_t *mtop) +{ + t_block *ins_mtype, *rest_mtype; + int i, j; + + snew(ins_mtype, 1); + snew(rest_mtype, 1); + ins_mtype->nr = get_mtype_list(ins_at, mtop, ins_mtype ); + rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype); + + for (i = 0; i < ins_mtype->nr; i++) + { + for (j = 0; j < rest_mtype->nr; j++) + { + if (ins_mtype->index[i] == rest_mtype->index[j]) + { + gmx_fatal(FARGS, "Moleculetype %s is found both in the group to insert and the rest of the system.\n" + "1. Your *.ndx and *.top do not match\n" + "2. You are inserting some molecules of type %s (for example xray-solvent), while\n" + "the same moleculetype is also used in the rest of the system (solvent box). Because\n" + "we need to exclude all interactions between the atoms in the group to\n" + "insert, the same moleculetype can not be used in both groups. Change the\n" + "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n" + "an appropriate *.itp file", *(mtop->moltype[rest_mtype->index[j]].name), + *(mtop->moltype[rest_mtype->index[j]].name), *(mtop->moltype[rest_mtype->index[j]].name)); + } + } + } + + done_block(ins_mtype); + done_block(rest_mtype); + sfree(ins_mtype); + sfree(rest_mtype); +} + +static void get_input(const char *membed_input, real *xy_fac, real *xy_max, real *z_fac, real *z_max, + int *it_xy, int *it_z, real *probe_rad, int *low_up_rm, int *maxwarn, + int *pieces, gmx_bool *bALLOW_ASYMMETRY) +{ + warninp_t wi; + std::vector inp; + + wi = init_warning(TRUE, 0); + + { + gmx::TextInputFile stream(membed_input); + inp = read_inpfile(&stream, membed_input, wi); + stream.close(); + } + *it_xy = get_eint(&inp, "nxy", 1000, wi); + *it_z = get_eint(&inp, "nz", 0, wi); + *xy_fac = get_ereal(&inp, "xyinit", 0.5, wi); + *xy_max = get_ereal(&inp, "xyend", 1.0, wi); + *z_fac = get_ereal(&inp, "zinit", 1.0, wi); + *z_max = get_ereal(&inp, "zend", 1.0, wi); + *probe_rad = get_ereal(&inp, "rad", 0.22, wi); + *low_up_rm = get_eint(&inp, "ndiff", 0, wi); + *maxwarn = get_eint(&inp, "maxwarn", 0, wi); + *pieces = get_eint(&inp, "pieces", 1, wi); + const char *yesno_names[BOOL_NR+1] = + { + "no", "yes", nullptr + }; + *bALLOW_ASYMMETRY = (get_eeenum(&inp, "asymmetry", yesno_names, wi) != 0); + + check_warning_error(wi, FARGS); + { + gmx::TextOutputFile stream(membed_input); + write_inpfile(&stream, membed_input, &inp, FALSE, WriteMdpHeader::yes, wi); + stream.close(); + } + done_warning(wi, FARGS); +} + +/* Obtain the maximum and minimum coordinates of the group to be embedded */ +static int init_ins_at(t_block *ins_at, t_block *rest_at, t_state *state, pos_ins_t *pos_ins, + gmx_groups_t *groups, int ins_grp_id, real xy_max) +{ + int i, gid, c = 0; + real xmin, xmax, ymin, ymax, zmin, zmax; + const real min_memthick = 6.0; /* minimum thickness of the bilayer that will be used to * + * determine the overlap between molecule to embed and membrane */ + const real fac_inp_size = 1.000001; /* scaling factor to obtain input_size + 0.000001 (comparing reals) */ + snew(rest_at->index, state->natoms); + auto x = makeArrayRef(state->x); + + xmin = xmax = x[ins_at->index[0]][XX]; + ymin = ymax = x[ins_at->index[0]][YY]; + zmin = zmax = x[ins_at->index[0]][ZZ]; + + for (i = 0; i < state->natoms; i++) + { + gid = groups->grpnr[egcFREEZE][i]; + if (groups->grps[egcFREEZE].nm_ind[gid] == ins_grp_id) + { + xmin = std::min(xmin, x[i][XX]); + xmax = std::max(xmax, x[i][XX]); + ymin = std::min(ymin, x[i][YY]); + ymax = std::max(ymax, x[i][YY]); + zmin = std::min(zmin, x[i][ZZ]); + zmax = std::max(zmax, x[i][ZZ]); + } + else + { + rest_at->index[c] = i; + c++; + } + } + + rest_at->nr = c; + srenew(rest_at->index, c); + + if (xy_max > fac_inp_size) + { + pos_ins->xmin[XX] = xmin-((xmax-xmin)*xy_max-(xmax-xmin))/2; + pos_ins->xmin[YY] = ymin-((ymax-ymin)*xy_max-(ymax-ymin))/2; + + pos_ins->xmax[XX] = xmax+((xmax-xmin)*xy_max-(xmax-xmin))/2; + pos_ins->xmax[YY] = ymax+((ymax-ymin)*xy_max-(ymax-ymin))/2; + } + else + { + pos_ins->xmin[XX] = xmin; + pos_ins->xmin[YY] = ymin; + + pos_ins->xmax[XX] = xmax; + pos_ins->xmax[YY] = ymax; + } + + if ( (zmax-zmin) < min_memthick) + { + pos_ins->xmin[ZZ] = zmin+(zmax-zmin)/2.0-0.5*min_memthick; + pos_ins->xmax[ZZ] = zmin+(zmax-zmin)/2.0+0.5*min_memthick; + } + else + { + pos_ins->xmin[ZZ] = zmin; + pos_ins->xmax[ZZ] = zmax; + } + + return c; +} + +/* Estimate the area of the embedded molecule by projecting all coordinates on a grid in the * + * xy-plane and counting the number of occupied grid points */ +static real est_prot_area(pos_ins_t *pos_ins, rvec *r, t_block *ins_at, mem_t *mem_p) +{ + real x, y, dx = 0.15, dy = 0.15, area = 0.0; + real add, memmin, memmax; + int c, at; + + /* min and max membrane coordinate are altered to reduce the influence of the * + * boundary region */ + memmin = mem_p->zmin+0.1*(mem_p->zmax-mem_p->zmin); + memmax = mem_p->zmax-0.1*(mem_p->zmax-mem_p->zmin); + + //NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter) + for (x = pos_ins->xmin[XX]; x < pos_ins->xmax[XX]; x += dx) + { + //NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter) + for (y = pos_ins->xmin[YY]; y < pos_ins->xmax[YY]; y += dy) + { + c = 0; + add = 0.0; + do + { + at = ins_at->index[c]; + if ( (r[at][XX] >= x) && (r[at][XX] < x+dx) && + (r[at][YY] >= y) && (r[at][YY] < y+dy) && + (r[at][ZZ] > memmin) && (r[at][ZZ] < memmax) ) + { + add = 1.0; + } + c++; + } + while ( (c < ins_at->nr) && (add < 0.5) ); + area += add; + } + } + area = area*dx*dy; + + return area; +} + +static int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins) +{ + int i, j, at, mol, nmol, nmolbox, count; + t_block *mem_a; + real z, zmin, zmax, mem_area; + gmx_bool bNew; + int *mol_id; + int type = 0, block = 0; + + nmol = count = 0; + mem_a = &(mem_p->mem_at); + snew(mol_id, mem_a->nr); + zmin = pos_ins->xmax[ZZ]; + zmax = pos_ins->xmin[ZZ]; + for (i = 0; i < mem_a->nr; i++) + { + at = mem_a->index[i]; + if ( (r[at][XX] > pos_ins->xmin[XX]) && (r[at][XX] < pos_ins->xmax[XX]) && + (r[at][YY] > pos_ins->xmin[YY]) && (r[at][YY] < pos_ins->xmax[YY]) && + (r[at][ZZ] > pos_ins->xmin[ZZ]) && (r[at][ZZ] < pos_ins->xmax[ZZ]) ) + { + mol = get_mol_id(at, mtop, &type, &block); + bNew = TRUE; + for (j = 0; j < nmol; j++) + { + if (mol == mol_id[j]) + { + bNew = FALSE; + } + } + + if (bNew) + { + mol_id[nmol] = mol; + nmol++; + } + + z = r[at][ZZ]; + if (z < zmin) + { + zmin = z; + } + + if (z > zmax) + { + zmax = z; + } + + count++; + } + } + + mem_p->nmol = nmol; + srenew(mol_id, nmol); + mem_p->mol_id = mol_id; + + if ((zmax-zmin) > (box[ZZ][ZZ]-0.5)) + { + gmx_fatal(FARGS, "Something is wrong with your membrane. Max and min z values are %f and %f.\n" + "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n" + "so that one membrane is distributed over two periodic box images. Another possibility is that\n" + "your water layer is not thick enough.\n", zmax, zmin); + } + mem_p->zmin = zmin; + mem_p->zmax = zmax; + mem_p->zmed = (zmax-zmin)/2+zmin; + + /*number of membrane molecules in protein box*/ + nmolbox = count/mtop->moltype[mtop->molblock[block].type].atoms.nr; + /*membrane area within the box defined by the min and max coordinates of the embedded molecule*/ + mem_area = (pos_ins->xmax[XX]-pos_ins->xmin[XX])*(pos_ins->xmax[YY]-pos_ins->xmin[YY]); + /*rough estimate of area per lipid, assuming there is only one type of lipid in the membrane*/ + mem_p->lip_area = 2.0*mem_area/static_cast(nmolbox); + + return mem_p->mem_at.nr; +} + +static void init_resize(t_block *ins_at, rvec *r_ins, pos_ins_t *pos_ins, mem_t *mem_p, rvec *r, + gmx_bool bALLOW_ASYMMETRY) +{ + int i, j, at, c, outsidesum, gctr = 0; + int idxsum = 0; + + /*sanity check*/ + for (i = 0; i < pos_ins->pieces; i++) + { + idxsum += pos_ins->nidx[i]; + } + + if (idxsum != ins_at->nr) + { + gmx_fatal(FARGS, "Piecewise sum of inserted atoms not same as size of group selected to insert."); + } + + snew(pos_ins->geom_cent, pos_ins->pieces); + for (i = 0; i < pos_ins->pieces; i++) + { + c = 0; + outsidesum = 0; + for (j = 0; j < DIM; j++) + { + pos_ins->geom_cent[i][j] = 0; + } + + for (j = 0; j < pos_ins->nidx[i]; j++) + { + at = pos_ins->subindex[i][j]; + copy_rvec(r[at], r_ins[gctr]); + if ( (r_ins[gctr][ZZ] < mem_p->zmax) && (r_ins[gctr][ZZ] > mem_p->zmin) ) + { + rvec_inc(pos_ins->geom_cent[i], r_ins[gctr]); + c++; + } + else + { + outsidesum++; + } + gctr++; + } + + if (c > 0) + { + svmul(1/static_cast(c), pos_ins->geom_cent[i], pos_ins->geom_cent[i]); + } + + if (!bALLOW_ASYMMETRY) + { + pos_ins->geom_cent[i][ZZ] = mem_p->zmed; + } + + fprintf(stderr, "Embedding piece %d with center of geometry: %f %f %f\n", + i, pos_ins->geom_cent[i][XX], pos_ins->geom_cent[i][YY], pos_ins->geom_cent[i][ZZ]); + } + fprintf(stderr, "\n"); +} + +/* resize performed in the md loop */ +static void resize(rvec *r_ins, rvec *r, pos_ins_t *pos_ins, const rvec fac) +{ + int i, j, k, at, c = 0; + for (k = 0; k < pos_ins->pieces; k++) + { + for (i = 0; i < pos_ins->nidx[k]; i++) + { + at = pos_ins->subindex[k][i]; + for (j = 0; j < DIM; j++) + { + r[at][j] = pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]); + } + c++; + } + } +} + +/* generate the list of membrane molecules that overlap with the molecule to be embedded. * + * The molecule to be embedded is already reduced in size. */ +static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc, gmx_mtop_t *mtop, + rvec *r, mem_t *mem_p, pos_ins_t *pos_ins, real probe_rad, + int low_up_rm, gmx_bool bALLOW_ASYMMETRY) +{ + int i, j, k, l, at, at2, mol_id; + int type = 0, block = 0; + int nrm, nupper, nlower; + real r_min_rad, z_lip, min_norm; + gmx_bool bRM; + rvec dr, dr_tmp; + real *dist; + int *order; + + r_min_rad = probe_rad*probe_rad; + gmx::RangePartitioning molecules = gmx_mtop_molecules(*mtop); + snew(rm_p->block, molecules.numBlocks()); + nrm = nupper = 0; + nlower = low_up_rm; + for (i = 0; i < ins_at->nr; i++) + { + at = ins_at->index[i]; + for (j = 0; j < rest_at->nr; j++) + { + at2 = rest_at->index[j]; + pbc_dx(pbc, r[at], r[at2], dr); + + if (norm2(dr) < r_min_rad) + { + mol_id = get_mol_id(at2, mtop, &type, &block); + bRM = TRUE; + for (l = 0; l < nrm; l++) + { + if (rm_p->mol[l] == mol_id) + { + bRM = FALSE; + } + } + + if (bRM) + { + rm_p->mol[nrm] = mol_id; + rm_p->block[nrm] = block; + nrm++; + z_lip = 0.0; + for (l = 0; l < mem_p->nmol; l++) + { + if (mol_id == mem_p->mol_id[l]) + { + for (int k : molecules.block(mol_id)) + { + z_lip += r[k][ZZ]; + } + z_lip /= molecules.block(mol_id).size(); + if (z_lip < mem_p->zmed) + { + nlower++; + } + else + { + nupper++; + } + } + } + } + } + } + } + + /*make sure equal number of lipids from upper and lower layer are removed */ + if ( (nupper != nlower) && (!bALLOW_ASYMMETRY) ) + { + snew(dist, mem_p->nmol); + snew(order, mem_p->nmol); + for (i = 0; i < mem_p->nmol; i++) + { + at = molecules.block(mem_p->mol_id[i]).begin(); + pbc_dx(pbc, r[at], pos_ins->geom_cent[0], dr); + if (pos_ins->pieces > 1) + { + /*minimum dr value*/ + min_norm = norm2(dr); + for (k = 1; k < pos_ins->pieces; k++) + { + pbc_dx(pbc, r[at], pos_ins->geom_cent[k], dr_tmp); + if (norm2(dr_tmp) < min_norm) + { + min_norm = norm2(dr_tmp); + copy_rvec(dr_tmp, dr); + } + } + } + dist[i] = dr[XX]*dr[XX]+dr[YY]*dr[YY]; + j = i-1; + while (j >= 0 && dist[i] < dist[order[j]]) + { + order[j+1] = order[j]; + j--; + } + order[j+1] = i; + } + + i = 0; + while (nupper != nlower) + { + mol_id = mem_p->mol_id[order[i]]; + block = get_molblock(mol_id, mtop->molblock); + bRM = TRUE; + for (l = 0; l < nrm; l++) + { + if (rm_p->mol[l] == mol_id) + { + bRM = FALSE; + } + } + + if (bRM) + { + z_lip = 0; + for (int k : molecules.block(mol_id)) + { + z_lip += r[k][ZZ]; + } + z_lip /= molecules.block(mol_id).size(); + if (nupper > nlower && z_lip < mem_p->zmed) + { + rm_p->mol[nrm] = mol_id; + rm_p->block[nrm] = block; + nrm++; + nlower++; + } + else if (nupper < nlower && z_lip > mem_p->zmed) + { + rm_p->mol[nrm] = mol_id; + rm_p->block[nrm] = block; + nrm++; + nupper++; + } + } + i++; + + if (i > mem_p->nmol) + { + gmx_fatal(FARGS, "Trying to remove more lipid molecules than there are in the membrane"); + } + } + sfree(dist); + sfree(order); + } + + rm_p->nr = nrm; + srenew(rm_p->mol, nrm); + srenew(rm_p->block, nrm); + + return nupper+nlower; +} + +/*remove all lipids and waters overlapping and update all important structures (e.g. state and mtop)*/ +static void rm_group(gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state *state, + t_block *ins_at, pos_ins_t *pos_ins) +{ + int j, k, n, rm, mol_id, at, block; + rvec *x_tmp, *v_tmp; + int *list; + unsigned char *new_egrp[egcNR]; + gmx_bool bRM; + int RMmolblock; + + /* Construct the molecule range information */ + gmx::RangePartitioning molecules = gmx_mtop_molecules(*mtop); + + snew(list, state->natoms); + n = 0; + for (int i = 0; i < rm_p->nr; i++) + { + mol_id = rm_p->mol[i]; + at = molecules.block(mol_id).size(); + block = rm_p->block[i]; + mtop->molblock[block].nmol--; + for (j = 0; j < mtop->moltype[mtop->molblock[block].type].atoms.nr; j++) + { + list[n] = at+j; + n++; + } + } + + mtop->natoms -= n; - state->natoms -= n; ++ state_change_natoms(state, state->natoms - n); + snew(x_tmp, state->natoms); + snew(v_tmp, state->natoms); + + for (int i = 0; i < egcNR; i++) + { + if (groups->grpnr[i] != nullptr) + { + groups->ngrpnr[i] = state->natoms; + snew(new_egrp[i], state->natoms); + } + } + + auto x = makeArrayRef(state->x); + auto v = makeArrayRef(state->v); + rm = 0; + for (int i = 0; i < state->natoms+n; i++) + { + bRM = FALSE; + for (j = 0; j < n; j++) + { + if (i == list[j]) + { + bRM = TRUE; + rm++; + } + } + + if (!bRM) + { + for (j = 0; j < egcNR; j++) + { + if (groups->grpnr[j] != nullptr) + { + new_egrp[j][i-rm] = groups->grpnr[j][i]; + } + } + copy_rvec(x[i], x_tmp[i-rm]); + copy_rvec(v[i], v_tmp[i-rm]); + for (j = 0; j < ins_at->nr; j++) + { + if (i == ins_at->index[j]) + { + ins_at->index[j] = i-rm; + } + } + + for (j = 0; j < pos_ins->pieces; j++) + { + for (k = 0; k < pos_ins->nidx[j]; k++) + { + if (i == pos_ins->subindex[j][k]) + { + pos_ins->subindex[j][k] = i-rm; + } + } + } + } + } + for (int i = 0; i < state->natoms; i++) + { + copy_rvec(x_tmp[i], x[i]); + } + sfree(x_tmp); + for (int i = 0; i < state->natoms; i++) + { + copy_rvec(v_tmp[i], v[i]); + } + sfree(v_tmp); + + for (int i = 0; i < egcNR; i++) + { + if (groups->grpnr[i] != nullptr) + { + sfree(groups->grpnr[i]); + groups->grpnr[i] = new_egrp[i]; + } + } + + /* remove empty molblocks */ + RMmolblock = 0; + for (size_t i = 0; i < mtop->molblock.size(); i++) + { + if (mtop->molblock[i].nmol == 0) + { + RMmolblock++; + } + else + { + mtop->molblock[i-RMmolblock] = mtop->molblock[i]; + } + } + mtop->molblock.resize(mtop->molblock.size() - RMmolblock); +} + +/* remove al bonded interactions from mtop for the molecule to be embedded */ +static int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop) +{ + int j, m; + int type, natom, nmol, at, atom1 = 0, rm_at = 0; + gmx_bool *bRM, bINS; + /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/ + /*this routine does not live as dangerously as it seems. There is namely a check in init_membed to make * + * sure that g_membed exits with a warning when there are molecules of the same type not in the * + * ins_at index group. MGWolf 050710 */ + + + snew(bRM, mtop->moltype.size()); + for (size_t i = 0; i < mtop->moltype.size(); i++) + { + bRM[i] = TRUE; + } + + for (size_t i = 0; i < mtop->molblock.size(); i++) + { + /*loop over molecule blocks*/ + type = mtop->molblock[i].type; + natom = mtop->moltype[type].atoms.nr; + nmol = mtop->molblock[i].nmol; + + for (j = 0; j < natom*nmol && bRM[type]; j++) + { + /*loop over atoms in the block*/ + at = j+atom1; /*atom index = block index + offset*/ + bINS = FALSE; + + for (m = 0; (m < ins_at->nr) && (!bINS); m++) + { + /*loop over atoms in insertion index group to determine if we're inserting one*/ + if (at == ins_at->index[m]) + { + bINS = TRUE; + } + } + bRM[type] = bINS; + } + atom1 += natom*nmol; /*update offset*/ + if (bRM[type]) + { + rm_at += natom*nmol; /*increment bonded removal counter by # atoms in block*/ + } + } + + for (size_t i = 0; i < mtop->moltype.size(); i++) + { + if (bRM[i]) + { + for (j = 0; j < F_LJ; j++) + { + mtop->moltype[i].ilist[j].iatoms.clear(); + } + + for (j = F_POSRES; j <= F_VSITEN; j++) + { + mtop->moltype[i].ilist[j].iatoms.clear(); + } + } + } + sfree(bRM); + + return rm_at; +} + +/* Write a topology where the number of molecules is correct for the system after embedding */ +static void top_update(const char *topfile, rm_t *rm_p, gmx_mtop_t *mtop) +{ + int bMolecules = 0; + FILE *fpin, *fpout; + char buf[STRLEN], buf2[STRLEN], *temp; + int i, *nmol_rm, nmol, line; + char temporary_filename[STRLEN]; + + fpin = gmx_ffopen(topfile, "r"); + strncpy(temporary_filename, "temp.topXXXXXX", STRLEN); + gmx_tmpnam(temporary_filename); + fpout = gmx_ffopen(temporary_filename, "w"); + + snew(nmol_rm, mtop->moltype.size()); + for (i = 0; i < rm_p->nr; i++) + { + nmol_rm[rm_p->block[i]]++; + } + + line = 0; + while (fgets(buf, STRLEN, fpin)) + { + line++; + if (buf[0] != ';') + { + strcpy(buf2, buf); + if ((temp = strchr(buf2, '\n')) != nullptr) + { + temp[0] = '\0'; + } + ltrim(buf2); + if (buf2[0] == '[') + { + buf2[0] = ' '; + if ((temp = strchr(buf2, '\n')) != nullptr) + { + temp[0] = '\0'; + } + rtrim(buf2); + if (buf2[strlen(buf2)-1] == ']') + { + buf2[strlen(buf2)-1] = '\0'; + ltrim(buf2); + rtrim(buf2); + if (gmx_strcasecmp(buf2, "molecules") == 0) + { + bMolecules = 1; + } + } + fprintf(fpout, "%s", buf); + } + else if (bMolecules == 1) + { + for (size_t i = 0; i < mtop->molblock.size(); i++) + { + nmol = mtop->molblock[i].nmol; + sprintf(buf, "%-15s %5d\n", *(mtop->moltype[mtop->molblock[i].type].name), nmol); + fprintf(fpout, "%s", buf); + } + bMolecules = 2; + } + else if (bMolecules == 2) + { + /* print nothing */ + } + else + { + fprintf(fpout, "%s", buf); + } + } + else + { + fprintf(fpout, "%s", buf); + } + } + + gmx_ffclose(fpout); + /* use gmx_ffopen to generate backup of topinout */ + fpout = gmx_ffopen(topfile, "w"); + gmx_ffclose(fpout); + rename(temporary_filename, topfile); +} + +void rescale_membed(int step_rel, gmx_membed_t *membed, rvec *x) +{ + /* Set new positions for the group to embed */ + if (step_rel <= membed->it_xy) + { + membed->fac[0] += membed->xy_step; + membed->fac[1] += membed->xy_step; + } + else if (step_rel <= (membed->it_xy+membed->it_z)) + { + membed->fac[2] += membed->z_step; + } + resize(membed->r_ins, x, membed->pos_ins, membed->fac); +} + +/* We would like gn to be const as well, but C doesn't allow this */ +/* TODO this is utility functionality (search for the index of a + string in a collection), so should be refactored and located more + centrally. Also, it nearly duplicates the same string in readir.c */ +static int search_string(const char *s, int ng, char *gn[]) +{ + int i; + + for (i = 0; (i < ng); i++) + { + if (gmx_strcasecmp(s, gn[i]) == 0) + { + return i; + } + } + + gmx_fatal(FARGS, + "Group %s selected for embedding was not found in the index file.\n" + "Group names must match either [moleculetype] names or custom index group\n" + "names, in which case you must supply an index file to the '-n' option\n" + "of grompp.", + s); +} + +gmx_membed_t *init_membed(FILE *fplog, int nfile, const t_filenm fnm[], gmx_mtop_t *mtop, + t_inputrec *inputrec, t_state *state, t_commrec *cr, real *cpt) +{ + char *ins, **gnames; + int i, rm_bonded_at, fr_id, fr_i = 0, tmp_id, warn = 0; + int ng, j, max_lip_rm, ins_grp_id, ntype, lip_rm; + real prot_area; + rvec *r_ins = nullptr; + t_block *ins_at, *rest_at; + pos_ins_t *pos_ins; + mem_t *mem_p; + rm_t *rm_p; + gmx_groups_t *groups; + gmx_bool bExcl = FALSE; + t_atoms atoms; + t_pbc *pbc; + char **piecename = nullptr; + gmx_membed_t *membed = nullptr; + + /* input variables */ + real xy_fac = 0.5; + real xy_max = 1.0; + real z_fac = 1.0; + real z_max = 1.0; + int it_xy = 1000; + int it_z = 0; + real probe_rad = 0.22; + int low_up_rm = 0; + int maxwarn = 0; + int pieces = 1; + gmx_bool bALLOW_ASYMMETRY = FALSE; + + /* sanity check constants */ /* Issue a warning when: */ + const real min_probe_rad = 0.2199999; /* A probe radius for overlap between embedded molecule * + * and rest smaller than this value is probably too small */ + const real min_xy_init = 0.0999999; /* the initial shrinking of the molecule to embed is smaller */ + const int min_it_xy = 1000; /* the number of steps to embed in xy-plane is smaller */ + const int min_it_z = 100; /* the number of steps to embed in z is smaller */ + const real prot_vs_box = 7.5; /* molecule to embed is large (more then prot_vs_box) with respect */ + const real box_vs_prot = 50; /* to the box size (less than box_vs_prot) */ + + snew(membed, 1); + snew(ins_at, 1); + snew(pos_ins, 1); + + if (MASTER(cr)) + { + fprintf(fplog, "Note: it is expected that in future gmx mdrun -membed will not be the " + "way to access this feature, perhaps changing to e.g. gmx membed."); + /* get input data out membed file */ + try + { + get_input(opt2fn("-membed", nfile, fnm), + &xy_fac, &xy_max, &z_fac, &z_max, &it_xy, &it_z, &probe_rad, &low_up_rm, + &maxwarn, &pieces, &bALLOW_ASYMMETRY); + } + GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR; + + if (!EI_DYNAMICS(inputrec->eI) ) + { + gmx_input("Change integrator to a dynamics integrator in mdp file (e.g. md or sd)."); + } + + if (PAR(cr)) + { + gmx_input("Sorry, parallel membed is not yet fully functional."); + } + + if (*cpt >= 0) + { + fprintf(stderr, "\nSetting -cpt to -1, because embedding cannot be restarted from cpt-files.\n"); + *cpt = -1; + } + groups = &(mtop->groups); + snew(gnames, groups->ngrpname); + for (i = 0; i < groups->ngrpname; i++) + { + gnames[i] = *(groups->grpname[i]); + } + + atoms = gmx_mtop_global_atoms(mtop); + snew(mem_p, 1); + fprintf(stderr, "\nSelect a group to embed in the membrane:\n"); + get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(ins_at->nr), &(ins_at->index), &ins); + ins_grp_id = search_string(ins, groups->ngrpname, gnames); + fprintf(stderr, "\nSelect a group to embed %s into (e.g. the membrane):\n", ins); + get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(mem_p->mem_at.nr), &(mem_p->mem_at.index), &(mem_p->name)); + + pos_ins->pieces = pieces; + snew(pos_ins->nidx, pieces); + snew(pos_ins->subindex, pieces); + snew(piecename, pieces); + if (pieces > 1) + { + fprintf(stderr, "\nSelect pieces to embed:\n"); + get_index(&atoms, opt2fn_null("-mn", nfile, fnm), pieces, pos_ins->nidx, pos_ins->subindex, piecename); + } + else + { + /*use whole embedded group*/ + snew(pos_ins->nidx, 1); + snew(pos_ins->subindex, 1); + pos_ins->nidx[0] = ins_at->nr; + pos_ins->subindex[0] = ins_at->index; + } + + if (probe_rad < min_probe_rad) + { + warn++; + fprintf(stderr, "\nWarning %d:\nA probe radius (-rad) smaller than 0.2 nm can result " + "in overlap between waters and the group to embed, which will result " + "in Lincs errors etc.\n\n", warn); + } + + if (xy_fac < min_xy_init) + { + warn++; + fprintf(stderr, "\nWarning %d:\nThe initial size of %s is probably too small.\n\n", warn, ins); + } + + if (it_xy < min_it_xy) + { + warn++; + fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d)" + " is probably too small.\nIncrease -nxy or.\n\n", warn, ins, it_xy); + } + + if ( (it_z < min_it_z) && ( z_fac < 0.99999999 || z_fac > 1.0000001) ) + { + warn++; + fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d)" + " is probably too small.\nIncrease -nz or the maxwarn setting in the membed input file.\n\n", warn, ins, it_z); + } + + if (it_xy+it_z > inputrec->nsteps) + { + warn++; + fprintf(stderr, "\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the " + "number of steps in the tpr.\n" + "(increase maxwarn in the membed input file to override)\n\n", warn); + } + + fr_id = -1; + if (inputrec->opts.ngfrz == 1) + { + gmx_fatal(FARGS, "You did not specify \"%s\" as a freezegroup.", ins); + } + + for (i = 0; i < inputrec->opts.ngfrz; i++) + { + tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i]; + if (ins_grp_id == tmp_id) + { + fr_id = tmp_id; + fr_i = i; + } + } + + if (fr_id == -1) + { + gmx_fatal(FARGS, "\"%s\" not as freezegroup defined in the mdp-file.", ins); + } + + for (i = 0; i < DIM; i++) + { + if (inputrec->opts.nFreeze[fr_i][i] != 1) + { + gmx_fatal(FARGS, "freeze dimensions for %s are not Y Y Y\n", ins); + } + } + + ng = groups->grps[egcENER].nr; + if (ng == 1) + { + gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group"); + } + + for (i = 0; i < ng; i++) + { + for (j = 0; j < ng; j++) + { + if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL) + { + bExcl = TRUE; + if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) || + (groups->grps[egcENER].nm_ind[j] != ins_grp_id) ) + { + gmx_fatal(FARGS, "Energy exclusions \"%s\" and \"%s\" do not match the group " + "to embed \"%s\"", *groups->grpname[groups->grps[egcENER].nm_ind[i]], + *groups->grpname[groups->grps[egcENER].nm_ind[j]], ins); + } + } + } + } + + if (!bExcl) + { + gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in " + "the freeze group"); + } + + /* Obtain the maximum and minimum coordinates of the group to be embedded */ + snew(rest_at, 1); + init_ins_at(ins_at, rest_at, state, pos_ins, groups, ins_grp_id, xy_max); + /* Check that moleculetypes in insertion group are not part of the rest of the system */ + check_types(ins_at, rest_at, mtop); + + init_mem_at(mem_p, mtop, state->x.rvec_array(), state->box, pos_ins); + + prot_area = est_prot_area(pos_ins, state->x.rvec_array(), ins_at, mem_p); + if ( (prot_area > prot_vs_box) && ( (state->box[XX][XX]*state->box[YY][YY]-state->box[XX][YY]*state->box[YY][XX]) < box_vs_prot) ) + { + warn++; + fprintf(stderr, "\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n" + "This might cause pressure problems during the growth phase. Just try with\n" + "current setup and increase 'maxwarn' in your membed settings file, but lower the\n" + "compressibility in the mdp-file or disable pressure coupling if problems occur.\n\n", warn); + } + + if (warn > maxwarn) + { + gmx_fatal(FARGS, "Too many warnings (override by setting maxwarn in the membed input file)\n"); + } + + printf("The estimated area of the protein in the membrane is %.3f nm^2\n", prot_area); + printf("\nThere are %d lipids in the membrane part that overlaps the protein.\n" + "The area per lipid is %.4f nm^2.\n", mem_p->nmol, mem_p->lip_area); + + /* Maximum number of lipids to be removed*/ + max_lip_rm = static_cast(2*prot_area/mem_p->lip_area); + printf("Maximum number of lipids that will be removed is %d.\n", max_lip_rm); + + printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n" + "This resizing will be done with respect to the geometrical center of all protein atoms\n" + "that span the membrane region, i.e. z between %.3f and %.3f\n\n", + xy_fac, z_fac, mem_p->zmin, mem_p->zmax); + + /* resize the protein by xy and by z if necessary*/ + snew(r_ins, ins_at->nr); + init_resize(ins_at, r_ins, pos_ins, mem_p, state->x.rvec_array(), bALLOW_ASYMMETRY); + membed->fac[0] = membed->fac[1] = xy_fac; + membed->fac[2] = z_fac; + + membed->xy_step = (xy_max-xy_fac)/static_cast(it_xy); + membed->z_step = (z_max-z_fac)/static_cast(it_z-1); + + resize(r_ins, state->x.rvec_array(), pos_ins, membed->fac); + + /* remove overlapping lipids and water from the membrane box*/ + /*mark molecules to be removed*/ + snew(pbc, 1); + set_pbc(pbc, inputrec->ePBC, state->box); + + snew(rm_p, 1); + lip_rm = gen_rm_list(rm_p, ins_at, rest_at, pbc, mtop, state->x.rvec_array(), mem_p, pos_ins, + probe_rad, low_up_rm, bALLOW_ASYMMETRY); + lip_rm -= low_up_rm; + + if (fplog) + { + for (i = 0; i < rm_p->nr; i++) + { + fprintf(fplog, "rm mol %d\n", rm_p->mol[i]); + } + } + + for (size_t i = 0; i < mtop->molblock.size(); i++) + { + ntype = 0; + for (j = 0; j < rm_p->nr; j++) + { + if (rm_p->block[j] == static_cast(i)) + { + ntype++; + } + } + printf("Will remove %d %s molecules\n", ntype, *(mtop->moltype[mtop->molblock[i].type].name)); + } + + if (lip_rm > max_lip_rm) + { + warn++; + fprintf(stderr, "\nWarning %d:\nTrying to remove a larger lipid area than the estimated " + "protein area\nTry making the -xyinit resize factor smaller or increase " + "maxwarn in the membed input file.\n\n", warn); + } + + /*remove all lipids and waters overlapping and update all important structures*/ + rm_group(groups, mtop, rm_p, state, ins_at, pos_ins); + + rm_bonded_at = rm_bonded(ins_at, mtop); + if (rm_bonded_at != ins_at->nr) + { + fprintf(stderr, "Warning: The number of atoms for which the bonded interactions are removed is %d, " + "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same" + "molecule type as atoms that are not to be embedded.\n", rm_bonded_at, ins_at->nr); + } + + if (warn > maxwarn) + { + gmx_fatal(FARGS, "Too many warnings.\nIf you are sure these warnings are harmless,\n" + "you can increase the maxwarn setting in the membed input file."); + } + + // Re-establish the invariants of the derived values within + // mtop. + gmx_mtop_finalize(mtop); + + if (ftp2bSet(efTOP, nfile, fnm)) + { + top_update(opt2fn("-mp", nfile, fnm), rm_p, mtop); + } + + sfree(pbc); + sfree(rest_at); + if (pieces > 1) + { + sfree(piecename); + } + + membed->it_xy = it_xy; + membed->it_z = it_z; + membed->pos_ins = pos_ins; + membed->r_ins = r_ins; + } + + return membed; +} + +void free_membed(gmx_membed_t *membed) +{ + sfree(membed); +} diff --cc src/gromacs/mdtypes/enerdata.h index eccbe71c94,0000000000..2087c8b44a mode 100644,000000..100644 --- a/src/gromacs/mdtypes/enerdata.h +++ b/src/gromacs/mdtypes/enerdata.h @@@ -1,66 -1,0 +1,71 @@@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2018, 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 GMX_MDTYPES_TYPES_ENERDATA_H +#define GMX_MDTYPES_TYPES_ENERDATA_H + +#include "gromacs/mdtypes/md_enums.h" +#include "gromacs/topology/idef.h" +#include "gromacs/utility/real.h" + +enum { + egCOULSR, egLJSR, egBHAMSR, + egCOUL14, egLJ14, egNR +}; + +struct gmx_grppairener_t +{ + int nener; /* The number of energy group pairs */ + real *ener[egNR]; /* Energy terms for each pair of groups */ +}; + +struct gmx_enerdata_t +{ + real term[F_NRE]; /* The energies for all different interaction types */ + struct gmx_grppairener_t grpp; + double dvdl_lin[efptNR]; /* Contributions to dvdl with linear lam-dependence */ + double dvdl_nonlin[efptNR]; /* Idem, but non-linear dependence */ ++ /* The idea is that dvdl terms with linear lambda dependence will be added ++ * automatically to enerpart_lambda. Terms with non-linear lambda dependence ++ * should explicitly determine the energies at foreign lambda points ++ * when n_lambda > 0. */ ++ + int n_lambda; + int fep_state; /*current fep state -- just for printing */ - double *enerpart_lambda; /* Partial energy for lambda and flambda[] */ ++ double *enerpart_lambda; /* Partial Hamiltonian for lambda and flambda[], includes at least all perturbed terms */ + real foreign_term[F_NRE]; /* alternate array for storing foreign lambda energies */ + struct gmx_grppairener_t foreign_grpp; /* alternate array for storing foreign lambda energies */ +}; + +#endif