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, 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 <string>
  43
  44 #include "gromacs/commandline/cmdlineoptionsmodule.h"
  45 #include "gromacs/fileio/confio.h"
  46 #include "gromacs/fileio/filenm.h"
  47 #include "gromacs/fileio/xvgr.h"
  48 #include "gromacs/gmxlib/conformation-utilities.h"
  49 #include "gromacs/gmxpreprocess/read-conformation.h"
  50 #include "gromacs/math/utilities.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/options/basicoptions.h"
  53 #include "gromacs/options/filenameoption.h"
  54 #include "gromacs/options/ioptionscontainer.h"
  55 #include "gromacs/pbcutil/pbc.h"
  56 #include "gromacs/random/random.h"
  57 #include "gromacs/selection/nbsearch.h"
  58 #include "gromacs/topology/atomprop.h"
  59 #include "gromacs/topology/atoms.h"
  60 #include "gromacs/utility/cstringutil.h"
  61 #include "gromacs/utility/fatalerror.h"
  62 #include "gromacs/utility/smalloc.h"
  63
  64 /* enum for random rotations of inserted solutes */
  65 enum {
  66     en_rotXYZ, en_rotZ, en_rotNone
  67 };
  68
  69 static void center_molecule(int atomCount, rvec x[])
  70 {
  71     rvec center;
  72     clear_rvec(center);
  73     for (int i = 0; i < atomCount; ++i)
  74     {
  75         rvec_inc(center, x[i]);
  76     }
  77     svmul(1.0/atomCount, center, center);
  78     for (int i = 0; i < atomCount; ++i)
  79     {
  80         rvec_dec(x[i], center);
  81     }
  82 }
  83
  84 static void generate_trial_conf(int atomCount, const rvec xin[],
  85                                 const rvec offset, int enum_rot, gmx_rng_t rng,
  86                                 rvec xout[])
  87 {
  88     for (int i = 0; i < atomCount; ++i)
  89     {
  90         copy_rvec(xin[i], xout[i]);
  91     }
  92     real alfa = 0.0, beta = 0.0, gamma = 0.0;
  93     switch (enum_rot)
  94     {
  95         case en_rotXYZ:
  96             alfa  = 2*M_PI * gmx_rng_uniform_real(rng);
  97             beta  = 2*M_PI * gmx_rng_uniform_real(rng);
  98             gamma = 2*M_PI * gmx_rng_uniform_real(rng);
  99             break;
 100         case en_rotZ:
 101             alfa  = beta = 0.;
 102             gamma = 2*M_PI * gmx_rng_uniform_real(rng);
 103             break;
 104         case en_rotNone:
 105             alfa = beta = gamma = 0.;
 106             break;
 107     }
 108     if (enum_rot == en_rotXYZ || (enum_rot == en_rotZ))
 109     {
 110         rotate_conf(atomCount, xout, NULL, alfa, beta, gamma);
 111     }
 112     for (int i = 0; i < atomCount; ++i)
 113     {
 114         rvec_inc(xout[i], offset);
 115     }
 116 }
 117
 118 static bool is_insertion_allowed(gmx::AnalysisNeighborhoodSearch *search,
 119                                  const real *exclusionDistances,
 120                                  int atomCount, const rvec *x,
 121                                  const real *exclusionDistances_insrt)
 122 {
 123     gmx::AnalysisNeighborhoodPositions  pos(x, atomCount);
 124     gmx::AnalysisNeighborhoodPairSearch pairSearch = search->startPairSearch(pos);
 125     gmx::AnalysisNeighborhoodPair       pair;
 126     while (pairSearch.findNextPair(&pair))
 127     {
 128         const real r1 = exclusionDistances[pair.refIndex()];
 129         const real r2 = exclusionDistances_insrt[pair.testIndex()];
 130         if (pair.distance2() < sqr(r1 + r2))
 131         {
 132             return false;
 133         }
 134     }
 135     return true;
 136 }
 137
 138 static void merge_atoms_noalloc(t_atoms *atoms, const t_atoms *atoms_add)
 139 {
 140     int resnr = 0;
 141     if (atoms->nr > 0)
 142     {
 143         resnr = atoms->resinfo[atoms->atom[atoms->nr-1].resind].nr;
 144     }
 145     int prevResInd = -1;
 146     for (int i = 0; i < atoms_add->nr; ++i)
 147     {
 148         if (atoms_add->atom[i].resind != prevResInd)
 149         {
 150             prevResInd = atoms_add->atom[i].resind;
 151             ++resnr;
 152             atoms->resinfo[atoms->nres]    = atoms_add->resinfo[prevResInd];
 153             atoms->resinfo[atoms->nres].nr = resnr;
 154             ++atoms->nres;
 155         }
 156         atoms->atom[atoms->nr]        = atoms_add->atom[i];
 157         atoms->atomname[atoms->nr]    = atoms_add->atomname[i];
 158         atoms->atom[atoms->nr].resind = atoms->nres-1;
 159         ++atoms->nr;
 160     }
 161 }
 162
 163 static void insert_mols(int nmol_insrt, int ntry, int seed,
 164                         real defaultDistance, real scaleFactor,
 165                         t_atoms *atoms, rvec **x,
 166                         const t_atoms *atoms_insrt, const rvec *x_insrt,
 167                         int ePBC, matrix box,
 168                         const std::string &posfn, const rvec deltaR, int enum_rot)
 169 {
 170     t_pbc            pbc;
 171     rvec            *x_n;
 172
 173     fprintf(stderr, "Initialising inter-atomic distances...\n");
 174     gmx_atomprop_t   aps = gmx_atomprop_init();
 175     real            *exclusionDistances
 176         = makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor);
 177     real            *exclusionDistances_insrt
 178         = makeExclusionDistances(atoms_insrt, aps, defaultDistance, scaleFactor);
 179     gmx_atomprop_destroy(aps);
 180
 181     const real       maxInsertRadius
 182         = *std::max_element(exclusionDistances_insrt,
 183                             exclusionDistances_insrt + atoms_insrt->nr);
 184     real             maxRadius = maxInsertRadius;
 185     if (atoms->nr > 0)
 186     {
 187         const real maxExistingRadius
 188             = *std::max_element(exclusionDistances,
 189                                 exclusionDistances + atoms->nr);
 190         maxRadius = std::max(maxInsertRadius, maxExistingRadius);
 191     }
 192
 193     // TODO: Make all of this exception-safe.
 194     gmx::AnalysisNeighborhood nb;
 195     nb.setCutoff(maxInsertRadius + maxRadius);
 196
 197     gmx_rng_t        rng = gmx_rng_init(seed);
 198     set_pbc(&pbc, ePBC, box);
 199
 200     snew(x_n, atoms_insrt->nr);
 201
 202     /* With -ip, take nmol_insrt from file posfn */
 203     double         **rpos = NULL;
 204     if (!posfn.empty())
 205     {
 206         int ncol;
 207         nmol_insrt = read_xvg(posfn.c_str(), &rpos, &ncol);
 208         if (ncol != 3)
 209         {
 210             gmx_fatal(FARGS, "Expected 3 columns (x/y/z coordinates) in file %s\n",
 211                       posfn.c_str());
 212         }
 213         fprintf(stderr, "Read %d positions from file %s\n\n",
 214                 nmol_insrt, posfn.c_str());
 215     }
 216
 217     {
 218         const int finalAtomCount    = atoms->nr + nmol_insrt * atoms_insrt->nr;
 219         const int finalResidueCount = atoms->nres + nmol_insrt * atoms_insrt->nres;
 220         srenew(atoms->resinfo,      finalResidueCount);
 221         srenew(atoms->atomname,     finalAtomCount);
 222         srenew(atoms->atom,         finalAtomCount);
 223         srenew(*x,                  finalAtomCount);
 224         srenew(exclusionDistances,  finalAtomCount);
 225     }
 226
 227     int mol        = 0;
 228     int trial      = 0;
 229     int firstTrial = 0;
 230     int failed     = 0;
 231     while ((mol < nmol_insrt) && (trial < ntry*nmol_insrt))
 232     {
 233         rvec offset_x;
 234         if (posfn.empty())
 235         {
 236             // Insert at random positions.
 237             offset_x[XX] = box[XX][XX] * gmx_rng_uniform_real(rng);
 238             offset_x[YY] = box[YY][YY] * gmx_rng_uniform_real(rng);
 239             offset_x[ZZ] = box[ZZ][ZZ] * gmx_rng_uniform_real(rng);
 240         }
 241         else
 242         {
 243             // Skip a position if ntry trials were not successful.
 244             if (trial >= firstTrial + ntry)
 245             {
 246                 fprintf(stderr, " skipped position (%.3f, %.3f, %.3f)\n",
 247                         rpos[XX][mol], rpos[YY][mol], rpos[ZZ][mol]);
 248                 ++mol;
 249                 ++failed;
 250             }
 251             // Insert at positions taken from option -ip file.
 252             offset_x[XX] = rpos[XX][mol] + deltaR[XX]*(2 * gmx_rng_uniform_real(rng)-1);
 253             offset_x[YY] = rpos[YY][mol] + deltaR[YY]*(2 * gmx_rng_uniform_real(rng)-1);
 254             offset_x[ZZ] = rpos[ZZ][mol] + deltaR[ZZ]*(2 * gmx_rng_uniform_real(rng)-1);
 255         }
 256         fprintf(stderr, "\rTry %d", ++trial);
 257         generate_trial_conf(atoms_insrt->nr, x_insrt, offset_x, enum_rot, rng,
 258                             x_n);
 259         gmx::AnalysisNeighborhoodPositions pos(*x, atoms->nr);
 260         gmx::AnalysisNeighborhoodSearch    search = nb.initSearch(&pbc, pos);
 261         if (is_insertion_allowed(&search, exclusionDistances, atoms_insrt->nr,
 262                                  x_n, exclusionDistances_insrt))
 263         {
 264             const int firstIndex = atoms->nr;
 265             for (int i = 0; i < atoms_insrt->nr; ++i)
 266             {
 267                 copy_rvec(x_n[i], (*x)[firstIndex + i]);
 268                 exclusionDistances[firstIndex + i] = exclusionDistances_insrt[i];
 269             }
 270             merge_atoms_noalloc(atoms, atoms_insrt);
 271             ++mol;
 272             firstTrial = trial;
 273             fprintf(stderr, " success (now %d atoms)!\n", atoms->nr);
 274         }
 275     }
 276     gmx_rng_destroy(rng);
 277     srenew(atoms->resinfo,  atoms->nres);
 278     srenew(atoms->atomname, atoms->nr);
 279     srenew(atoms->atom,     atoms->nr);
 280     srenew(*x,              atoms->nr);
 281
 282     fprintf(stderr, "\n");
 283     /* print number of molecules added */
 284     fprintf(stderr, "Added %d molecules (out of %d requested)\n",
 285             mol - failed, nmol_insrt);
 286
 287     sfree(x_n);
 288     sfree(exclusionDistances);
 289     sfree(exclusionDistances_insrt);
 290     if (rpos != NULL)
 291     {
 292         for (int i = 0; i < DIM; ++i)
 293         {
 294             sfree(rpos[i]);
 295         }
 296         sfree(rpos);
 297     }
 298 }
 299
 300 namespace gmx
 301 {
 302
 303 namespace
 304 {
 305
 306 class InsertMolecules : public ICommandLineOptionsModule
 307 {
 308     public:
 309         InsertMolecules()
 310             : bBox_(false), nmolIns_(0), nmolTry_(10), seed_(1997),
 311               defaultDistance_(0.105), scaleFactor_(0.57), enumRot_(en_rotXYZ)
 312         {
 313             clear_rvec(newBox_);
 314             clear_rvec(deltaR_);
 315         }
 316
 317         virtual void init(CommandLineModuleSettings * /*settings*/)
 318         {
 319         }
 320
 321         virtual void initOptions(IOptionsContainer                 *options,
 322                                  ICommandLineOptionsModuleSettings *settings);
 323         virtual void optionsFinished() {}
 324
 325         virtual int run();
 326
 327     private:
 328         std::string inputConfFile_;
 329         std::string insertConfFile_;
 330         std::string positionFile_;
 331         std::string outputConfFile_;
 332         rvec        newBox_;
 333         bool        bBox_;
 334         int         nmolIns_;
 335         int         nmolTry_;
 336         int         seed_;
 337         real        defaultDistance_;
 338         real        scaleFactor_;
 339         rvec        deltaR_;
 340         int         enumRot_;
 341 };
 342
 343 void InsertMolecules::initOptions(IOptionsContainer                 *options,
 344                                   ICommandLineOptionsModuleSettings *settings)
 345 {
 346     const char *const desc[] = {
 347         "[THISMODULE] inserts [TT]-nmol[tt] copies of the system specified in",
 348         "the [TT]-ci[tt] input file. The insertions take place either into",
 349         "vacant space in the solute conformation given with [TT]-f[tt], or",
 350         "into an empty box given by [TT]-box[tt]. Specifying both [TT]-f[tt]",
 351         "and [TT]-box[tt] behaves like [TT]-f[tt], but places a new box",
 352         "around the solute before insertions. Any velocities present are",
 353         "discarded.[PAR]",
 354
 355         "By default, the insertion positions are random (with initial seed",
 356         "specified by [TT]-seed[tt]). The program iterates until [TT]-nmol[tt]",
 357         "molecules have been inserted in the box. Molecules are not inserted",
 358         "where the distance between any existing atom and any atom of the",
 359         "inserted molecule is less than the sum based on the van der Waals",
 360         "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
 361         "Waals radii is read by the program, and the resulting radii scaled",
 362         "by [TT]-scale[tt]. If radii are not found in the database, those"
 363         "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].[PAR]",
 364
 365         "A total of [TT]-nmol[tt] * [TT]-try[tt] insertion attempts are made",
 366         "before giving up. Increase [TT]-try[tt] if you have several small",
 367         "holes to fill. Option [TT]-rot[tt] specifies whether the insertion",
 368         "molecules are randomly oriented before insertion attempts.[PAR]",
 369
 370         "Alternatively, the molecules can be inserted only at positions defined in",
 371         "positions.dat ([TT]-ip[tt]). That file should have 3 columns (x,y,z),",
 372         "that give the displacements compared to the input molecule position",
 373         "([TT]-ci[tt]). Hence, if that file should contain the absolute",
 374         "positions, the molecule must be centered on (0,0,0) before using",
 375         "[THISMODULE] (e.g. from [gmx-editconf] [TT]-center[tt]).",
 376         "Comments in that file starting with # are ignored. Option [TT]-dr[tt]",
 377         "defines the maximally allowed displacements during insertial trials.",
 378         "[TT]-try[tt] and [TT]-rot[tt] work as in the default mode (see above)."
 379     };
 380
 381     settings->setHelpText(desc);
 382
 383     // TODO: Replace use of legacyType.
 384     options->addOption(FileNameOption("f")
 385                            .legacyType(efSTX).inputFile()
 386                            .store(&inputConfFile_)
 387                            .defaultBasename("protein")
 388                            .description("Existing configuration to insert into"));
 389     options->addOption(FileNameOption("ci")
 390                            .legacyType(efSTX).inputFile().required()
 391                            .store(&insertConfFile_)
 392                            .defaultBasename("insert")
 393                            .description("Configuration to insert"));
 394     options->addOption(FileNameOption("ip")
 395                            .filetype(eftGenericData).inputFile()
 396                            .store(&positionFile_)
 397                            .defaultBasename("positions")
 398                            .description("Predefined insertion trial positions"));
 399     options->addOption(FileNameOption("o")
 400                            .legacyType(efSTO).outputFile().required()
 401                            .store(&outputConfFile_)
 402                            .defaultBasename("out")
 403                            .description("Output configuration after insertion"));
 404
 405     options->addOption(RealOption("box").vector()
 406                            .store(newBox_).storeIsSet(&bBox_)
 407                            .description("Box size (in nm)"));
 408     options->addOption(IntegerOption("nmol")
 409                            .store(&nmolIns_)
 410                            .description("Number of extra molecules to insert"));
 411     options->addOption(IntegerOption("try")
 412                            .store(&nmolTry_)
 413                            .description("Try inserting [TT]-nmol[tt] times [TT]-try[tt] times"));
 414     options->addOption(IntegerOption("seed")
 415                            .store(&seed_)
 416                            .description("Random generator seed"));
 417     options->addOption(RealOption("radius")
 418                            .store(&defaultDistance_)
 419                            .description("Default van der Waals distance"));
 420     options->addOption(RealOption("scale")
 421                            .store(&scaleFactor_)
 422                            .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."));
 423     options->addOption(RealOption("dr").vector()
 424                            .store(deltaR_)
 425                            .description("Allowed displacement in x/y/z from positions in [TT]-ip[tt] file"));
 426     const char *const cRotationEnum[] = {"xyz", "z", "none"};
 427     options->addOption(StringOption("rot").enumValue(cRotationEnum)
 428                            .storeEnumIndex(&enumRot_)
 429                            .description("Rotate inserted molecules randomly"));
 430 }
 431
 432 int InsertMolecules::run()
 433 {
 434     const bool bProt = !inputConfFile_.empty();
 435
 436     /* check input */
 437     if (nmolIns_ <= 0 && positionFile_.empty())
 438     {
 439         gmx_fatal(FARGS, "Either -nmol must be larger than 0, "
 440                   "or positions must be given with -ip");
 441     }
 442     if (!bProt && !bBox_)
 443     {
 444         gmx_fatal(FARGS, "When no solute (-f) is specified, "
 445                   "a box size (-box) must be specified");
 446     }
 447
 448     char    *title = NULL;
 449     t_atoms *atoms;
 450     rvec    *x = NULL;
 451     matrix   box;
 452     int      ePBC = -1;
 453     snew(atoms, 1);
 454     init_t_atoms(atoms, 0, FALSE);
 455     if (bProt)
 456     {
 457         /* Generate a solute configuration */
 458         title = readConformation(inputConfFile_.c_str(), atoms, &x, NULL,
 459                                  &ePBC, box, "solute");
 460         if (atoms->nr == 0)
 461         {
 462             fprintf(stderr, "Note: no atoms in %s\n", inputConfFile_.c_str());
 463             sfree(title);
 464             title = NULL;
 465         }
 466     }
 467     if (bBox_)
 468     {
 469         ePBC = epbcXYZ;
 470         clear_mat(box);
 471         box[XX][XX] = newBox_[XX];
 472         box[YY][YY] = newBox_[YY];
 473         box[ZZ][ZZ] = newBox_[ZZ];
 474     }
 475     if (det(box) == 0)
 476     {
 477         gmx_fatal(FARGS, "Undefined solute box.\nCreate one with gmx editconf "
 478                   "or give explicit -box command line option");
 479     }
 480
 481     t_atoms *atoms_insrt;
 482     rvec    *x_insrt = NULL;
 483     snew(atoms_insrt, 1);
 484     init_t_atoms(atoms_insrt, 0, FALSE);
 485     {
 486         int         ePBC_dummy;
 487         matrix      box_dummy;
 488         char       *title_ins
 489             = readConformation(insertConfFile_.c_str(), atoms_insrt, &x_insrt,
 490                                NULL, &ePBC_dummy, box_dummy, "molecule");
 491         if (atoms_insrt->nr == 0)
 492         {
 493             gmx_fatal(FARGS, "No molecule in %s, please check your input",
 494                       insertConfFile_.c_str());
 495         }
 496         if (title == NULL)
 497         {
 498             title = title_ins;
 499         }
 500         else
 501         {
 502             sfree(title_ins);
 503         }
 504         if (positionFile_.empty())
 505         {
 506             center_molecule(atoms_insrt->nr, x_insrt);
 507         }
 508     }
 509
 510     /* add nmol_ins molecules of atoms_ins
 511        in random orientation at random place */
 512     insert_mols(nmolIns_, nmolTry_, seed_, defaultDistance_, scaleFactor_,
 513                 atoms, &x, atoms_insrt, x_insrt,
 514                 ePBC, box, positionFile_, deltaR_, enumRot_);
 515
 516     /* write new configuration to file confout */
 517     fprintf(stderr, "Writing generated configuration to %s\n",
 518             outputConfFile_.c_str());
 519     write_sto_conf(outputConfFile_.c_str(), title, atoms, x, NULL, ePBC, box);
 520
 521     /* print size of generated configuration */
 522     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
 523             atoms->nr, atoms->nres);
 524
 525     sfree(x);
 526     sfree(x_insrt);
 527     done_atom(atoms);
 528     done_atom(atoms_insrt);
 529     sfree(atoms);
 530     sfree(atoms_insrt);
 531     sfree(title);
 532
 533     return 0;
 534 }
 535
 536 }   // namespace
 537
 538 const char InsertMoleculesInfo::name[]             = "insert-molecules";
 539 const char InsertMoleculesInfo::shortDescription[] =
 540     "Insert molecules into existing vacancies";
 541 ICommandLineOptionsModule *InsertMoleculesInfo::create()
 542 {
 543     return new InsertMolecules();
 544 }
 545
 546 } // namespace gmx