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