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,2017,2018, 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()), nullptr, 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              = nullptr;
 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 != nullptr)
 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               ePBC_(-1)
 333         {
 334             clear_rvec(newBox_);
 335             clear_rvec(deltaR_);
 336             clear_mat(box_);
 337         }
 338
 339         // From ITopologyProvider
 340         virtual gmx_mtop_t *getTopology(bool /*required*/) { return &top_; }
 341         virtual int getAtomCount() { return 0; }
 342
 343         // From ICommandLineOptionsModule
 344         virtual void init(CommandLineModuleSettings * /*settings*/)
 345         {
 346         }
 347         virtual void initOptions(IOptionsContainer                 *options,
 348                                  ICommandLineOptionsModuleSettings *settings);
 349         virtual void optionsFinished();
 350         virtual int run();
 351
 352     private:
 353         SelectionCollection selections_;
 354
 355         std::string         inputConfFile_;
 356         std::string         insertConfFile_;
 357         std::string         positionFile_;
 358         std::string         outputConfFile_;
 359         rvec                newBox_;
 360         bool                bBox_;
 361         int                 nmolIns_;
 362         int                 nmolTry_;
 363         int                 seed_;
 364         real                defaultDistance_;
 365         real                scaleFactor_;
 366         rvec                deltaR_;
 367         RotationType        enumRot_;
 368         Selection           replaceSel_;
 369
 370         gmx_mtop_t          top_;
 371         std::vector<RVec>   x_;
 372         matrix              box_;
 373         int                 ePBC_;
 374 };
 375
 376 void InsertMolecules::initOptions(IOptionsContainer                 *options,
 377                                   ICommandLineOptionsModuleSettings *settings)
 378 {
 379     const char *const desc[] = {
 380         "[THISMODULE] inserts [TT]-nmol[tt] copies of the system specified in",
 381         "the [TT]-ci[tt] input file. The insertions take place either into",
 382         "vacant space in the solute conformation given with [TT]-f[tt], or",
 383         "into an empty box given by [TT]-box[tt]. Specifying both [TT]-f[tt]",
 384         "and [TT]-box[tt] behaves like [TT]-f[tt], but places a new box",
 385         "around the solute before insertions. Any velocities present are",
 386         "discarded.",
 387         "",
 388         "It is possible to also insert into a solvated configuration and",
 389         "replace solvent atoms with the inserted atoms. To do this, use",
 390         "[TT]-replace[tt] to specify a selection that identifies the atoms",
 391         "that can be replaced. The tool assumes that all molecules in this",
 392         "selection consist of single residues: each residue from this",
 393         "selection that overlaps with the inserted molecules will be removed",
 394         "instead of preventing insertion.",
 395         "",
 396         "By default, the insertion positions are random (with initial seed",
 397         "specified by [TT]-seed[tt]). The program iterates until [TT]-nmol[tt]",
 398         "molecules have been inserted in the box. Molecules are not inserted",
 399         "where the distance between any existing atom and any atom of the",
 400         "inserted molecule is less than the sum based on the van der Waals",
 401         "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
 402         "Waals radii is read by the program, and the resulting radii scaled",
 403         "by [TT]-scale[tt]. If radii are not found in the database, those",
 404         "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].",
 405         "Note that the usefulness of those radii depends on the atom names,",
 406         "and thus varies widely with force field.",
 407         "",
 408         "A total of [TT]-nmol[tt] * [TT]-try[tt] insertion attempts are made",
 409         "before giving up. Increase [TT]-try[tt] if you have several small",
 410         "holes to fill. Option [TT]-rot[tt] specifies whether the insertion",
 411         "molecules are randomly oriented before insertion attempts.",
 412         "",
 413         "Alternatively, the molecules can be inserted only at positions defined in",
 414         "positions.dat ([TT]-ip[tt]). That file should have 3 columns (x,y,z),",
 415         "that give the displacements compared to the input molecule position",
 416         "([TT]-ci[tt]). Hence, if that file should contain the absolute",
 417         "positions, the molecule must be centered on (0,0,0) before using",
 418         "[THISMODULE] (e.g. from [gmx-editconf] [TT]-center[tt]).",
 419         "Comments in that file starting with # are ignored. Option [TT]-dr[tt]",
 420         "defines the maximally allowed displacements during insertial trials.",
 421         "[TT]-try[tt] and [TT]-rot[tt] work as in the default mode (see above)."
 422     };
 423
 424     settings->setHelpText(desc);
 425
 426     std::shared_ptr<SelectionOptionBehavior> selectionOptionBehavior(
 427             new SelectionOptionBehavior(&selections_, this));
 428     settings->addOptionsBehavior(selectionOptionBehavior);
 429
 430     // TODO: Replace use of legacyType.
 431     options->addOption(FileNameOption("f")
 432                            .legacyType(efSTX).inputFile()
 433                            .store(&inputConfFile_)
 434                            .defaultBasename("protein")
 435                            .description("Existing configuration to insert into"));
 436     options->addOption(FileNameOption("ci")
 437                            .legacyType(efSTX).inputFile().required()
 438                            .store(&insertConfFile_)
 439                            .defaultBasename("insert")
 440                            .description("Configuration to insert"));
 441     options->addOption(FileNameOption("ip")
 442                            .filetype(eftGenericData).inputFile()
 443                            .store(&positionFile_)
 444                            .defaultBasename("positions")
 445                            .description("Predefined insertion trial positions"));
 446     options->addOption(FileNameOption("o")
 447                            .legacyType(efSTO).outputFile().required()
 448                            .store(&outputConfFile_)
 449                            .defaultBasename("out")
 450                            .description("Output configuration after insertion"));
 451
 452     options->addOption(SelectionOption("replace").onlyAtoms()
 453                            .store(&replaceSel_)
 454                            .description("Atoms that can be removed if overlapping"));
 455     selectionOptionBehavior->initOptions(options);
 456
 457     options->addOption(RealOption("box").vector()
 458                            .store(newBox_).storeIsSet(&bBox_)
 459                            .description("Box size (in nm)"));
 460     options->addOption(IntegerOption("nmol")
 461                            .store(&nmolIns_)
 462                            .description("Number of extra molecules to insert"));
 463     options->addOption(IntegerOption("try")
 464                            .store(&nmolTry_)
 465                            .description("Try inserting [TT]-nmol[tt] times [TT]-try[tt] times"));
 466     options->addOption(IntegerOption("seed")
 467                            .store(&seed_)
 468                            .description("Random generator seed (0 means generate)"));
 469     options->addOption(RealOption("radius")
 470                            .store(&defaultDistance_)
 471                            .description("Default van der Waals distance"));
 472     options->addOption(RealOption("scale")
 473                            .store(&scaleFactor_)
 474                            .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."));
 475     options->addOption(RealOption("dr").vector()
 476                            .store(deltaR_)
 477                            .description("Allowed displacement in x/y/z from positions in [TT]-ip[tt] file"));
 478     options->addOption(EnumOption<RotationType>("rot").enumValue(cRotationEnum)
 479                            .store(&enumRot_)
 480                            .description("Rotate inserted molecules randomly"));
 481 }
 482
 483 void InsertMolecules::optionsFinished()
 484 {
 485     if (nmolIns_ <= 0 && positionFile_.empty())
 486     {
 487         GMX_THROW(InconsistentInputError("Either -nmol must be larger than 0, "
 488                                          "or positions must be given with -ip."));
 489     }
 490     if (inputConfFile_.empty() && !bBox_)
 491     {
 492         GMX_THROW(InconsistentInputError("When no solute (-f) is specified, "
 493                                          "a box size (-box) must be specified."));
 494     }
 495     if (replaceSel_.isValid() && inputConfFile_.empty())
 496     {
 497         GMX_THROW(InconsistentInputError("Replacement (-replace) only makes sense "
 498                                          "together with an existing configuration (-f)."));
 499     }
 500
 501     if (!inputConfFile_.empty())
 502     {
 503         readConformation(inputConfFile_.c_str(), &top_, &x_, nullptr,
 504                          &ePBC_, box_, "solute");
 505         if (top_.natoms == 0)
 506         {
 507             fprintf(stderr, "Note: no atoms in %s\n", inputConfFile_.c_str());
 508         }
 509     }
 510 }
 511
 512 int InsertMolecules::run()
 513 {
 514     std::set<int> removableAtoms;
 515     if (replaceSel_.isValid())
 516     {
 517         t_pbc       pbc;
 518         set_pbc(&pbc, ePBC_, box_);
 519         t_trxframe *fr;
 520         snew(fr, 1);
 521         fr->natoms = x_.size();
 522         fr->bX     = TRUE;
 523         fr->x      = as_rvec_array(x_.data());
 524         selections_.evaluate(fr, &pbc);
 525         sfree(fr);
 526         removableAtoms.insert(replaceSel_.atomIndices().begin(),
 527                               replaceSel_.atomIndices().end());
 528         // TODO: It could be nice to check that removableAtoms contains full
 529         // residues, since we anyways remove whole residues instead of
 530         // individual atoms.
 531     }
 532
 533     if (bBox_)
 534     {
 535         ePBC_ = epbcXYZ;
 536         clear_mat(box_);
 537         box_[XX][XX] = newBox_[XX];
 538         box_[YY][YY] = newBox_[YY];
 539         box_[ZZ][ZZ] = newBox_[ZZ];
 540     }
 541     if (det(box_) == 0)
 542     {
 543         gmx_fatal(FARGS, "Undefined solute box.\nCreate one with gmx editconf "
 544                   "or give explicit -box command line option");
 545     }
 546
 547     t_topology        *top_insrt;
 548     std::vector<RVec>  x_insrt;
 549     snew(top_insrt, 1);
 550     {
 551         int         ePBC_dummy;
 552         matrix      box_dummy;
 553         readConformation(insertConfFile_.c_str(), top_insrt, &x_insrt,
 554                          nullptr, &ePBC_dummy, box_dummy, "molecule");
 555         if (top_insrt->atoms.nr == 0)
 556         {
 557             gmx_fatal(FARGS, "No molecule in %s, please check your input",
 558                       insertConfFile_.c_str());
 559         }
 560         if (top_.name == nullptr)
 561         {
 562             top_.name = top_insrt->name;
 563         }
 564         if (positionFile_.empty())
 565         {
 566             center_molecule(&x_insrt);
 567         }
 568     }
 569
 570     // TODO: Adapt to use mtop throughout.
 571     t_atoms atoms = gmx_mtop_global_atoms(&top_);
 572
 573     /* add nmol_ins molecules of atoms_ins
 574        in random orientation at random place */
 575     insert_mols(nmolIns_, nmolTry_, seed_, defaultDistance_, scaleFactor_,
 576                 &atoms, &top_.symtab, &x_, removableAtoms, top_insrt->atoms, x_insrt,
 577                 ePBC_, box_, positionFile_, deltaR_, enumRot_);
 578
 579     /* write new configuration to file confout */
 580     fprintf(stderr, "Writing generated configuration to %s\n",
 581             outputConfFile_.c_str());
 582     write_sto_conf(outputConfFile_.c_str(), *top_.name, &atoms,
 583                    as_rvec_array(x_.data()), nullptr, ePBC_, box_);
 584
 585     /* print size of generated configuration */
 586     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
 587             atoms.nr, atoms.nres);
 588
 589     done_atom(&atoms);
 590     done_top(top_insrt);
 591     sfree(top_insrt);
 592
 593     return 0;
 594 }
 595
 596 }   // namespace
 597
 598 const char InsertMoleculesInfo::name[]             = "insert-molecules";
 599 const char InsertMoleculesInfo::shortDescription[] =
 600     "Insert molecules into existing vacancies";
 601 ICommandLineOptionsModulePointer InsertMoleculesInfo::create()
 602 {
 603     return ICommandLineOptionsModulePointer(new InsertMolecules());
 604 }
 605
 606 } // namespace gmx