src/gromacs/gmxpreprocess/genconf.c

   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 "genconf.h"
  40
  41 #include "gromacs/commandline/pargs.h"
  42 #include "gromacs/fileio/confio.h"
  43 #include "gromacs/fileio/trxio.h"
  44 #include "gromacs/gmxpreprocess/sortwater.h"
  45 #include "gromacs/legacyheaders/macros.h"
  46 #include "gromacs/legacyheaders/names.h"
  47 #include "gromacs/legacyheaders/readinp.h"
  48 #include "gromacs/legacyheaders/txtdump.h"
  49 #include "gromacs/math/3dtransforms.h"
  50 #include "gromacs/math/utilities.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/random/random.h"
  53 #include "gromacs/utility/fatalerror.h"
  54 #include "gromacs/utility/smalloc.h"
  55
  56 static void rand_rot(int natoms, rvec x[], rvec v[], vec4 xrot[], vec4 vrot[],
  57                      gmx_rng_t rng, rvec max_rot)
  58 {
  59     mat4 mt1, mt2, mr[DIM], mtemp1, mtemp2, mtemp3, mxtot, mvtot;
  60     rvec xcm;
  61     real phi;
  62     int  i, m;
  63
  64     clear_rvec(xcm);
  65     for (i = 0; (i < natoms); i++)
  66     {
  67         for (m = 0; (m < DIM); m++)
  68         {
  69             xcm[m] += x[i][m]/natoms; /* get center of mass of one molecule  */
  70         }
  71     }
  72     fprintf(stderr, "center of geometry: %f, %f, %f\n", xcm[0], xcm[1], xcm[2]);
  73
  74     /* move c.o.ma to origin */
  75     gmx_mat4_init_translation(-xcm[XX], -xcm[YY], -xcm[ZZ], mt1);
  76     for (m = 0; (m < DIM); m++)
  77     {
  78         phi = M_PI*max_rot[m]*(2*gmx_rng_uniform_real(rng) - 1)/180;
  79         gmx_mat4_init_rotation(m, phi, mr[m]);
  80     }
  81     gmx_mat4_init_translation(xcm[XX], xcm[YY], xcm[ZZ], mt2);
  82
  83     /* For gmx_mat4_mmul() we need to multiply in the opposite order
  84      * compared to normal mathematical notation.
  85      */
  86     gmx_mat4_mmul(mtemp1, mt1, mr[XX]);
  87     gmx_mat4_mmul(mtemp2, mr[YY], mr[ZZ]);
  88     gmx_mat4_mmul(mtemp3, mtemp1, mtemp2);
  89     gmx_mat4_mmul(mxtot, mtemp3, mt2);
  90     gmx_mat4_mmul(mvtot, mr[XX], mtemp2);
  91
  92     for (i = 0; (i < natoms); i++)
  93     {
  94         gmx_mat4_transform_point(mxtot, x[i], xrot[i]);
  95         gmx_mat4_transform_point(mvtot, v[i], vrot[i]);
  96     }
  97 }
  98
  99 static void move_x(int natoms, rvec x[], matrix box)
 100 {
 101     int  i, m;
 102     rvec xcm;
 103
 104     clear_rvec(xcm);
 105     for (i = 0; (i < natoms); i++)
 106     {
 107         for (m = 0; (m < DIM); m++)
 108         {
 109             xcm[m] += x[i][m];
 110         }
 111     }
 112     for (m = 0; (m < DIM); m++)
 113     {
 114         xcm[m] = 0.5*box[m][m]-xcm[m]/natoms;
 115     }
 116     for (i = 0; (i < natoms); i++)
 117     {
 118         for (m = 0; (m < DIM); m++)
 119         {
 120             x[i][m] += xcm[m];
 121         }
 122     }
 123 }
 124
 125 int gmx_genconf(int argc, char *argv[])
 126 {
 127     const char     *desc[] = {
 128         "[THISMODULE] multiplies a given coordinate file by simply stacking them",
 129         "on top of each other, like a small child playing with wooden blocks.",
 130         "The program makes a grid of [IT]user-defined[it]",
 131         "proportions ([TT]-nbox[tt]), ",
 132         "and interspaces the grid point with an extra space [TT]-dist[tt].[PAR]",
 133         "When option [TT]-rot[tt] is used the program does not check for overlap",
 134         "between molecules on grid points. It is recommended to make the box in",
 135         "the input file at least as big as the coordinates + ",
 136         "van der Waals radius.[PAR]",
 137         "If the optional trajectory file is given, conformations are not",
 138         "generated, but read from this file and translated appropriately to",
 139         "build the grid."
 140
 141     };
 142     const char     *bugs[] = {
 143         "The program should allow for random displacement of lattice points."
 144     };
 145
 146     int             vol;
 147     t_atoms        *atoms;      /* list with all atoms */
 148     char            title[STRLEN];
 149     rvec           *x, *xx, *v; /* coordinates? */
 150     real            t;
 151     vec4           *xrot, *vrot;
 152     int             ePBC;
 153     matrix          box, boxx; /* box length matrix */
 154     rvec            shift;
 155     int             natoms;    /* number of atoms in one molecule  */
 156     int             nres;      /* number of molecules? */
 157     int             i, j, k, l, m, ndx, nrdx, nx, ny, nz;
 158     t_trxstatus    *status;
 159     gmx_bool        bTRX;
 160     output_env_t    oenv;
 161     gmx_rng_t       rng;
 162
 163     t_filenm        fnm[] = {
 164         { efSTX, "-f", "conf", ffREAD  },
 165         { efSTO, "-o", "out",  ffWRITE },
 166         { efTRX, "-trj", NULL,  ffOPTRD }
 167     };
 168 #define NFILE asize(fnm)
 169     static rvec     nrbox    = {1, 1, 1};
 170     static int      seed     = 0;    /* seed for random number generator */
 171     static int      nmolat   = 3;
 172     static int      nblock   = 1;
 173     static gmx_bool bShuffle = FALSE;
 174     static gmx_bool bSort    = FALSE;
 175     static gmx_bool bRandom  = FALSE;           /* False: no random rotations */
 176     static gmx_bool bRenum   = TRUE;            /* renumber residues */
 177     static rvec     dist     = {0, 0, 0};       /* space added between molecules ? */
 178     static rvec     max_rot  = {180, 180, 180}; /* maximum rotation */
 179     t_pargs         pa[]     = {
 180         { "-nbox",   FALSE, etRVEC, {nrbox},   "Number of boxes" },
 181         { "-dist",   FALSE, etRVEC, {dist},    "Distance between boxes" },
 182         { "-seed",   FALSE, etINT,  {&seed},
 183           "Random generator seed, if 0 generated from the time" },
 184         { "-rot",    FALSE, etBOOL, {&bRandom}, "Randomly rotate conformations" },
 185         { "-shuffle", FALSE, etBOOL, {&bShuffle}, "Random shuffling of molecules" },
 186         { "-sort",   FALSE, etBOOL, {&bSort},   "Sort molecules on X coord" },
 187         { "-block",  FALSE, etINT,  {&nblock},
 188           "Divide the box in blocks on this number of cpus" },
 189         { "-nmolat", FALSE, etINT,  {&nmolat},
 190           "Number of atoms per molecule, assumed to start from 0. "
 191           "If you set this wrong, it will screw up your system!" },
 192         { "-maxrot", FALSE, etRVEC, {max_rot}, "Maximum random rotation" },
 193         { "-renumber", FALSE, etBOOL, {&bRenum},  "Renumber residues" }
 194     };
 195
 196     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 197                            asize(desc), desc, asize(bugs), bugs, &oenv))
 198     {
 199         return 0;
 200     }
 201
 202     if (seed == 0)
 203     {
 204         rng = gmx_rng_init(gmx_rng_make_seed());
 205     }
 206     else
 207     {
 208         rng = gmx_rng_init(seed);
 209     }
 210
 211     bTRX = ftp2bSet(efTRX, NFILE, fnm);
 212     nx   = (int)(nrbox[XX]+0.5);
 213     ny   = (int)(nrbox[YY]+0.5);
 214     nz   = (int)(nrbox[ZZ]+0.5);
 215
 216     if ((nx <= 0) || (ny <= 0) || (nz <= 0))
 217     {
 218         gmx_fatal(FARGS, "Number of boxes (-nbox) should be larger than zero");
 219     }
 220     if ((nmolat <= 0) && bShuffle)
 221     {
 222         gmx_fatal(FARGS, "Can not shuffle if the molecules only have %d atoms",
 223                   nmolat);
 224     }
 225
 226     vol = nx*ny*nz; /* calculate volume in grid points (= nr. molecules) */
 227
 228     get_stx_coordnum(opt2fn("-f", NFILE, fnm), &natoms);
 229     snew(atoms, 1);
 230     /* make space for all the atoms */
 231     init_t_atoms(atoms, natoms*vol, FALSE);
 232     snew(x, natoms*vol);           /* get space for coordinates of all atoms */
 233     snew(xrot, natoms);            /* get space for rotation matrix? */
 234     snew(v, natoms*vol);           /* velocities. not really needed? */
 235     snew(vrot, natoms);
 236     /* set atoms->nr to the number in one box *
 237      * to avoid complaints in read_stx_conf   *
 238      */
 239     atoms->nr = natoms;
 240     read_stx_conf(opt2fn("-f", NFILE, fnm), title, atoms, x, v, &ePBC, box);
 241
 242     nres = atoms->nres;            /* nr of residues in one element? */
 243
 244     if (bTRX)
 245     {
 246         if (!read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &xx, boxx))
 247         {
 248             gmx_fatal(FARGS, "No atoms in trajectory %s", ftp2fn(efTRX, NFILE, fnm));
 249         }
 250     }
 251     else
 252     {
 253         snew(xx, natoms);
 254         for (i = 0; i < natoms; i++)
 255         {
 256             copy_rvec(x[i], xx[i]);
 257         }
 258     }
 259
 260
 261     for (k = 0; (k < nz); k++)     /* loop over all gridpositions    */
 262     {
 263         shift[ZZ] = k*(dist[ZZ]+box[ZZ][ZZ]);
 264
 265         for (j = 0; (j < ny); j++)
 266         {
 267             shift[YY] = j*(dist[YY]+box[YY][YY])+k*box[ZZ][YY];
 268
 269             for (i = 0; (i < nx); i++)
 270             {
 271                 shift[XX] = i*(dist[XX]+box[XX][XX])+j*box[YY][XX]+k*box[ZZ][XX];
 272
 273                 ndx  = (i*ny*nz+j*nz+k)*natoms;
 274                 nrdx = (i*ny*nz+j*nz+k)*nres;
 275
 276                 /* Random rotation on input coords */
 277                 if (bRandom)
 278                 {
 279                     rand_rot(natoms, xx, v, xrot, vrot, rng, max_rot);
 280                 }
 281
 282                 for (l = 0; (l < natoms); l++)
 283                 {
 284                     for (m = 0; (m < DIM); m++)
 285                     {
 286                         if (bRandom)
 287                         {
 288                             x[ndx+l][m] = xrot[l][m];
 289                             v[ndx+l][m] = vrot[l][m];
 290                         }
 291                         else
 292                         {
 293                             x[ndx+l][m] = xx[l][m];
 294                             v[ndx+l][m] = v[l][m];
 295                         }
 296                     }
 297                     if (ePBC == epbcSCREW && i % 2 == 1)
 298                     {
 299                         /* Rotate around x axis */
 300                         for (m = YY; m <= ZZ; m++)
 301                         {
 302                             x[ndx+l][m] = box[YY][m] + box[ZZ][m] - x[ndx+l][m];
 303                             v[ndx+l][m] = -v[ndx+l][m];
 304                         }
 305                     }
 306                     for (m = 0; (m < DIM); m++)
 307                     {
 308                         x[ndx+l][m] += shift[m];
 309                     }
 310                     atoms->atom[ndx+l].resind = nrdx + atoms->atom[l].resind;
 311                     atoms->atomname[ndx+l]    = atoms->atomname[l];
 312                 }
 313
 314                 for (l = 0; (l < nres); l++)
 315                 {
 316                     atoms->resinfo[nrdx+l] = atoms->resinfo[l];
 317                     if (bRenum)
 318                     {
 319                         atoms->resinfo[nrdx+l].nr += nrdx;
 320                     }
 321                 }
 322                 if (bTRX)
 323                 {
 324                     if (!read_next_x(oenv, status, &t, xx, boxx) &&
 325                         ((i+1)*(j+1)*(k+1) < vol))
 326                     {
 327                         gmx_fatal(FARGS, "Not enough frames in trajectory");
 328                     }
 329                 }
 330             }
 331         }
 332     }
 333     if (bTRX)
 334     {
 335         close_trj(status);
 336     }
 337
 338     /* make box bigger */
 339     for (m = 0; (m < DIM); m++)
 340     {
 341         box[m][m] += dist[m];
 342     }
 343     svmul(nx, box[XX], box[XX]);
 344     svmul(ny, box[YY], box[YY]);
 345     svmul(nz, box[ZZ], box[ZZ]);
 346     if (ePBC == epbcSCREW && nx % 2 == 0)
 347     {
 348         /* With an even number of boxes in x we can forgot about the screw */
 349         ePBC = epbcXYZ;
 350     }
 351
 352     /* move_x(natoms*vol,x,box); */          /* put atoms in box? */
 353
 354     atoms->nr   *= vol;
 355     atoms->nres *= vol;
 356
 357     /*depending on how you look at it, this is either a nasty hack or the way it should work*/
 358     if (bRenum)
 359     {
 360         for (i = 0; i < atoms->nres; i++)
 361         {
 362             atoms->resinfo[i].nr = i+1;
 363         }
 364     }
 365
 366
 367     if (bShuffle)
 368     {
 369         randwater(0, atoms->nr/nmolat, nmolat, x, v, rng);
 370     }
 371     else if (bSort)
 372     {
 373         sortwater(0, atoms->nr/nmolat, nmolat, x, v);
 374     }
 375     else if (opt2parg_bSet("-block", asize(pa), pa))
 376     {
 377         mkcompact(0, atoms->nr/nmolat, nmolat, x, v, nblock, box);
 378     }
 379     gmx_rng_destroy(rng);
 380
 381     write_sto_conf(opt2fn("-o", NFILE, fnm), title, atoms, x, v, ePBC, box);
 382
 383     return 0;
 384 }