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