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 "smalloc.h"
  47 #include "string2.h"
  48 #include "gromacs/math/utilities.h"
  49 #include "gromacs/fileio/confio.h"
  50 #include "macros.h"
  51 #include "gromacs/random/random.h"
  52 #include "gromacs/fileio/futil.h"
  53 #include "atomprop.h"
  54 #include "names.h"
  55 #include "vec.h"
  56 #include "gmx_fatal.h"
  57 #include "gromacs/commandline/pargs.h"
  58 #include "gromacs/gmxlib/conformation-utilities.h"
  59 #include "addconf.h"
  60 #include "read-conformation.h"
  61 #include "pbc.h"
  62 #include "xvgr.h"
  63
  64 static gmx_bool in_box(t_pbc *pbc, rvec x)
  65 {
  66     rvec box_center, dx;
  67     int  shift;
  68
  69     /* pbc_dx_aiuc only works correctly with the rectangular box center */
  70     calc_box_center(ecenterRECT, pbc->box, box_center);
  71
  72     shift = pbc_dx_aiuc(pbc, x, box_center, dx);
  73
  74     return (shift == CENTRAL);
  75 }
  76
  77 /* This is a (maybe) slow workaround to avoid the neighbor searching in addconf.c, which
  78  * leaks memory (May 2012). The function could be deleted as soon as the memory leaks
  79  * there are fixed.
  80  * However, when inserting a small molecule in a system containing not too many atoms,
  81  * allPairsDistOk is probably even faster than the other code.
  82  */
  83 static gmx_bool
  84 allPairsDistOk(t_atoms *atoms, rvec *x, real *r,
  85                int ePBC, matrix box,
  86                t_atoms *atoms_insrt, rvec *x_n, real *r_insrt)
  87 {
  88     int   i, j;
  89     rvec  dx;
  90     real  n2, r2;
  91     t_pbc pbc;
  92
  93     set_pbc(&pbc, ePBC, box);
  94     for (i = 0; i < atoms->nr; i++)
  95     {
  96         for (j = 0; j < atoms_insrt->nr; j++)
  97         {
  98             pbc_dx(&pbc, x[i], x_n[j], dx);
  99             n2 = norm2(dx);
 100             r2 = sqr(r[i]+r_insrt[j]);
 101             if (n2 < r2)
 102             {
 103                 return FALSE;
 104             }
 105         }
 106     }
 107     return TRUE;
 108 }
 109
 110 /* enum for random rotations of inserted solutes */
 111 enum {
 112     en_rot, en_rotXYZ, en_rotZ, en_rotNone, en_NR
 113 };
 114
 115 static char *insert_mols(const char *mol_insrt, int nmol_insrt, int ntry, int seed,
 116                          t_atoms *atoms, rvec **x, real **r, int ePBC, matrix box,
 117                          gmx_atomprop_t aps,
 118                          real r_distance, real r_scale, real rshell,
 119                          const output_env_t oenv,
 120                          const char* posfn, const rvec deltaR, int enum_rot,
 121                          gmx_bool bCheckAllPairDist)
 122 {
 123     t_pbc            pbc;
 124     static  char    *title_insrt;
 125     t_atoms          atoms_insrt;
 126     rvec            *x_insrt, *x_n;
 127     real            *r_insrt;
 128     int              ePBC_insrt;
 129     matrix           box_insrt;
 130     int              i, mol, onr, ncol;
 131     real             alfa = 0., beta = 0., gamma = 0.;
 132     rvec             offset_x;
 133     int              trial;
 134     double         **rpos;
 135     gmx_rng_t        rng;
 136
 137     rng = gmx_rng_init(seed);
 138     set_pbc(&pbc, ePBC, box);
 139
 140     /* read number of atoms of insert molecules */
 141     get_stx_coordnum(mol_insrt, &atoms_insrt.nr);
 142     if (atoms_insrt.nr == 0)
 143     {
 144         gmx_fatal(FARGS, "No molecule in %s, please check your input\n", mol_insrt);
 145     }
 146     /* allocate memory for atom coordinates of insert molecules */
 147     snew(x_insrt, atoms_insrt.nr);
 148     snew(r_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 van der waals arrays for inserted molecules */
 165     mk_vdw(&atoms_insrt, r_insrt, aps, r_distance, r_scale);
 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(*r, (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, *r, ePBC, box, &atoms_insrt, x_n, r_insrt))
 244         {
 245             continue;
 246         }
 247
 248         add_conf(atoms, x, NULL, r, FALSE, ePBC, box, TRUE,
 249                  &atoms_insrt, x_n, NULL, r_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(*r,              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     return title_insrt;
 270 }
 271
 272 int gmx_insert_molecules(int argc, char *argv[])
 273 {
 274     const char *desc[] = {
 275         "[THISMODULE] inserts [TT]-nmol[tt] copies of the system specified in",
 276         "the [TT]-ci[tt] input file. The insertions take place either into",
 277         "vacant space in the solute conformation given with [TT]-f[tt], or",
 278         "into an empty box given by [TT]-box[tt]. Specifying both [TT]-f[tt]",
 279         "and [TT]-box[tt] behaves like [TT]-f[tt], but places a new box",
 280         "around the solute before insertions. Any velocities present are",
 281         "discarded.[PAR]",
 282
 283         "By default, the insertion positions are random (with initial seed",
 284         "specified by [TT]-seed[tt]). The program iterates until [TT]-nmol[tt]",
 285         "molecules have been inserted in the box. Molecules are not inserted",
 286         "where the distance between any existing atom and any atom of ",
 287         "the inserted molecule is less than the sum of the van der Waals radii of ",
 288         "both atoms. A database ([TT]vdwradii.dat[tt]) of van der Waals radii is ",
 289         "read by the program, and atoms not in the database are ",
 290         "assigned a default distance [TT]-vdwd[tt].[PAR]",
 291
 292         "A total of [TT]-nmol[tt] * [TT]-try[tt] insertion attempts are made",
 293         "before giving up. Increase [TT]-try[tt] if you have several small",
 294         "holes to fill. Option [TT]-rot[tt] specifies whether the insertion",
 295         "molecules are randomly oriented before insertion attempts.[PAR]",
 296
 297         "Alternatively, the molecules can be inserted only at positions defined in",
 298         "positions.dat ([TT]-ip[tt]). That file should have 3 columns (x,y,z),",
 299         "that give the displacements compared to the input molecule position",
 300         "([TT]-ci[tt]). Hence, if that file should contain the absolute",
 301         "positions, the molecule must be centered on (0,0,0) before using",
 302         "[THISMODULE] (e.g. from [gmx-editconf] [TT]-center[tt]).",
 303         "Comments in that file starting with # are ignored. Option [TT]-dr[tt]",
 304         "defines the maximally allowed displacements during insertial trials.",
 305         "[TT]-try[tt] and [TT]-rot[tt] work as in the default mode (see above).",
 306         "[PAR]",
 307     };
 308
 309     const char *bugs[] = {
 310         "Molecules must be whole in the initial configurations.",
 311         "Many repeated neighbor searchings with -ci blows up the allocated memory. "
 312         "Option -allpair avoids this using all-to-all distance checks (slow for large systems)"
 313     };
 314
 315     /* parameter data */
 316     gmx_bool       bProt, bBox;
 317     const char    *conf_prot, *confout;
 318     real          *r;
 319     char          *title_ins = NULL;
 320     gmx_atomprop_t aps;
 321
 322     /* protein configuration data */
 323     char          *title = NULL;
 324     t_atoms        atoms;
 325     rvec          *x    = NULL;
 326     int            ePBC = -1;
 327     matrix         box;
 328
 329     t_filenm       fnm[] = {
 330         { efSTX, "-f", "protein", ffOPTRD },
 331         { efSTX, "-ci", "insert",  ffREAD},
 332         { efDAT, "-ip", "positions",  ffOPTRD},
 333         { efSTO, NULL,  NULL,      ffWRITE},
 334     };
 335 #define NFILE asize(fnm)
 336
 337     static int      nmol_ins          = 0, nmol_try = 10, seed = 1997, enum_rot;
 338     static real     r_distance        = 0.105, r_scale = 0.57;
 339     static rvec     new_box           = {0.0, 0.0, 0.0}, deltaR = {0.0, 0.0, 0.0};
 340     static gmx_bool bCheckAllPairDist = FALSE;
 341     output_env_t    oenv;
 342     const char     *enum_rot_string[] = {NULL, "xyz", "z", "none", NULL};
 343     t_pargs         pa[]              = {
 344         { "-box",    FALSE, etRVEC, {new_box},
 345           "Box size (in nm)" },
 346         { "-nmol",   FALSE, etINT, {&nmol_ins},
 347           "Number of extra molecules to insert" },
 348         { "-try",    FALSE, etINT, {&nmol_try},
 349           "Try inserting [TT]-nmol[tt] times [TT]-try[tt] times" },
 350         { "-seed",   FALSE, etINT, {&seed},
 351           "Random generator seed"},
 352         { "-vdwd",   FALSE, etREAL, {&r_distance},
 353           "Default van der Waals distance"},
 354         { "-vdwscale", FALSE, etREAL, {&r_scale},
 355           "HIDDENScale 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." },
 356         { "-dr",    FALSE, etRVEC, {deltaR},
 357           "Allowed displacement in x/y/z from positions in [TT]-ip[tt] file" },
 358         { "-rot", FALSE,  etENUM, {enum_rot_string},
 359           "rotate inserted molecules randomly" },
 360         { "-allpair",    FALSE, etBOOL, {&bCheckAllPairDist},
 361           "Avoid momory leaks during neighbor searching with option -ci. May be slow for large systems." },
 362     };
 363
 364     if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
 365                            asize(desc), desc, asize(bugs), bugs, &oenv))
 366     {
 367         return 0;
 368     }
 369
 370     bProt     = opt2bSet("-f", NFILE, fnm);
 371     bBox      = opt2parg_bSet("-box", asize(pa), pa);
 372     enum_rot  = nenum(enum_rot_string);
 373
 374     /* check input */
 375     const char *insertionMoleculeFileName = opt2fn("-ci", NFILE, fnm);
 376     if (!gmx_fexist(insertionMoleculeFileName))
 377     {
 378         gmx_fatal(FARGS,
 379                   "A molecule conformation to insert is required in -ci. %s was not found!",
 380                   insertionMoleculeFileName);
 381     }
 382     if (nmol_ins <= 0 && !opt2bSet("-ip", NFILE, fnm))
 383     {
 384         gmx_fatal(FARGS, "Either -nmol must be larger than 0, "
 385                   "or positions must be given with -ip");
 386     }
 387     if (!bProt && !bBox)
 388     {
 389         gmx_fatal(FARGS, "When no solute (-f) is specified, "
 390                   "a box size (-box) must be specified");
 391     }
 392
 393     aps = gmx_atomprop_init();
 394
 395     if (bProt)
 396     {
 397         /*generate a solute configuration */
 398         conf_prot = opt2fn("-f", NFILE, fnm);
 399         title     = read_conformation(conf_prot, &atoms, &x, NULL, &r, &ePBC, box,
 400                                       aps, r_distance, r_scale);
 401         if (atoms.nr == 0)
 402         {
 403             fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
 404             bProt = FALSE;
 405         }
 406     }
 407     else
 408     {
 409         atoms.nr       = 0;
 410         atoms.nres     = 0;
 411         atoms.resinfo  = NULL;
 412         atoms.atomname = NULL;
 413         atoms.atom     = NULL;
 414         atoms.pdbinfo  = NULL;
 415         x              = NULL;
 416         r              = NULL;
 417     }
 418     if (bBox)
 419     {
 420         ePBC = epbcXYZ;
 421         clear_mat(box);
 422         box[XX][XX] = new_box[XX];
 423         box[YY][YY] = new_box[YY];
 424         box[ZZ][ZZ] = new_box[ZZ];
 425     }
 426     if (det(box) == 0)
 427     {
 428         gmx_fatal(FARGS, "Undefined solute box.\nCreate one with gmx editconf "
 429                   "or give explicit -box command line option");
 430     }
 431
 432     /* add nmol_ins molecules of atoms_ins
 433        in random orientation at random place */
 434     title_ins = insert_mols(insertionMoleculeFileName, nmol_ins, nmol_try, seed,
 435                             &atoms, &x, &r, ePBC, box, aps,
 436                             r_distance, r_scale, 0,
 437                             oenv, opt2fn_null("-ip", NFILE, fnm), deltaR, enum_rot,
 438                             bCheckAllPairDist);
 439
 440     /* write new configuration to file confout */
 441     confout = ftp2fn(efSTO, NFILE, fnm);
 442     fprintf(stderr, "Writing generated configuration to %s\n", confout);
 443     if (bProt)
 444     {
 445         write_sto_conf(confout, title, &atoms, x, NULL, ePBC, box);
 446         /* print box sizes and box type to stderr */
 447         fprintf(stderr, "%s\n", title);
 448     }
 449     else
 450     {
 451         write_sto_conf(confout, title_ins, &atoms, x, NULL, ePBC, box);
 452     }
 453
 454     /* print size of generated configuration */
 455     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
 456             atoms.nr, atoms.nres);
 457
 458     gmx_atomprop_destroy(aps);
 459
 460     return 0;
 461 }