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