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