src/gromacs/gmxpreprocess/insert-molecules.cpp

   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,2015,2016, 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 "insert-molecules.h"
  40
  41 #include <algorithm>
  42 #include <memory>
  43 #include <set>
  44 #include <string>
  45 #include <vector>
  46
  47 #include "gromacs/commandline/cmdlineoptionsmodule.h"
  48 #include "gromacs/fileio/confio.h"
  49 #include "gromacs/fileio/filetypes.h"
  50 #include "gromacs/fileio/xvgr.h"
  51 #include "gromacs/gmxlib/conformation-utilities.h"
  52 #include "gromacs/gmxpreprocess/read-conformation.h"
  53 #include "gromacs/math/functions.h"
  54 #include "gromacs/math/utilities.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/options/basicoptions.h"
  57 #include "gromacs/options/filenameoption.h"
  58 #include "gromacs/options/ioptionscontainer.h"
  59 #include "gromacs/pbcutil/pbc.h"
  60 #include "gromacs/random/threefry.h"
  61 #include "gromacs/random/uniformrealdistribution.h"
  62 #include "gromacs/selection/nbsearch.h"
  63 #include "gromacs/selection/selection.h"
  64 #include "gromacs/selection/selectioncollection.h"
  65 #include "gromacs/selection/selectionoption.h"
  66 #include "gromacs/selection/selectionoptionbehavior.h"
  67 #include "gromacs/topology/atomprop.h"
  68 #include "gromacs/topology/atoms.h"
  69 #include "gromacs/topology/atomsbuilder.h"
  70 #include "gromacs/topology/mtop_util.h"
  71 #include "gromacs/topology/topology.h"
  72 #include "gromacs/trajectory/trajectoryframe.h"
  73 #include "gromacs/utility/cstringutil.h"
  74 #include "gromacs/utility/exceptions.h"
  75 #include "gromacs/utility/fatalerror.h"
  76 #include "gromacs/utility/smalloc.h"
  77
  78 using gmx::RVec;
  79
  80 /* enum for random rotations of inserted solutes */
  81 enum RotationType {
  82     en_rotXYZ, en_rotZ, en_rotNone
  83 };
  84 const char *const cRotationEnum[] = {"xyz", "z", "none"};
  85
  86 static void center_molecule(std::vector<RVec> *x)
  87 {
  88     const size_t atomCount = x->size();
  89     rvec         center;
  90     clear_rvec(center);
  91     for (size_t i = 0; i < atomCount; ++i)
  92     {
  93         rvec_inc(center, (*x)[i]);
  94     }
  95     svmul(1.0/atomCount, center, center);
  96     for (size_t i = 0; i < atomCount; ++i)
  97     {
  98         rvec_dec((*x)[i], center);
  99     }
 100 }
 101
 102 static void generate_trial_conf(const std::vector<RVec> &xin,
 103                                 const rvec offset, RotationType enum_rot,
 104                                 gmx::DefaultRandomEngine * rng,
 105                                 std::vector<RVec> *xout)
 106 {
 107     gmx::UniformRealDistribution<real> dist(0, 2.0*M_PI);
 108     *xout = xin;
 109
 110     real alfa = 0.0, beta = 0.0, gamma = 0.0;
 111     switch (enum_rot)
 112     {
 113         case en_rotXYZ:
 114             alfa  = dist(*rng);
 115             beta  = dist(*rng);
 116             gamma = dist(*rng);
 117             break;
 118         case en_rotZ:
 119             alfa  = beta = 0.;
 120             gamma = dist(*rng);
 121             break;
 122         case en_rotNone:
 123             alfa = beta = gamma = 0.;
 124             break;
 125     }
 126     if (enum_rot == en_rotXYZ || enum_rot == en_rotZ)
 127     {
 128         rotate_conf(xout->size(), as_rvec_array(xout->data()), NULL, alfa, beta, gamma);
 129     }
 130     for (size_t i = 0; i < xout->size(); ++i)
 131     {
 132         rvec_inc((*xout)[i], offset);
 133     }
 134 }
 135
 136 static bool isInsertionAllowed(gmx::AnalysisNeighborhoodSearch *search,
 137                                const std::vector<real>         &exclusionDistances,
 138                                const std::vector<RVec>         &x,
 139                                const std::vector<real>         &exclusionDistances_insrt,
 140                                const t_atoms                   &atoms,
 141                                const std::set<int>             &removableAtoms,
 142                                gmx::AtomsRemover               *remover)
 143 {
 144     gmx::AnalysisNeighborhoodPositions  pos(x);
 145     gmx::AnalysisNeighborhoodPairSearch pairSearch = search->startPairSearch(pos);
 146     gmx::AnalysisNeighborhoodPair       pair;
 147     while (pairSearch.findNextPair(&pair))
 148     {
 149         const real r1 = exclusionDistances[pair.refIndex()];
 150         const real r2 = exclusionDistances_insrt[pair.testIndex()];
 151         if (pair.distance2() < gmx::square(r1 + r2))
 152         {
 153             if (removableAtoms.count(pair.refIndex()) == 0)
 154             {
 155                 return false;
 156             }
 157             // TODO: If molecule information is available, this should ideally
 158             // use it to remove whole molecules.
 159             remover->markResidue(atoms, pair.refIndex(), true);
 160         }
 161     }
 162     return true;
 163 }
 164
 165 static void insert_mols(int nmol_insrt, int ntry, int seed,
 166                         real defaultDistance, real scaleFactor,
 167                         t_atoms *atoms, t_symtab *symtab, std::vector<RVec> *x,
 168                         const std::set<int> &removableAtoms,
 169                         const t_atoms &atoms_insrt, const std::vector<RVec> &x_insrt,
 170                         int ePBC, matrix box,
 171                         const std::string &posfn, const rvec deltaR,
 172                         RotationType enum_rot)
 173 {
 174     fprintf(stderr, "Initialising inter-atomic distances...\n");
 175     gmx_atomprop_t          aps = gmx_atomprop_init();
 176     std::vector<real>       exclusionDistances(
 177             makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor));
 178     const std::vector<real> exclusionDistances_insrt(
 179             makeExclusionDistances(&atoms_insrt, aps, defaultDistance, scaleFactor));
 180     gmx_atomprop_destroy(aps);
 181
 182     const real       maxInsertRadius
 183         = *std::max_element(exclusionDistances_insrt.begin(),
 184                             exclusionDistances_insrt.end());
 185     real             maxRadius = maxInsertRadius;
 186     if (!exclusionDistances.empty())
 187     {
 188         const real maxExistingRadius
 189             = *std::max_element(exclusionDistances.begin(),
 190                                 exclusionDistances.end());
 191         maxRadius = std::max(maxInsertRadius, maxExistingRadius);
 192     }
 193
 194     // TODO: Make all of this exception-safe.
 195     gmx::AnalysisNeighborhood nb;
 196     nb.setCutoff(maxInsertRadius + maxRadius);
 197
 198
 199     if (seed == 0)
 200     {
 201         seed = static_cast<int>(gmx::makeRandomSeed());
 202     }
 203     fprintf(stderr, "Using random seed %d\n", seed);
 204
 205     gmx::DefaultRandomEngine rng(seed);
 206
 207     t_pbc                    pbc;
 208     set_pbc(&pbc, ePBC, box);
 209
 210     /* With -ip, take nmol_insrt from file posfn */
 211     double     **rpos              = NULL;
 212     const bool   insertAtPositions = !posfn.empty();
 213     if (insertAtPositions)
 214     {
 215         int ncol;
 216         nmol_insrt = read_xvg(posfn.c_str(), &rpos, &ncol);
 217         if (ncol != 3)
 218         {
 219             gmx_fatal(FARGS, "Expected 3 columns (x/y/z coordinates) in file %s\n",
 220                       posfn.c_str());
 221         }
 222         fprintf(stderr, "Read %d positions from file %s\n\n",
 223                 nmol_insrt, posfn.c_str());
 224     }
 225
 226     gmx::AtomsBuilder builder(atoms, symtab);
 227     gmx::AtomsRemover remover(*atoms);
 228     {
 229         const int finalAtomCount    = atoms->nr + nmol_insrt * atoms_insrt.nr;
 230         const int finalResidueCount = atoms->nres + nmol_insrt * atoms_insrt.nres;
 231         builder.reserve(finalAtomCount, finalResidueCount);
 232         x->reserve(finalAtomCount);
 233         exclusionDistances.reserve(finalAtomCount);
 234     }
 235
 236     std::vector<RVec>                    x_n(x_insrt.size());
 237
 238     int                                  mol        = 0;
 239     int                                  trial      = 0;
 240     int                                  firstTrial = 0;
 241     int                                  failed     = 0;
 242     gmx::UniformRealDistribution<real>   dist;
 243
 244     while (mol < nmol_insrt && trial < ntry*nmol_insrt)
 245     {
 246         // cppcheck 1.72 complains about uninitialized variables in the
 247         // assignments below otherwise...
 248         rvec offset_x = {0};
 249         if (!insertAtPositions)
 250         {
 251             // Insert at random positions.
 252             offset_x[XX] = box[XX][XX] * dist(rng);
 253             offset_x[YY] = box[YY][YY] * dist(rng);
 254             offset_x[ZZ] = box[ZZ][ZZ] * dist(rng);
 255         }
 256         else
 257         {
 258             // Skip a position if ntry trials were not successful.
 259             if (trial >= firstTrial + ntry)
 260             {
 261                 fprintf(stderr, " skipped position (%.3f, %.3f, %.3f)\n",
 262                         rpos[XX][mol], rpos[YY][mol], rpos[ZZ][mol]);
 263                 ++mol;
 264                 ++failed;
 265                 firstTrial = trial;
 266                 continue;
 267             }
 268             // Insert at positions taken from option -ip file.
 269             offset_x[XX] = rpos[XX][mol] + deltaR[XX]*(2 * dist(rng)-1);
 270             offset_x[YY] = rpos[YY][mol] + deltaR[YY]*(2 * dist(rng)-1);
 271             offset_x[ZZ] = rpos[ZZ][mol] + deltaR[ZZ]*(2 * dist(rng)-1);
 272         }
 273         fprintf(stderr, "\rTry %d", ++trial);
 274         fflush(stderr);
 275
 276         generate_trial_conf(x_insrt, offset_x, enum_rot, &rng, &x_n);
 277         gmx::AnalysisNeighborhoodPositions pos(*x);
 278         gmx::AnalysisNeighborhoodSearch    search = nb.initSearch(&pbc, pos);
 279         if (isInsertionAllowed(&search, exclusionDistances, x_n, exclusionDistances_insrt,
 280                                *atoms, removableAtoms, &remover))
 281         {
 282             x->insert(x->end(), x_n.begin(), x_n.end());
 283             exclusionDistances.insert(exclusionDistances.end(),
 284                                       exclusionDistances_insrt.begin(),
 285                                       exclusionDistances_insrt.end());
 286             builder.mergeAtoms(atoms_insrt);
 287             ++mol;
 288             firstTrial = trial;
 289             fprintf(stderr, " success (now %d atoms)!\n", builder.currentAtomCount());
 290         }
 291     }
 292
 293     fprintf(stderr, "\n");
 294     /* print number of molecules added */
 295     fprintf(stderr, "Added %d molecules (out of %d requested)\n",
 296             mol - failed, nmol_insrt);
 297
 298     const int originalAtomCount    = atoms->nr;
 299     const int originalResidueCount = atoms->nres;
 300     remover.refreshAtomCount(*atoms);
 301     remover.removeMarkedElements(x);
 302     remover.removeMarkedAtoms(atoms);
 303     if (atoms->nr < originalAtomCount)
 304     {
 305         fprintf(stderr, "Replaced %d residues (%d atoms)\n",
 306                 originalResidueCount - atoms->nres,
 307                 originalAtomCount - atoms->nr);
 308     }
 309
 310     if (rpos != NULL)
 311     {
 312         for (int i = 0; i < DIM; ++i)
 313         {
 314             sfree(rpos[i]);
 315         }
 316         sfree(rpos);
 317     }
 318 }
 319
 320 namespace gmx
 321 {
 322
 323 namespace
 324 {
 325
 326 class InsertMolecules : public ICommandLineOptionsModule, public ITopologyProvider
 327 {
 328     public:
 329         InsertMolecules()
 330             : bBox_(false), nmolIns_(0), nmolTry_(10), seed_(0),
 331               defaultDistance_(0.105), scaleFactor_(0.57), enumRot_(en_rotXYZ),
 332               top_(NULL), ePBC_(-1)
 333         {
 334             clear_rvec(newBox_);
 335             clear_rvec(deltaR_);
 336             clear_mat(box_);
 337         }
 338         virtual ~InsertMolecules()
 339         {
 340             if (top_ != NULL)
 341             {
 342                 done_mtop(top_);
 343                 sfree(top_);
 344             }
 345         }
 346
 347         // From ITopologyProvider
 348         virtual gmx_mtop_t *getTopology(bool /*required*/) { return top_; }
 349         virtual int getAtomCount() { return 0; }
 350
 351         // From ICommandLineOptionsModule
 352         virtual void init(CommandLineModuleSettings * /*settings*/)
 353         {
 354         }
 355         virtual void initOptions(IOptionsContainer                 *options,
 356                                  ICommandLineOptionsModuleSettings *settings);
 357         virtual void optionsFinished();
 358         virtual int run();
 359
 360     private:
 361         void loadSolute();
 362
 363         SelectionCollection selections_;
 364
 365         std::string         inputConfFile_;
 366         std::string         insertConfFile_;
 367         std::string         positionFile_;
 368         std::string         outputConfFile_;
 369         rvec                newBox_;
 370         bool                bBox_;
 371         int                 nmolIns_;
 372         int                 nmolTry_;
 373         int                 seed_;
 374         real                defaultDistance_;
 375         real                scaleFactor_;
 376         rvec                deltaR_;
 377         RotationType        enumRot_;
 378         Selection           replaceSel_;
 379
 380         gmx_mtop_t         *top_;
 381         std::vector<RVec>   x_;
 382         matrix              box_;
 383         int                 ePBC_;
 384 };
 385
 386 void InsertMolecules::initOptions(IOptionsContainer                 *options,
 387                                   ICommandLineOptionsModuleSettings *settings)
 388 {
 389     const char *const desc[] = {
 390         "[THISMODULE] inserts [TT]-nmol[tt] copies of the system specified in",
 391         "the [TT]-ci[tt] input file. The insertions take place either into",
 392         "vacant space in the solute conformation given with [TT]-f[tt], or",
 393         "into an empty box given by [TT]-box[tt]. Specifying both [TT]-f[tt]",
 394         "and [TT]-box[tt] behaves like [TT]-f[tt], but places a new box",
 395         "around the solute before insertions. Any velocities present are",
 396         "discarded.",
 397         "",
 398         "It is possible to also insert into a solvated configuration and",
 399         "replace solvent atoms with the inserted atoms. To do this, use",
 400         "[TT]-replace[tt] to specify a selection that identifies the atoms",
 401         "that can be replaced. The tool assumes that all molecules in this",
 402         "selection consist of single residues: each residue from this",
 403         "selection that overlaps with the inserted molecules will be removed",
 404         "instead of preventing insertion.",
 405         "",
 406         "By default, the insertion positions are random (with initial seed",
 407         "specified by [TT]-seed[tt]). The program iterates until [TT]-nmol[tt]",
 408         "molecules have been inserted in the box. Molecules are not inserted",
 409         "where the distance between any existing atom and any atom of the",
 410         "inserted molecule is less than the sum based on the van der Waals",
 411         "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
 412         "Waals radii is read by the program, and the resulting radii scaled",
 413         "by [TT]-scale[tt]. If radii are not found in the database, those"
 414         "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].",
 415         "",
 416         "A total of [TT]-nmol[tt] * [TT]-try[tt] insertion attempts are made",
 417         "before giving up. Increase [TT]-try[tt] if you have several small",
 418         "holes to fill. Option [TT]-rot[tt] specifies whether the insertion",
 419         "molecules are randomly oriented before insertion attempts.",
 420         "",
 421         "Alternatively, the molecules can be inserted only at positions defined in",
 422         "positions.dat ([TT]-ip[tt]). That file should have 3 columns (x,y,z),",
 423         "that give the displacements compared to the input molecule position",
 424         "([TT]-ci[tt]). Hence, if that file should contain the absolute",
 425         "positions, the molecule must be centered on (0,0,0) before using",
 426         "[THISMODULE] (e.g. from [gmx-editconf] [TT]-center[tt]).",
 427         "Comments in that file starting with # are ignored. Option [TT]-dr[tt]",
 428         "defines the maximally allowed displacements during insertial trials.",
 429         "[TT]-try[tt] and [TT]-rot[tt] work as in the default mode (see above)."
 430     };
 431
 432     settings->setHelpText(desc);
 433
 434     std::shared_ptr<SelectionOptionBehavior> selectionOptionBehavior(
 435             new SelectionOptionBehavior(&selections_, this));
 436     settings->addOptionsBehavior(selectionOptionBehavior);
 437
 438     // TODO: Replace use of legacyType.
 439     options->addOption(FileNameOption("f")
 440                            .legacyType(efSTX).inputFile()
 441                            .store(&inputConfFile_)
 442                            .defaultBasename("protein")
 443                            .description("Existing configuration to insert into"));
 444     options->addOption(FileNameOption("ci")
 445                            .legacyType(efSTX).inputFile().required()
 446                            .store(&insertConfFile_)
 447                            .defaultBasename("insert")
 448                            .description("Configuration to insert"));
 449     options->addOption(FileNameOption("ip")
 450                            .filetype(eftGenericData).inputFile()
 451                            .store(&positionFile_)
 452                            .defaultBasename("positions")
 453                            .description("Predefined insertion trial positions"));
 454     options->addOption(FileNameOption("o")
 455                            .legacyType(efSTO).outputFile().required()
 456                            .store(&outputConfFile_)
 457                            .defaultBasename("out")
 458                            .description("Output configuration after insertion"));
 459
 460     options->addOption(SelectionOption("replace").onlyAtoms()
 461                            .store(&replaceSel_)
 462                            .description("Atoms that can be removed if overlapping"));
 463     selectionOptionBehavior->initOptions(options);
 464
 465     options->addOption(RealOption("box").vector()
 466                            .store(newBox_).storeIsSet(&bBox_)
 467                            .description("Box size (in nm)"));
 468     options->addOption(IntegerOption("nmol")
 469                            .store(&nmolIns_)
 470                            .description("Number of extra molecules to insert"));
 471     options->addOption(IntegerOption("try")
 472                            .store(&nmolTry_)
 473                            .description("Try inserting [TT]-nmol[tt] times [TT]-try[tt] times"));
 474     options->addOption(IntegerOption("seed")
 475                            .store(&seed_)
 476                            .description("Random generator seed (0 means generate)"));
 477     options->addOption(RealOption("radius")
 478                            .store(&defaultDistance_)
 479                            .description("Default van der Waals distance"));
 480     options->addOption(RealOption("scale")
 481                            .store(&scaleFactor_)
 482                            .description("Scale factor to multiply Van der Waals radii from the database in share/gromacs/top/vdwradii.dat. The default value of 0.57 yields density close to 1000 g/l for proteins in water."));
 483     options->addOption(RealOption("dr").vector()
 484                            .store(deltaR_)
 485                            .description("Allowed displacement in x/y/z from positions in [TT]-ip[tt] file"));
 486     options->addOption(EnumOption<RotationType>("rot").enumValue(cRotationEnum)
 487                            .store(&enumRot_)
 488                            .description("Rotate inserted molecules randomly"));
 489 }
 490
 491 void InsertMolecules::optionsFinished()
 492 {
 493     if (nmolIns_ <= 0 && positionFile_.empty())
 494     {
 495         GMX_THROW(InconsistentInputError("Either -nmol must be larger than 0, "
 496                                          "or positions must be given with -ip."));
 497     }
 498     if (inputConfFile_.empty() && !bBox_)
 499     {
 500         GMX_THROW(InconsistentInputError("When no solute (-f) is specified, "
 501                                          "a box size (-box) must be specified."));
 502     }
 503     if (replaceSel_.isValid() && inputConfFile_.empty())
 504     {
 505         GMX_THROW(InconsistentInputError("Replacement (-replace) only makes sense "
 506                                          "together with an existing configuration (-f)."));
 507     }
 508
 509     snew(top_, 1);
 510     if (!inputConfFile_.empty())
 511     {
 512         readConformation(inputConfFile_.c_str(), top_, &x_, NULL,
 513                          &ePBC_, box_, "solute");
 514         if (top_->natoms == 0)
 515         {
 516             fprintf(stderr, "Note: no atoms in %s\n", inputConfFile_.c_str());
 517         }
 518     }
 519 }
 520
 521 int InsertMolecules::run()
 522 {
 523     std::set<int> removableAtoms;
 524     if (replaceSel_.isValid())
 525     {
 526         t_pbc       pbc;
 527         set_pbc(&pbc, ePBC_, box_);
 528         t_trxframe *fr;
 529         snew(fr, 1);
 530         fr->natoms = x_.size();
 531         fr->bX     = TRUE;
 532         fr->x      = as_rvec_array(x_.data());
 533         selections_.evaluate(fr, &pbc);
 534         sfree(fr);
 535         removableAtoms.insert(replaceSel_.atomIndices().begin(),
 536                               replaceSel_.atomIndices().end());
 537         // TODO: It could be nice to check that removableAtoms contains full
 538         // residues, since we anyways remove whole residues instead of
 539         // individual atoms.
 540     }
 541
 542     if (bBox_)
 543     {
 544         ePBC_ = epbcXYZ;
 545         clear_mat(box_);
 546         box_[XX][XX] = newBox_[XX];
 547         box_[YY][YY] = newBox_[YY];
 548         box_[ZZ][ZZ] = newBox_[ZZ];
 549     }
 550     if (det(box_) == 0)
 551     {
 552         gmx_fatal(FARGS, "Undefined solute box.\nCreate one with gmx editconf "
 553                   "or give explicit -box command line option");
 554     }
 555
 556     t_topology        *top_insrt;
 557     std::vector<RVec>  x_insrt;
 558     snew(top_insrt, 1);
 559     {
 560         int         ePBC_dummy;
 561         matrix      box_dummy;
 562         readConformation(insertConfFile_.c_str(), top_insrt, &x_insrt,
 563                          NULL, &ePBC_dummy, box_dummy, "molecule");
 564         if (top_insrt->atoms.nr == 0)
 565         {
 566             gmx_fatal(FARGS, "No molecule in %s, please check your input",
 567                       insertConfFile_.c_str());
 568         }
 569         if (top_->name == NULL)
 570         {
 571             top_->name = top_insrt->name;
 572         }
 573         if (positionFile_.empty())
 574         {
 575             center_molecule(&x_insrt);
 576         }
 577     }
 578
 579     // TODO: Adapt to use mtop throughout.
 580     t_atoms atoms = gmx_mtop_global_atoms(top_);
 581
 582     /* add nmol_ins molecules of atoms_ins
 583        in random orientation at random place */
 584     insert_mols(nmolIns_, nmolTry_, seed_, defaultDistance_, scaleFactor_,
 585                 &atoms, &top_->symtab, &x_, removableAtoms, top_insrt->atoms, x_insrt,
 586                 ePBC_, box_, positionFile_, deltaR_, enumRot_);
 587
 588     /* write new configuration to file confout */
 589     fprintf(stderr, "Writing generated configuration to %s\n",
 590             outputConfFile_.c_str());
 591     write_sto_conf(outputConfFile_.c_str(), *top_->name, &atoms,
 592                    as_rvec_array(x_.data()), NULL, ePBC_, box_);
 593
 594     /* print size of generated configuration */
 595     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
 596             atoms.nr, atoms.nres);
 597
 598     done_atom(&atoms);
 599     done_top(top_insrt);
 600     sfree(top_insrt);
 601
 602     return 0;
 603 }
 604
 605 }   // namespace
 606
 607 const char InsertMoleculesInfo::name[]             = "insert-molecules";
 608 const char InsertMoleculesInfo::shortDescription[] =
 609     "Insert molecules into existing vacancies";
 610 ICommandLineOptionsModulePointer InsertMoleculesInfo::create()
 611 {
 612     return ICommandLineOptionsModulePointer(new InsertMolecules());
 613 }
 614
 615 } // namespace gmx