2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
44 #include "gromacs/legacyheaders/chargegroup.h"
45 #include "gromacs/legacyheaders/force.h"
46 #include "gromacs/legacyheaders/macros.h"
47 #include "gromacs/legacyheaders/mdatoms.h"
48 #include "gromacs/legacyheaders/names.h"
49 #include "gromacs/legacyheaders/nrnb.h"
50 #include "gromacs/legacyheaders/ns.h"
51 #include "gromacs/legacyheaders/nsgrid.h"
52 #include "gromacs/legacyheaders/types/commrec.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/pbcutil/pbc.h"
55 #include "gromacs/topology/block.h"
56 #include "gromacs/topology/mtop_util.h"
57 #include "gromacs/topology/topology.h"
58 #include "gromacs/utility/smalloc.h"
60 static real box_margin;
62 static real max_dist(rvec *x, real *r, int start, int end)
68 for (i = start; i < end; i++)
70 for (j = i+1; j < end; j++)
72 maxd = max(maxd, sqrt(distance2(x[i], x[j]))+0.5*(r[i]+r[j]));
79 static gmx_bool outside_box_minus_margin2(rvec x, matrix box)
81 return ( (x[XX] < 2*box_margin) || (x[XX] > box[XX][XX]-2*box_margin) ||
82 (x[YY] < 2*box_margin) || (x[YY] > box[YY][YY]-2*box_margin) ||
83 (x[ZZ] < 2*box_margin) || (x[ZZ] > box[ZZ][ZZ]-2*box_margin) );
86 static gmx_bool outside_box_plus_margin(rvec x, matrix box)
88 return ( (x[XX] < -box_margin) || (x[XX] > box[XX][XX]+box_margin) ||
89 (x[YY] < -box_margin) || (x[YY] > box[YY][YY]+box_margin) ||
90 (x[ZZ] < -box_margin) || (x[ZZ] > box[ZZ][ZZ]+box_margin) );
93 static int mark_res(int at, gmx_bool *mark, int natoms, t_atom *atom, int *nmark)
97 resind = atom[at].resind;
98 while ( (at > 0) && (resind == atom[at-1].resind) )
102 while ( (at < natoms) && (resind == atom[at].resind) )
115 static real find_max_real(int n, real radius[])
124 for (i = 1; (i < n); i++)
126 rmax = max(rmax, radius[i]);
132 static void combine_atoms(t_atoms *ap, t_atoms *as,
133 rvec xp[], rvec *vp, rvec xs[], rvec *vs,
134 t_atoms **a_comb, rvec **x_comb, rvec **v_comb)
137 rvec *xc, *vc = NULL;
138 int i, j, natot, res0;
140 /* Total number of atoms */
141 natot = ap->nr+as->nr;
144 init_t_atoms(ac, natot, FALSE);
152 /* Fill the new structures */
153 for (i = j = 0; (i < ap->nr); i++, j++)
155 copy_rvec(xp[i], xc[j]);
158 copy_rvec(vp[i], vc[j]);
160 memcpy(&(ac->atom[j]), &(ap->atom[i]), sizeof(ap->atom[i]));
161 ac->atom[j].type = 0;
164 for (i = 0; (i < as->nr); i++, j++)
166 copy_rvec(xs[i], xc[j]);
169 copy_rvec(vs[i], vc[j]);
171 memcpy(&(ac->atom[j]), &(as->atom[i]), sizeof(as->atom[i]));
172 ac->atom[j].type = 0;
173 ac->atom[j].resind += res0;
176 ac->nres = ac->atom[j-1].resind+1;
177 /* Fill all elements to prevent uninitialized memory */
178 for (i = 0; i < ac->nr; i++)
184 ac->atom[i].type = 0;
185 ac->atom[i].typeB = 0;
186 ac->atom[i].ptype = eptAtom;
195 static t_forcerec *fr = NULL;
197 static void do_nsgrid(FILE *fp, gmx_bool bVerbose,
198 matrix box, rvec x[], t_atoms *atoms, real rlong,
199 const output_env_t oenv)
218 /* Charge group index */
219 snew(cg_index, natoms);
220 for (i = 0; (i < natoms); i++)
225 /* Topology needs charge groups and exclusions */
228 mtop->natoms = natoms;
229 /* Make one moltype that contains the whol system */
231 snew(mtop->moltype, mtop->nmoltype);
232 molt = &mtop->moltype[0];
233 molt->name = mtop->name;
234 molt->atoms = *atoms;
235 stupid_fill_block(&molt->cgs, mtop->natoms, FALSE);
236 stupid_fill_blocka(&molt->excls, natoms);
237 /* Make one molblock for the whole system */
239 snew(mtop->molblock, mtop->nmolblock);
240 mtop->molblock[0].type = 0;
241 mtop->molblock[0].nmol = 1;
242 mtop->molblock[0].natoms_mol = natoms;
243 /* Initialize a single energy group */
244 mtop->groups.grps[egcENER].nr = 1;
245 mtop->groups.ngrpnr[egcENER] = 0;
246 mtop->groups.grpnr[egcENER] = NULL;
248 ffp = &mtop->ffparams;
253 snew(ffp->functype, 1);
254 snew(ffp->iparams, 1);
255 ffp->iparams[0].lj.c6 = 1;
256 ffp->iparams[0].lj.c12 = 1;
258 /* inputrec structure */
260 ir->cutoff_scheme = ecutsGROUP;
261 ir->coulomb_modifier = eintmodNONE;
262 ir->vdw_modifier = eintmodNONE;
263 ir->coulombtype = eelCUT;
264 ir->vdwtype = evdwCUT;
265 ir->ns_type = ensGRID;
266 snew(ir->opts.egp_flags, 1);
268 top = gmx_mtop_generate_local_top(mtop, ir);
270 /* Some nasty shortcuts */
273 /* mdatoms structure */
276 md = init_mdatoms(fp, mtop, FALSE);
277 atoms2md(mtop, ir, 0, NULL, mtop->natoms, md);
280 /* forcerec structure */
287 /* cr->nthreads = 1; */
289 /* ir->rlist = ir->rcoulomb = ir->rvdw = rlong;
290 printf("Neighborsearching with a cut-off of %g\n",rlong);
291 init_forcerec(stdout,fr,ir,top,cr,md,box,FALSE,NULL,NULL,NULL,TRUE);*/
293 fr->hcg = top->cgs.nr;
296 /* Prepare for neighboursearching */
299 /* Init things dependent on parameters */
300 ir->rlistlong = ir->rlist = ir->rcoulomb = ir->rvdw = rlong;
301 /* create free energy data to avoid NULLs */
302 snew(ir->fepvals, 1);
303 printf("Neighborsearching with a cut-off of %g\n", rlong);
304 init_forcerec(stdout, oenv, fr, NULL, ir, mtop, cr, box,
305 NULL, NULL, NULL, NULL, NULL, TRUE, -1);
308 pr_forcerec(debug, fr);
311 /* Calculate new stuff dependent on coords and box */
312 for (m = 0; (m < DIM); m++)
314 box_size[m] = box[m][m];
316 calc_shifts(box, fr->shift_vec);
317 put_charge_groups_in_box(fp, 0, cgs->nr, fr->ePBC, box, cgs, x, fr->cg_cm);
319 /* Do the actual neighboursearching */
320 snew(lambda, efptNR);
322 init_neighbor_list(fp, fr, md->homenr);
323 search_neighbours(fp, fr, box, top,
324 &mtop->groups, cr, &nrnb, md, TRUE, FALSE);
328 dump_nblist(debug, cr, fr, 0);
333 fprintf(stderr, "Successfully made neighbourlist\n");
337 static gmx_bool bXor(gmx_bool b1, gmx_bool b2)
339 return (b1 && !b2) || (b2 && !b1);
342 void add_conf(t_atoms *atoms, rvec **x, rvec **v, real **r, gmx_bool bSrenew,
343 int ePBC, matrix box, gmx_bool bInsert,
344 t_atoms *atoms_solvt, rvec *x_solvt, rvec *v_solvt, real *r_solvt,
345 gmx_bool bVerbose, real rshell, int max_sol, const output_env_t oenv)
349 real max_vdw, *r_prot, *r_all, n2, r2, ib1, ib2;
350 int natoms_prot, natoms_solvt;
351 int i, j, jj, m, j0, j1, jjj, jnres, jnr, inr, iprot, is1, is2;
352 int prev, resnr, nresadd, d, k, ncells, maxincell;
353 int dx0, dx1, dy0, dy1, dz0, dz1;
354 int ntest, nremove, nkeep;
355 rvec dx, xi, xj, xpp, *x_all, *v_all;
356 gmx_bool *remove, *keep;
359 natoms_prot = atoms->nr;
360 natoms_solvt = atoms_solvt->nr;
361 if (natoms_solvt <= 0)
363 fprintf(stderr, "WARNING: Nothing to add\n");
367 if (ePBC == epbcSCREW)
369 gmx_fatal(FARGS, "Sorry, %s pbc is not yet supported", epbc_names[ePBC]);
374 fprintf(stderr, "Calculating Overlap...\n");
377 /* Set margin around box edges to largest solvent dimension.
378 * The maximum distance between atoms in a solvent molecule should
379 * be calculated. At the moment a fudge factor of 3 is used.
382 box_margin = 3*find_max_real(natoms_solvt, r_solvt);
383 max_vdw = max(3*find_max_real(natoms_prot, r_prot), box_margin);
384 fprintf(stderr, "box_margin = %g\n", box_margin);
386 snew(remove, natoms_solvt);
391 for (i = 0; i < atoms_solvt->nr; i++)
393 if (outside_box_plus_margin(x_solvt[i], box) )
395 i = mark_res(i, remove, atoms_solvt->nr, atoms_solvt->atom, &nremove);
398 fprintf(stderr, "Removed %d atoms that were outside the box\n", nremove);
401 /* Define grid stuff */
402 /* Largest VDW radius */
403 snew(r_all, natoms_prot+natoms_solvt);
404 for (i = j = 0; i < natoms_prot; i++, j++)
406 r_all[j] = r_prot[i];
408 for (i = 0; i < natoms_solvt; i++, j++)
410 r_all[j] = r_solvt[i];
414 combine_atoms(atoms, atoms_solvt, *x, v ? *v : NULL, x_solvt, v_solvt,
415 &atoms_all, &x_all, &v_all);
417 /* Do neighboursearching step */
418 do_nsgrid(stdout, bVerbose, box, x_all, atoms_all, max_vdw, oenv);
420 /* check solvent with solute */
421 nlist = &(fr->nblists[0].nlist_sr[eNL_VDW]);
422 fprintf(stderr, "nri = %d, nrj = %d\n", nlist->nri, nlist->nrj);
423 for (bSolSol = 0; (bSolSol <= (bInsert ? 0 : 1)); bSolSol++)
426 fprintf(stderr, "Checking %s-Solvent overlap:",
427 bSolSol ? "Solvent" : "Protein");
428 for (i = 0; (i < nlist->nri && nremove < natoms_solvt); i++)
430 inr = nlist->iinr[i];
431 j0 = nlist->jindex[i];
432 j1 = nlist->jindex[i+1];
433 rvec_add(x_all[inr], fr->shift_vec[nlist->shift[i]], xi);
435 for (j = j0; (j < j1 && nremove < natoms_solvt); j++)
437 jnr = nlist->jjnr[j];
443 copy_rvec(x_all[jnr], xj);
445 /* Check solvent-protein and solvent-solvent */
446 is1 = inr-natoms_prot;
447 is2 = jnr-natoms_prot;
449 /* Check if at least one of the atoms is a solvent that is not yet
450 * listed for removal, and if both are solvent, that they are not in the
454 bXor((is1 >= 0), (is2 >= 0)) && /* One atom is protein */
455 ((is1 < 0) || ((is1 >= 0) && !remove[is1])) &&
456 ((is2 < 0) || ((is2 >= 0) && !remove[is2]))) ||
459 (is1 >= 0) && (!remove[is1]) && /* is1 is solvent */
460 (is2 >= 0) && (!remove[is2]) && /* is2 is solvent */
461 (bInsert || /* when inserting also check inside the box */
462 (outside_box_minus_margin2(x_solvt[is1], box) && /* is1 on edge */
463 outside_box_minus_margin2(x_solvt[is2], box)) /* is2 on edge */
465 (atoms_solvt->atom[is1].resind != /* Not the same residue */
466 atoms_solvt->atom[is2].resind)))
470 rvec_sub(xi, xj, dx);
472 r2 = sqr(r_all[inr]+r_all[jnr]);
477 nremove = natoms_solvt;
478 for (k = 0; k < nremove; k++)
483 /* Need only remove one of the solvents... */
486 (void) mark_res(is2, remove, natoms_solvt, atoms_solvt->atom,
491 (void) mark_res(is1, remove, natoms_solvt, atoms_solvt->atom,
496 fprintf(stderr, "Neither atom is solvent%d %d\n", is1, is2);
504 fprintf(stderr, " tested %d pairs, removed %d atoms.\n", ntest, nremove);
509 for (i = 0; i < natoms_solvt; i++)
511 fprintf(debug, "remove[%5d] = %s\n", i, bool_names[remove[i]]);
515 /* Search again, now with another cut-off */
518 do_nsgrid(stdout, bVerbose, box, x_all, atoms_all, rshell, oenv);
519 nlist = &(fr->nblists[0].nlist_sr[eNL_VDW]);
520 fprintf(stderr, "nri = %d, nrj = %d\n", nlist->nri, nlist->nrj);
522 snew(keep, natoms_solvt);
523 for (i = 0; i < nlist->nri; i++)
525 inr = nlist->iinr[i];
526 j0 = nlist->jindex[i];
527 j1 = nlist->jindex[i+1];
529 for (j = j0; j < j1; j++)
531 jnr = nlist->jjnr[j];
533 /* Check solvent-protein and solvent-solvent */
534 is1 = inr-natoms_prot;
535 is2 = jnr-natoms_prot;
537 /* Check if at least one of the atoms is a solvent that is not yet
538 * listed for removal, and if both are solvent, that they are not in the
541 if (is1 >= 0 && is2 < 0)
543 mark_res(is1, keep, natoms_solvt, atoms_solvt->atom, &nkeep);
545 else if (is1 < 0 && is2 >= 0)
547 mark_res(is2, keep, natoms_solvt, atoms_solvt->atom, &nkeep);
551 fprintf(stderr, "Keeping %d solvent atoms after proximity check\n",
553 for (i = 0; i < natoms_solvt; i++)
555 remove[i] = remove[i] || !keep[i];
559 /* count how many atoms and residues will be added and make space */
563 jnres = atoms_solvt->nres;
569 for (i = 0; ((i < atoms_solvt->nr) &&
570 ((max_sol == 0) || (jnres < max_sol))); i++)
576 (atoms_solvt->atom[i].resind != atoms_solvt->atom[i-1].resind))
585 fprintf(debug, "Will add %d atoms in %d residues\n", j, jnres);
589 /* Flag the remaing solvent atoms to be removed */
590 jjj = atoms_solvt->atom[i-1].resind;
591 for (; (i < atoms_solvt->nr); i++)
593 if (atoms_solvt->atom[i].resind > jjj)
606 srenew(atoms->resinfo, atoms->nres+jnres);
607 srenew(atoms->atomname, atoms->nr+j);
608 srenew(atoms->atom, atoms->nr+j);
609 srenew(*x, atoms->nr+j);
612 srenew(*v, atoms->nr+j);
614 srenew(*r, atoms->nr+j);
617 /* add the selected atoms_solvt to atoms */
620 resnr = atoms->resinfo[atoms->atom[atoms->nr-1].resind].nr;
628 for (i = 0; i < atoms_solvt->nr; i++)
633 atoms_solvt->atom[i].resind != atoms_solvt->atom[prev].resind)
638 atoms->resinfo[atoms->nres-1] =
639 atoms_solvt->resinfo[atoms_solvt->atom[i].resind];
640 atoms->resinfo[atoms->nres-1].nr = resnr;
641 /* calculate shift of the solvent molecule using the first atom */
642 copy_rvec(x_solvt[i], dx);
643 put_atoms_in_box(ePBC, box, 1, &dx);
644 rvec_dec(dx, x_solvt[i]);
646 atoms->atom[atoms->nr] = atoms_solvt->atom[i];
647 atoms->atomname[atoms->nr] = atoms_solvt->atomname[i];
648 rvec_add(x_solvt[i], dx, (*x)[atoms->nr]);
651 copy_rvec(v_solvt[i], (*v)[atoms->nr]);
653 (*r)[atoms->nr] = r_solvt[i];
654 atoms->atom[atoms->nr].resind = atoms->nres-1;
661 srenew(atoms->resinfo, atoms->nres+nresadd);
666 fprintf(stderr, "Added %d molecules\n", nresadd);
670 done_atom(atoms_all);