src/gromacs/gmxpreprocess/addconf.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "addconf.h"
  40
  41 #include "config.h"
  42
  43 #include <stdlib.h>
  44 #include <string.h>
  45
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/legacyheaders/macros.h"
  48 #include "gromacs/legacyheaders/types/commrec.h"
  49 #include "gromacs/legacyheaders/force.h"
  50 #include "gromacs/legacyheaders/names.h"
  51 #include "gromacs/legacyheaders/nsgrid.h"
  52 #include "gromacs/legacyheaders/mdatoms.h"
  53 #include "gromacs/legacyheaders/nrnb.h"
  54 #include "gromacs/legacyheaders/ns.h"
  55 #include "gromacs/topology/mtop_util.h"
  56 #include "gromacs/legacyheaders/chargegroup.h"
  57
  58 #include "gromacs/pbcutil/pbc.h"
  59 #include "gromacs/topology/block.h"
  60 #include "gromacs/topology/topology.h"
  61 #include "gromacs/utility/smalloc.h"
  62
  63 static real box_margin;
  64
  65 static real max_dist(rvec *x, real *r, int start, int end)
  66 {
  67     real maxd;
  68     int  i, j;
  69
  70     maxd = 0;
  71     for (i = start; i < end; i++)
  72     {
  73         for (j = i+1; j < end; j++)
  74         {
  75             maxd = max(maxd, sqrt(distance2(x[i], x[j]))+0.5*(r[i]+r[j]));
  76         }
  77     }
  78
  79     return 0.5*maxd;
  80 }
  81
  82 static gmx_bool outside_box_minus_margin2(rvec x, matrix box)
  83 {
  84     return ( (x[XX] < 2*box_margin) || (x[XX] > box[XX][XX]-2*box_margin) ||
  85              (x[YY] < 2*box_margin) || (x[YY] > box[YY][YY]-2*box_margin) ||
  86              (x[ZZ] < 2*box_margin) || (x[ZZ] > box[ZZ][ZZ]-2*box_margin) );
  87 }
  88
  89 static gmx_bool outside_box_plus_margin(rvec x, matrix box)
  90 {
  91     return ( (x[XX] < -box_margin) || (x[XX] > box[XX][XX]+box_margin) ||
  92              (x[YY] < -box_margin) || (x[YY] > box[YY][YY]+box_margin) ||
  93              (x[ZZ] < -box_margin) || (x[ZZ] > box[ZZ][ZZ]+box_margin) );
  94 }
  95
  96 static int mark_res(int at, gmx_bool *mark, int natoms, t_atom *atom, int *nmark)
  97 {
  98     int resind;
  99
 100     resind = atom[at].resind;
 101     while ( (at > 0) && (resind == atom[at-1].resind) )
 102     {
 103         at--;
 104     }
 105     while ( (at < natoms) && (resind == atom[at].resind) )
 106     {
 107         if (!mark[at])
 108         {
 109             mark[at] = TRUE;
 110             (*nmark)++;
 111         }
 112         at++;
 113     }
 114
 115     return at;
 116 }
 117
 118 static real find_max_real(int n, real radius[])
 119 {
 120     int  i;
 121     real rmax;
 122
 123     rmax = 0;
 124     if (n > 0)
 125     {
 126         rmax = radius[0];
 127         for (i = 1; (i < n); i++)
 128         {
 129             rmax = max(rmax, radius[i]);
 130         }
 131     }
 132     return rmax;
 133 }
 134
 135 static void combine_atoms(t_atoms *ap, t_atoms *as,
 136                           rvec xp[], rvec *vp, rvec xs[], rvec *vs,
 137                           t_atoms **a_comb, rvec **x_comb, rvec **v_comb)
 138 {
 139     t_atoms *ac;
 140     rvec    *xc, *vc = NULL;
 141     int      i, j, natot, res0;
 142
 143     /* Total number of atoms */
 144     natot = ap->nr+as->nr;
 145
 146     snew(ac, 1);
 147     init_t_atoms(ac, natot, FALSE);
 148
 149     snew(xc, natot);
 150     if (vp && vs)
 151     {
 152         snew(vc, natot);
 153     }
 154
 155     /* Fill the new structures */
 156     for (i = j = 0; (i < ap->nr); i++, j++)
 157     {
 158         copy_rvec(xp[i], xc[j]);
 159         if (vc)
 160         {
 161             copy_rvec(vp[i], vc[j]);
 162         }
 163         memcpy(&(ac->atom[j]), &(ap->atom[i]), sizeof(ap->atom[i]));
 164         ac->atom[j].type = 0;
 165     }
 166     res0 = ap->nres;
 167     for (i = 0; (i < as->nr); i++, j++)
 168     {
 169         copy_rvec(xs[i], xc[j]);
 170         if (vc)
 171         {
 172             copy_rvec(vs[i], vc[j]);
 173         }
 174         memcpy(&(ac->atom[j]), &(as->atom[i]), sizeof(as->atom[i]));
 175         ac->atom[j].type    = 0;
 176         ac->atom[j].resind += res0;
 177     }
 178     ac->nr   = j;
 179     ac->nres = ac->atom[j-1].resind+1;
 180     /* Fill all elements to prevent uninitialized memory */
 181     for (i = 0; i < ac->nr; i++)
 182     {
 183         ac->atom[i].m     = 1;
 184         ac->atom[i].q     = 0;
 185         ac->atom[i].mB    = 1;
 186         ac->atom[i].qB    = 0;
 187         ac->atom[i].type  = 0;
 188         ac->atom[i].typeB = 0;
 189         ac->atom[i].ptype = eptAtom;
 190     }
 191
 192     /* Return values */
 193     *a_comb = ac;
 194     *x_comb = xc;
 195     *v_comb = vc;
 196 }
 197
 198 static t_forcerec *fr = NULL;
 199
 200 static void do_nsgrid(FILE *fp, gmx_bool bVerbose,
 201                       matrix box, rvec x[], t_atoms *atoms, real rlong,
 202                       const output_env_t oenv)
 203 {
 204     gmx_mtop_t     *mtop;
 205     gmx_localtop_t *top;
 206     t_mdatoms      *md;
 207     t_block        *cgs;
 208     t_inputrec     *ir;
 209     t_nrnb          nrnb;
 210     t_commrec      *cr;
 211     int            *cg_index;
 212     gmx_moltype_t  *molt;
 213     gmx_ffparams_t *ffp;
 214     ivec           *nFreeze;
 215     int             i, m, natoms;
 216     rvec            box_size;
 217     real           *lambda, *dvdl;
 218
 219     natoms = atoms->nr;
 220
 221     /* Charge group index */
 222     snew(cg_index, natoms);
 223     for (i = 0; (i < natoms); i++)
 224     {
 225         cg_index[i] = i;
 226     }
 227
 228     /* Topology needs charge groups and exclusions */
 229     snew(mtop, 1);
 230     init_mtop(mtop);
 231     mtop->natoms = natoms;
 232     /* Make one moltype that contains the whol system */
 233     mtop->nmoltype = 1;
 234     snew(mtop->moltype, mtop->nmoltype);
 235     molt        = &mtop->moltype[0];
 236     molt->name  = mtop->name;
 237     molt->atoms = *atoms;
 238     stupid_fill_block(&molt->cgs, mtop->natoms, FALSE);
 239     stupid_fill_blocka(&molt->excls, natoms);
 240     /* Make one molblock for the whole system */
 241     mtop->nmolblock = 1;
 242     snew(mtop->molblock, mtop->nmolblock);
 243     mtop->molblock[0].type       = 0;
 244     mtop->molblock[0].nmol       = 1;
 245     mtop->molblock[0].natoms_mol = natoms;
 246     /* Initialize a single energy group */
 247     mtop->groups.grps[egcENER].nr = 1;
 248     mtop->groups.ngrpnr[egcENER]  = 0;
 249     mtop->groups.grpnr[egcENER]   = NULL;
 250
 251     ffp = &mtop->ffparams;
 252
 253     ffp->ntypes = 1;
 254     ffp->atnr   = 1;
 255     ffp->reppow = 12;
 256     snew(ffp->functype, 1);
 257     snew(ffp->iparams, 1);
 258     ffp->iparams[0].lj.c6  = 1;
 259     ffp->iparams[0].lj.c12 = 1;
 260
 261     /* inputrec structure */
 262     snew(ir, 1);
 263     ir->cutoff_scheme    = ecutsGROUP;
 264     ir->coulomb_modifier = eintmodNONE;
 265     ir->vdw_modifier     = eintmodNONE;
 266     ir->coulombtype      = eelCUT;
 267     ir->vdwtype          = evdwCUT;
 268     ir->ns_type          = ensGRID;
 269     snew(ir->opts.egp_flags, 1);
 270
 271     top = gmx_mtop_generate_local_top(mtop, ir);
 272
 273     /* Some nasty shortcuts */
 274     cgs  = &(top->cgs);
 275
 276     /* mdatoms structure */
 277     snew(nFreeze, 2);
 278     snew(md, 1);
 279     md = init_mdatoms(fp, mtop, FALSE);
 280     atoms2md(mtop, ir, 0, NULL, mtop->natoms, md);
 281     sfree(nFreeze);
 282
 283     /* forcerec structure */
 284     if (fr == NULL)
 285     {
 286         fr = mk_forcerec();
 287     }
 288     snew(cr, 1);
 289     cr->nnodes   = 1;
 290     /* cr->nthreads = 1; */
 291
 292     /*    ir->rlist       = ir->rcoulomb = ir->rvdw = rlong;
 293        printf("Neighborsearching with a cut-off of %g\n",rlong);
 294        init_forcerec(stdout,fr,ir,top,cr,md,box,FALSE,NULL,NULL,NULL,TRUE);*/
 295     fr->cg0     = 0;
 296     fr->hcg     = top->cgs.nr;
 297     fr->nWatMol = 0;
 298
 299     /* Prepare for neighboursearching */
 300     init_nrnb(&nrnb);
 301
 302     /* Init things dependent on parameters */
 303     ir->rlistlong = ir->rlist = ir->rcoulomb = ir->rvdw = rlong;
 304     /* create free energy data to avoid NULLs */
 305     snew(ir->fepvals, 1);
 306     printf("Neighborsearching with a cut-off of %g\n", rlong);
 307     init_forcerec(stdout, oenv, fr, NULL, ir, mtop, cr, box,
 308                   NULL, NULL, NULL, NULL, NULL, TRUE, -1);
 309     if (debug)
 310     {
 311         pr_forcerec(debug, fr);
 312     }
 313
 314     /* Calculate new stuff dependent on coords and box */
 315     for (m = 0; (m < DIM); m++)
 316     {
 317         box_size[m] = box[m][m];
 318     }
 319     calc_shifts(box, fr->shift_vec);
 320     put_charge_groups_in_box(fp, 0, cgs->nr, fr->ePBC, box, cgs, x, fr->cg_cm);
 321
 322     /* Do the actual neighboursearching */
 323     snew(lambda, efptNR);
 324     snew(dvdl, efptNR);
 325     init_neighbor_list(fp, fr, md->homenr);
 326     search_neighbours(fp, fr, box, top,
 327                       &mtop->groups, cr, &nrnb, md, TRUE, FALSE);
 328
 329     if (debug)
 330     {
 331         dump_nblist(debug, cr, fr, 0);
 332     }
 333
 334     if (bVerbose)
 335     {
 336         fprintf(stderr, "Successfully made neighbourlist\n");
 337     }
 338 }
 339
 340 static gmx_bool bXor(gmx_bool b1, gmx_bool b2)
 341 {
 342     return (b1 && !b2) || (b2 && !b1);
 343 }
 344
 345 void add_conf(t_atoms *atoms, rvec **x, rvec **v, real **r, gmx_bool bSrenew,
 346               int ePBC, matrix box, gmx_bool bInsert,
 347               t_atoms *atoms_solvt, rvec *x_solvt, rvec *v_solvt, real *r_solvt,
 348               gmx_bool bVerbose, real rshell, int max_sol, const output_env_t oenv)
 349 {
 350     t_nblist       *nlist;
 351     t_atoms        *atoms_all;
 352     real            max_vdw, *r_prot, *r_all, n2, r2, ib1, ib2;
 353     int             natoms_prot, natoms_solvt;
 354     int             i, j, jj, m, j0, j1, jjj, jnres, jnr, inr, iprot, is1, is2;
 355     int             prev, resnr, nresadd, d, k, ncells, maxincell;
 356     int             dx0, dx1, dy0, dy1, dz0, dz1;
 357     int             ntest, nremove, nkeep;
 358     rvec            dx, xi, xj, xpp, *x_all, *v_all;
 359     gmx_bool       *remove, *keep;
 360     int             bSolSol;
 361
 362     natoms_prot  = atoms->nr;
 363     natoms_solvt = atoms_solvt->nr;
 364     if (natoms_solvt <= 0)
 365     {
 366         fprintf(stderr, "WARNING: Nothing to add\n");
 367         return;
 368     }
 369
 370     if (ePBC == epbcSCREW)
 371     {
 372         gmx_fatal(FARGS, "Sorry, %s pbc is not yet supported", epbc_names[ePBC]);
 373     }
 374
 375     if (bVerbose)
 376     {
 377         fprintf(stderr, "Calculating Overlap...\n");
 378     }
 379
 380     /* Set margin around box edges to largest solvent dimension.
 381      * The maximum distance between atoms in a solvent molecule should
 382      * be calculated. At the moment a fudge factor of 3 is used.
 383      */
 384     r_prot     = *r;
 385     box_margin = 3*find_max_real(natoms_solvt, r_solvt);
 386     max_vdw    = max(3*find_max_real(natoms_prot, r_prot), box_margin);
 387     fprintf(stderr, "box_margin = %g\n", box_margin);
 388
 389     snew(remove, natoms_solvt);
 390
 391     nremove = 0;
 392     if (!bInsert)
 393     {
 394         for (i = 0; i < atoms_solvt->nr; i++)
 395         {
 396             if (outside_box_plus_margin(x_solvt[i], box) )
 397             {
 398                 i = mark_res(i, remove, atoms_solvt->nr, atoms_solvt->atom, &nremove);
 399             }
 400         }
 401         fprintf(stderr, "Removed %d atoms that were outside the box\n", nremove);
 402     }
 403
 404     /* Define grid stuff */
 405     /* Largest VDW radius */
 406     snew(r_all, natoms_prot+natoms_solvt);
 407     for (i = j = 0; i < natoms_prot; i++, j++)
 408     {
 409         r_all[j] = r_prot[i];
 410     }
 411     for (i = 0; i < natoms_solvt; i++, j++)
 412     {
 413         r_all[j] = r_solvt[i];
 414     }
 415
 416     /* Combine arrays */
 417     combine_atoms(atoms, atoms_solvt, *x, v ? *v : NULL, x_solvt, v_solvt,
 418                   &atoms_all, &x_all, &v_all);
 419
 420     /* Do neighboursearching step */
 421     do_nsgrid(stdout, bVerbose, box, x_all, atoms_all, max_vdw, oenv);
 422
 423     /* check solvent with solute */
 424     nlist = &(fr->nblists[0].nlist_sr[eNL_VDW]);
 425     fprintf(stderr, "nri = %d, nrj = %d\n", nlist->nri, nlist->nrj);
 426     for (bSolSol = 0; (bSolSol <= (bInsert ? 0 : 1)); bSolSol++)
 427     {
 428         ntest = nremove = 0;
 429         fprintf(stderr, "Checking %s-Solvent overlap:",
 430                 bSolSol ? "Solvent" : "Protein");
 431         for (i = 0; (i < nlist->nri && nremove < natoms_solvt); i++)
 432         {
 433             inr = nlist->iinr[i];
 434             j0  = nlist->jindex[i];
 435             j1  = nlist->jindex[i+1];
 436             rvec_add(x_all[inr], fr->shift_vec[nlist->shift[i]], xi);
 437
 438             for (j = j0; (j < j1 && nremove < natoms_solvt); j++)
 439             {
 440                 jnr = nlist->jjnr[j];
 441                 if (jnr < 0)
 442                 {
 443                     /* skip padding */
 444                     continue;
 445                 }
 446                 copy_rvec(x_all[jnr], xj);
 447
 448                 /* Check solvent-protein and solvent-solvent */
 449                 is1 = inr-natoms_prot;
 450                 is2 = jnr-natoms_prot;
 451
 452                 /* Check if at least one of the atoms is a solvent that is not yet
 453                  * listed for removal, and if both are solvent, that they are not in the
 454                  * same residue.
 455                  */
 456                 if ((!bSolSol &&
 457                      bXor((is1 >= 0), (is2 >= 0)) && /* One atom is protein */
 458                      ((is1 < 0) || ((is1 >= 0) && !remove[is1])) &&
 459                      ((is2 < 0) || ((is2 >= 0) && !remove[is2]))) ||
 460
 461                     (bSolSol  &&
 462                      (is1 >= 0) && (!remove[is1]) &&                   /* is1 is solvent */
 463                      (is2 >= 0) && (!remove[is2]) &&                   /* is2 is solvent */
 464                      (bInsert ||                                       /* when inserting also check inside the box */
 465                       (outside_box_minus_margin2(x_solvt[is1], box) && /* is1 on edge */
 466                        outside_box_minus_margin2(x_solvt[is2], box))   /* is2 on edge */
 467                      ) &&
 468                      (atoms_solvt->atom[is1].resind !=                 /* Not the same residue */
 469                       atoms_solvt->atom[is2].resind)))
 470                 {
 471
 472                     ntest++;
 473                     rvec_sub(xi, xj, dx);
 474                     n2 = norm2(dx);
 475                     r2 = sqr(r_all[inr]+r_all[jnr]);
 476                     if (n2 < r2)
 477                     {
 478                         if (bInsert)
 479                         {
 480                             nremove = natoms_solvt;
 481                             for (k = 0; k < nremove; k++)
 482                             {
 483                                 remove[k] = TRUE;
 484                             }
 485                         }
 486                         /* Need only remove one of the solvents... */
 487                         if (is2 >= 0)
 488                         {
 489                             (void) mark_res(is2, remove, natoms_solvt, atoms_solvt->atom,
 490                                             &nremove);
 491                         }
 492                         else if (is1 >= 0)
 493                         {
 494                             (void) mark_res(is1, remove, natoms_solvt, atoms_solvt->atom,
 495                                             &nremove);
 496                         }
 497                         else
 498                         {
 499                             fprintf(stderr, "Neither atom is solvent%d %d\n", is1, is2);
 500                         }
 501                     }
 502                 }
 503             }
 504         }
 505         if (!bInsert)
 506         {
 507             fprintf(stderr, " tested %d pairs, removed %d atoms.\n", ntest, nremove);
 508         }
 509     }
 510     if (debug)
 511     {
 512         for (i = 0; i < natoms_solvt; i++)
 513         {
 514             fprintf(debug, "remove[%5d] = %s\n", i, bool_names[remove[i]]);
 515         }
 516     }
 517
 518     /* Search again, now with another cut-off */
 519     if (rshell > 0)
 520     {
 521         do_nsgrid(stdout, bVerbose, box, x_all, atoms_all, rshell, oenv);
 522         nlist = &(fr->nblists[0].nlist_sr[eNL_VDW]);
 523         fprintf(stderr, "nri = %d, nrj = %d\n", nlist->nri, nlist->nrj);
 524         nkeep = 0;
 525         snew(keep, natoms_solvt);
 526         for (i = 0; i < nlist->nri; i++)
 527         {
 528             inr = nlist->iinr[i];
 529             j0  = nlist->jindex[i];
 530             j1  = nlist->jindex[i+1];
 531
 532             for (j = j0; j < j1; j++)
 533             {
 534                 jnr = nlist->jjnr[j];
 535
 536                 /* Check solvent-protein and solvent-solvent */
 537                 is1 = inr-natoms_prot;
 538                 is2 = jnr-natoms_prot;
 539
 540                 /* Check if at least one of the atoms is a solvent that is not yet
 541                  * listed for removal, and if both are solvent, that they are not in the
 542                  * same residue.
 543                  */
 544                 if (is1 >= 0 && is2 < 0)
 545                 {
 546                     mark_res(is1, keep, natoms_solvt, atoms_solvt->atom, &nkeep);
 547                 }
 548                 else if (is1 < 0 && is2 >= 0)
 549                 {
 550                     mark_res(is2, keep, natoms_solvt, atoms_solvt->atom, &nkeep);
 551                 }
 552             }
 553         }
 554         fprintf(stderr, "Keeping %d solvent atoms after proximity check\n",
 555                 nkeep);
 556         for (i = 0; i < natoms_solvt; i++)
 557         {
 558             remove[i] = remove[i] || !keep[i];
 559         }
 560         sfree(keep);
 561     }
 562     /* count how many atoms and residues will be added and make space */
 563     if (bInsert)
 564     {
 565         j     = atoms_solvt->nr;
 566         jnres = atoms_solvt->nres;
 567     }
 568     else
 569     {
 570         j     = 0;
 571         jnres = 0;
 572         for (i = 0; ((i < atoms_solvt->nr) &&
 573                      ((max_sol == 0) || (jnres < max_sol))); i++)
 574         {
 575             if (!remove[i])
 576             {
 577                 j++;
 578                 if ((i == 0) ||
 579                     (atoms_solvt->atom[i].resind != atoms_solvt->atom[i-1].resind))
 580                 {
 581                     jnres++;
 582                 }
 583             }
 584         }
 585     }
 586     if (debug)
 587     {
 588         fprintf(debug, "Will add %d atoms in %d residues\n", j, jnres);
 589     }
 590     if (!bInsert)
 591     {
 592         /* Flag the remaing solvent atoms to be removed */
 593         jjj = atoms_solvt->atom[i-1].resind;
 594         for (; (i < atoms_solvt->nr); i++)
 595         {
 596             if (atoms_solvt->atom[i].resind > jjj)
 597             {
 598                 remove[i] = TRUE;
 599             }
 600             else
 601             {
 602                 j++;
 603             }
 604         }
 605     }
 606
 607     if (bSrenew)
 608     {
 609         srenew(atoms->resinfo,  atoms->nres+jnres);
 610         srenew(atoms->atomname, atoms->nr+j);
 611         srenew(atoms->atom,     atoms->nr+j);
 612         srenew(*x,              atoms->nr+j);
 613         if (v)
 614         {
 615             srenew(*v,       atoms->nr+j);
 616         }
 617         srenew(*r,              atoms->nr+j);
 618     }
 619
 620     /* add the selected atoms_solvt to atoms */
 621     if (atoms->nr > 0)
 622     {
 623         resnr = atoms->resinfo[atoms->atom[atoms->nr-1].resind].nr;
 624     }
 625     else
 626     {
 627         resnr = 0;
 628     }
 629     prev    = -1;
 630     nresadd = 0;
 631     for (i = 0; i < atoms_solvt->nr; i++)
 632     {
 633         if (!remove[i])
 634         {
 635             if (prev == -1 ||
 636                 atoms_solvt->atom[i].resind != atoms_solvt->atom[prev].resind)
 637             {
 638                 nresadd++;
 639                 atoms->nres++;
 640                 resnr++;
 641                 atoms->resinfo[atoms->nres-1] =
 642                     atoms_solvt->resinfo[atoms_solvt->atom[i].resind];
 643                 atoms->resinfo[atoms->nres-1].nr = resnr;
 644                 /* calculate shift of the solvent molecule using the first atom */
 645                 copy_rvec(x_solvt[i], dx);
 646                 put_atoms_in_box(ePBC, box, 1, &dx);
 647                 rvec_dec(dx, x_solvt[i]);
 648             }
 649             atoms->atom[atoms->nr]     = atoms_solvt->atom[i];
 650             atoms->atomname[atoms->nr] = atoms_solvt->atomname[i];
 651             rvec_add(x_solvt[i], dx, (*x)[atoms->nr]);
 652             if (v)
 653             {
 654                 copy_rvec(v_solvt[i], (*v)[atoms->nr]);
 655             }
 656             (*r)[atoms->nr]               = r_solvt[i];
 657             atoms->atom[atoms->nr].resind = atoms->nres-1;
 658             atoms->nr++;
 659             prev = i;
 660         }
 661     }
 662     if (bSrenew)
 663     {
 664         srenew(atoms->resinfo,  atoms->nres+nresadd);
 665     }
 666
 667     if (bVerbose)
 668     {
 669         fprintf(stderr, "Added %d molecules\n", nresadd);
 670     }
 671
 672     sfree(remove);
 673     done_atom(atoms_all);
 674     sfree(x_all);
 675     sfree(v_all);
 676 }