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