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, 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 "gromacs/commandline/pargs.h"
  42 #include "gromacs/fileio/confio.h"
  43 #include "gromacs/fileio/xvgr.h"
  44 #include "gromacs/gmxlib/conformation-utilities.h"
  45 #include "gromacs/gmxpreprocess/addconf.h"
  46 #include "gromacs/gmxpreprocess/read-conformation.h"
  47 #include "gromacs/legacyheaders/macros.h"
  48 #include "gromacs/legacyheaders/names.h"
  49 #include "gromacs/legacyheaders/typedefs.h"
  50 #include "gromacs/math/utilities.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/pbcutil/ishift.h"
  53 #include "gromacs/pbcutil/pbc.h"
  54 #include "gromacs/random/random.h"
  55 #include "gromacs/topology/atomprop.h"
  56 #include "gromacs/utility/cstringutil.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/futil.h"
  59 #include "gromacs/utility/smalloc.h"
  60
  61 static gmx_bool in_box(t_pbc *pbc, rvec x)
  62 {
  63     rvec box_center, dx;
  64     int  shift;
  65
  66     /* pbc_dx_aiuc only works correctly with the rectangular box center */
  67     calc_box_center(ecenterRECT, pbc->box, box_center);
  68
  69     shift = pbc_dx_aiuc(pbc, x, box_center, dx);
  70
  71     return (shift == CENTRAL);
  72 }
  73
  74 /* This is a (maybe) slow workaround to avoid the neighbor searching in addconf.c, which
  75  * leaks memory (May 2012). The function could be deleted as soon as the memory leaks
  76  * there are fixed.
  77  * However, when inserting a small molecule in a system containing not too many atoms,
  78  * allPairsDistOk is probably even faster than the other code.
  79  */
  80 static gmx_bool
  81 allPairsDistOk(t_atoms *atoms, rvec *x, real *exclusionDistances,
  82                int ePBC, matrix box,
  83                t_atoms *atoms_insrt, rvec *x_n, real *exclusionDistances_insrt)
  84 {
  85     int   i, j;
  86     rvec  dx;
  87     real  n2, r2;
  88     t_pbc pbc;
  89
  90     set_pbc(&pbc, ePBC, box);
  91     for (i = 0; i < atoms->nr; i++)
  92     {
  93         for (j = 0; j < atoms_insrt->nr; j++)
  94         {
  95             pbc_dx(&pbc, x[i], x_n[j], dx);
  96             n2 = norm2(dx);
  97             r2 = sqr(exclusionDistances[i]+exclusionDistances_insrt[j]);
  98             if (n2 < r2)
  99             {
 100                 return FALSE;
 101             }
 102         }
 103     }
 104     return TRUE;
 105 }
 106
 107 /* enum for random rotations of inserted solutes */
 108 enum {
 109     en_rot, en_rotXYZ, en_rotZ, en_rotNone, en_NR
 110 };
 111
 112 static char *insert_mols(const char *mol_insrt, int nmol_insrt, int ntry, int seed,
 113                          t_atoms *atoms, rvec **x, real **exclusionDistances, int ePBC, matrix box,
 114                          gmx_atomprop_t aps,
 115                          real defaultDistance, real scaleFactor, real rshell,
 116                          const output_env_t oenv,
 117                          const char* posfn, const rvec deltaR, int enum_rot,
 118                          gmx_bool bCheckAllPairDist)
 119 {
 120     t_pbc            pbc;
 121     static  char    *title_insrt;
 122     t_atoms          atoms_insrt;
 123     rvec            *x_insrt, *x_n;
 124     real            *exclusionDistances_insrt;
 125     int              ePBC_insrt;
 126     matrix           box_insrt;
 127     int              i, mol, onr, ncol;
 128     real             alfa = 0., beta = 0., gamma = 0.;
 129     rvec             offset_x;
 130     int              trial;
 131     double         **rpos;
 132     gmx_rng_t        rng;
 133
 134     rng = gmx_rng_init(seed);
 135     set_pbc(&pbc, ePBC, box);
 136
 137     /* read number of atoms of insert molecules */
 138     {
 139         int natoms;
 140         get_stx_coordnum(mol_insrt, &natoms);
 141         if (natoms == 0)
 142         {
 143             gmx_fatal(FARGS, "No molecule in %s, please check your input\n", mol_insrt);
 144         }
 145         init_t_atoms(&atoms_insrt, natoms, FALSE);
 146     }
 147     /* allocate memory for atom coordinates of insert molecules */
 148     snew(x_insrt, atoms_insrt.nr);
 149     snew(atoms_insrt.resinfo, atoms_insrt.nr);
 150     snew(atoms_insrt.atomname, atoms_insrt.nr);
 151     snew(atoms_insrt.atom, atoms_insrt.nr);
 152     atoms_insrt.pdbinfo = NULL;
 153     snew(x_n, atoms_insrt.nr);
 154     snew(title_insrt, STRLEN);
 155
 156     /* read residue number, residue names, atomnames, coordinates etc. */
 157     fprintf(stderr, "Reading molecule configuration \n");
 158     read_stx_conf(mol_insrt, title_insrt, &atoms_insrt, x_insrt, NULL,
 159                   &ePBC_insrt, box_insrt);
 160     fprintf(stderr, "%s\nContaining %d atoms in %d residue\n",
 161             title_insrt, atoms_insrt.nr, atoms_insrt.nres);
 162     srenew(atoms_insrt.resinfo, atoms_insrt.nres);
 163
 164     /* initialise distance arrays for inserted molecules */
 165     exclusionDistances_insrt = makeExclusionDistances(&atoms_insrt, aps, defaultDistance, scaleFactor);
 166
 167     /* With -ip, take nmol_insrt from file posfn */
 168     if (posfn != NULL)
 169     {
 170         nmol_insrt = read_xvg(posfn, &rpos, &ncol);
 171         if (ncol != 3)
 172         {
 173             gmx_fatal(FARGS, "Expected 3 columns (x/y/z coordinates) in file %s\n", ncol, posfn);
 174         }
 175         fprintf(stderr, "Read %d positions from file %s\n\n", nmol_insrt, posfn);
 176     }
 177
 178     srenew(atoms->resinfo, (atoms->nres+nmol_insrt*atoms_insrt.nres));
 179     srenew(atoms->atomname, (atoms->nr+atoms_insrt.nr*nmol_insrt));
 180     srenew(atoms->atom, (atoms->nr+atoms_insrt.nr*nmol_insrt));
 181     srenew(*x, (atoms->nr+atoms_insrt.nr*nmol_insrt));
 182     srenew(*exclusionDistances, (atoms->nr+atoms_insrt.nr*nmol_insrt));
 183
 184     trial = mol = 0;
 185     while ((mol < nmol_insrt) && (trial < ntry*nmol_insrt))
 186     {
 187         fprintf(stderr, "\rTry %d", trial++);
 188         for (i = 0; (i < atoms_insrt.nr); i++)
 189         {
 190             copy_rvec(x_insrt[i], x_n[i]);
 191         }
 192         switch (enum_rot)
 193         {
 194             case en_rotXYZ:
 195                 alfa  = 2*M_PI * gmx_rng_uniform_real(rng);
 196                 beta  = 2*M_PI * gmx_rng_uniform_real(rng);
 197                 gamma = 2*M_PI * gmx_rng_uniform_real(rng);
 198                 break;
 199             case en_rotZ:
 200                 alfa  = beta = 0.;
 201                 gamma = 2*M_PI * gmx_rng_uniform_real(rng);
 202                 break;
 203             case en_rotNone:
 204                 alfa = beta = gamma = 0.;
 205                 break;
 206         }
 207         if (enum_rot == en_rotXYZ || (enum_rot == en_rotZ))
 208         {
 209             rotate_conf(atoms_insrt.nr, x_n, NULL, alfa, beta, gamma);
 210         }
 211         if (posfn == NULL)
 212         {
 213             /* insert at random positions */
 214             offset_x[XX] = box[XX][XX] * gmx_rng_uniform_real(rng);
 215             offset_x[YY] = box[YY][YY] * gmx_rng_uniform_real(rng);
 216             offset_x[ZZ] = box[ZZ][ZZ] * gmx_rng_uniform_real(rng);
 217             make_new_box(atoms_insrt.nr, x_n, box_insrt, offset_x, TRUE);
 218             if (!in_box(&pbc, x_n[0]) || !in_box(&pbc, x_n[atoms_insrt.nr-1]))
 219             {
 220                 continue;
 221             }
 222         }
 223         else
 224         {
 225             /* Insert at positions taken from option -ip file */
 226             offset_x[XX] = rpos[XX][mol] + deltaR[XX]*(2 * gmx_rng_uniform_real(rng)-1);
 227             offset_x[YY] = rpos[YY][mol] + deltaR[YY]*(2 * gmx_rng_uniform_real(rng)-1);
 228             offset_x[ZZ] = rpos[ZZ][mol] + deltaR[ZZ]*(2 * gmx_rng_uniform_real(rng)-1);
 229             for (i = 0; i < atoms_insrt.nr; i++)
 230             {
 231                 rvec_inc(x_n[i], offset_x);
 232             }
 233         }
 234
 235         onr = atoms->nr;
 236
 237         /* This is a (maybe) slow workaround to avoid too many calls of add_conf, which
 238          * leaks memory (status May 2012). If the momory leaks in add_conf() are fixed,
 239          * this check could be removed. Note, however, that allPairsDistOk is probably
 240          * even faster than add_conf() when inserting a small molecule into a moderately
 241          * small system.
 242          */
 243         if (bCheckAllPairDist && !allPairsDistOk(atoms, *x, *exclusionDistances, ePBC, box, &atoms_insrt, x_n, exclusionDistances_insrt))
 244         {
 245             continue;
 246         }
 247
 248         add_conf(atoms, x, NULL, exclusionDistances, FALSE, ePBC, box, TRUE,
 249                  &atoms_insrt, x_n, NULL, exclusionDistances_insrt, FALSE, rshell, 0, oenv);
 250
 251         if (atoms->nr == (atoms_insrt.nr+onr))
 252         {
 253             mol++;
 254             fprintf(stderr, " success (now %d atoms)!\n", atoms->nr);
 255         }
 256     }
 257     gmx_rng_destroy(rng);
 258     srenew(atoms->resinfo,  atoms->nres);
 259     srenew(atoms->atomname, atoms->nr);
 260     srenew(atoms->atom,     atoms->nr);
 261     srenew(*x,              atoms->nr);
 262     srenew(*exclusionDistances, atoms->nr);
 263
 264     fprintf(stderr, "\n");
 265     /* print number of molecules added */
 266     fprintf(stderr, "Added %d molecules (out of %d requested) of %s\n",
 267             mol, nmol_insrt, *atoms_insrt.resinfo[0].name);
 268
 269     sfree(x_n);
 270     sfree(exclusionDistances_insrt);
 271     done_atom(&atoms_insrt);
 272
 273     return title_insrt;
 274 }
 275
 276 int gmx_insert_molecules(int argc, char *argv[])
 277 {
 278     const char *desc[] = {
 279         "[THISMODULE] inserts [TT]-nmol[tt] copies of the system specified in",
 280         "the [TT]-ci[tt] input file. The insertions take place either into",
 281         "vacant space in the solute conformation given with [TT]-f[tt], or",
 282         "into an empty box given by [TT]-box[tt]. Specifying both [TT]-f[tt]",
 283         "and [TT]-box[tt] behaves like [TT]-f[tt], but places a new box",
 284         "around the solute before insertions. Any velocities present are",
 285         "discarded.[PAR]",
 286
 287         "By default, the insertion positions are random (with initial seed",
 288         "specified by [TT]-seed[tt]). The program iterates until [TT]-nmol[tt]",
 289         "molecules have been inserted in the box. Molecules are not inserted",
 290         "where the distance between any existing atom and any atom of the",
 291         "inserted molecule is less than the sum based on the van der Waals",
 292         "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
 293         "Waals radii is read by the program, and the resulting radii scaled",
 294         "by [TT]-scale[tt]. If radii are not found in the database, those"
 295         "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].[PAR]",
 296
 297         "A total of [TT]-nmol[tt] * [TT]-try[tt] insertion attempts are made",
 298         "before giving up. Increase [TT]-try[tt] if you have several small",
 299         "holes to fill. Option [TT]-rot[tt] specifies whether the insertion",
 300         "molecules are randomly oriented before insertion attempts.[PAR]",
 301
 302         "Alternatively, the molecules can be inserted only at positions defined in",
 303         "positions.dat ([TT]-ip[tt]). That file should have 3 columns (x,y,z),",
 304         "that give the displacements compared to the input molecule position",
 305         "([TT]-ci[tt]). Hence, if that file should contain the absolute",
 306         "positions, the molecule must be centered on (0,0,0) before using",
 307         "[THISMODULE] (e.g. from [gmx-editconf] [TT]-center[tt]).",
 308         "Comments in that file starting with # are ignored. Option [TT]-dr[tt]",
 309         "defines the maximally allowed displacements during insertial trials.",
 310         "[TT]-try[tt] and [TT]-rot[tt] work as in the default mode (see above).",
 311         "[PAR]",
 312     };
 313
 314     const char *bugs[] = {
 315         "Molecules must be whole in the initial configurations.",
 316         "Many repeated neighbor searchings with -ci blows up the allocated memory. "
 317         "Option -allpair avoids this using all-to-all distance checks (slow for large systems)"
 318     };
 319
 320     /* parameter data */
 321     gmx_bool       bProt, bBox;
 322     const char    *conf_prot, *confout;
 323     real          *exclusionDistances = NULL;
 324     char          *title_ins          = NULL;
 325     gmx_atomprop_t aps;
 326
 327     /* protein configuration data */
 328     char          *title = NULL;
 329     t_atoms       *atoms;
 330     rvec          *x    = NULL;
 331     int            ePBC = -1;
 332     matrix         box;
 333
 334     t_filenm       fnm[] = {
 335         { efSTX, "-f", "protein", ffOPTRD },
 336         { efSTX, "-ci", "insert",  ffREAD},
 337         { efDAT, "-ip", "positions",  ffOPTRD},
 338         { efSTO, NULL,  NULL,      ffWRITE},
 339     };
 340 #define NFILE asize(fnm)
 341
 342     static int      nmol_ins               = 0, nmol_try = 10, seed = 1997, enum_rot;
 343     static real     defaultDistance        = 0.105, scaleFactor = 0.57;
 344     static rvec     new_box                = {0.0, 0.0, 0.0}, deltaR = {0.0, 0.0, 0.0};
 345     static gmx_bool bCheckAllPairDist      = FALSE;
 346     output_env_t    oenv;
 347     const char     *enum_rot_string[] = {NULL, "xyz", "z", "none", NULL};
 348     t_pargs         pa[]              = {
 349         { "-box",    FALSE, etRVEC, {new_box},
 350           "Box size (in nm)" },
 351         { "-nmol",   FALSE, etINT, {&nmol_ins},
 352           "Number of extra molecules to insert" },
 353         { "-try",    FALSE, etINT, {&nmol_try},
 354           "Try inserting [TT]-nmol[tt] times [TT]-try[tt] times" },
 355         { "-seed",   FALSE, etINT, {&seed},
 356           "Random generator seed"},
 357         { "-radius",   FALSE, etREAL, {&defaultDistance},
 358           "Default van der Waals distance"},
 359         { "-scale", FALSE, etREAL, {&scaleFactor},
 360           "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." },
 361         { "-dr",    FALSE, etRVEC, {deltaR},
 362           "Allowed displacement in x/y/z from positions in [TT]-ip[tt] file" },
 363         { "-rot", FALSE,  etENUM, {enum_rot_string},
 364           "rotate inserted molecules randomly" },
 365         { "-allpair",    FALSE, etBOOL, {&bCheckAllPairDist},
 366           "Avoid momory leaks during neighbor searching with option -ci. May be slow for large systems." },
 367     };
 368
 369     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 370                            asize(desc), desc, asize(bugs), bugs, &oenv))
 371     {
 372         return 0;
 373     }
 374
 375     bProt     = opt2bSet("-f", NFILE, fnm);
 376     bBox      = opt2parg_bSet("-box", asize(pa), pa);
 377     enum_rot  = nenum(enum_rot_string);
 378
 379     /* check input */
 380     const char *insertionMoleculeFileName = opt2fn("-ci", NFILE, fnm);
 381     if (!gmx_fexist(insertionMoleculeFileName))
 382     {
 383         gmx_fatal(FARGS,
 384                   "A molecule conformation to insert is required in -ci. %s was not found!",
 385                   insertionMoleculeFileName);
 386     }
 387     if (nmol_ins <= 0 && !opt2bSet("-ip", NFILE, fnm))
 388     {
 389         gmx_fatal(FARGS, "Either -nmol must be larger than 0, "
 390                   "or positions must be given with -ip");
 391     }
 392     if (!bProt && !bBox)
 393     {
 394         gmx_fatal(FARGS, "When no solute (-f) is specified, "
 395                   "a box size (-box) must be specified");
 396     }
 397
 398     aps = gmx_atomprop_init();
 399
 400     snew(atoms, 1);
 401     init_t_atoms(atoms, 0, FALSE);
 402     if (bProt)
 403     {
 404         /* Generate a solute configuration */
 405         conf_prot = opt2fn("-f", NFILE, fnm);
 406         title     = readConformation(conf_prot, atoms, &x, NULL,
 407                                      &ePBC, box);
 408         exclusionDistances = makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor);
 409         if (atoms->nr == 0)
 410         {
 411             fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
 412             bProt = FALSE;
 413         }
 414     }
 415     if (bBox)
 416     {
 417         ePBC = epbcXYZ;
 418         clear_mat(box);
 419         box[XX][XX] = new_box[XX];
 420         box[YY][YY] = new_box[YY];
 421         box[ZZ][ZZ] = new_box[ZZ];
 422     }
 423     if (det(box) == 0)
 424     {
 425         gmx_fatal(FARGS, "Undefined solute box.\nCreate one with gmx editconf "
 426                   "or give explicit -box command line option");
 427     }
 428
 429     /* add nmol_ins molecules of atoms_ins
 430        in random orientation at random place */
 431     title_ins = insert_mols(insertionMoleculeFileName, nmol_ins, nmol_try, seed,
 432                             atoms, &x, &exclusionDistances, ePBC, box, aps,
 433                             defaultDistance, scaleFactor, 0,
 434                             oenv, opt2fn_null("-ip", NFILE, fnm), deltaR, enum_rot,
 435                             bCheckAllPairDist);
 436
 437     /* write new configuration to file confout */
 438     confout = ftp2fn(efSTO, NFILE, fnm);
 439     fprintf(stderr, "Writing generated configuration to %s\n", confout);
 440     if (bProt)
 441     {
 442         write_sto_conf(confout, title, atoms, x, NULL, ePBC, box);
 443         /* print box sizes and box type to stderr */
 444         fprintf(stderr, "%s\n", title);
 445     }
 446     else
 447     {
 448         write_sto_conf(confout, title_ins, atoms, x, NULL, ePBC, box);
 449     }
 450
 451     /* print size of generated configuration */
 452     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
 453             atoms->nr, atoms->nres);
 454
 455     gmx_atomprop_destroy(aps);
 456     sfree(exclusionDistances);
 457     sfree(x);
 458     done_atom(atoms);
 459     sfree(atoms);
 460
 461     return 0;
 462 }