src/gromacs/gmxana/gmx_nmens.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,2015,2016,2017,2019, 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 <cmath>
  40 #include <cstring>
  41
  42 #include "gromacs/commandline/pargs.h"
  43 #include "gromacs/fileio/confio.h"
  44 #include "gromacs/fileio/pdbio.h"
  45 #include "gromacs/fileio/trxio.h"
  46 #include "gromacs/gmxana/eigio.h"
  47 #include "gromacs/gmxana/gmx_ana.h"
  48 #include "gromacs/math/functions.h"
  49 #include "gromacs/math/units.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/random/threefry.h"
  52 #include "gromacs/random/uniformintdistribution.h"
  53 #include "gromacs/topology/index.h"
  54 #include "gromacs/topology/topology.h"
  55 #include "gromacs/utility/arraysize.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
  62 int gmx_nmens(int argc, char* argv[])
  63 {
  64     const char* desc[] = {
  65         "[THISMODULE] generates an ensemble around an average structure",
  66         "in a subspace that is defined by a set of normal modes (eigenvectors).",
  67         "The eigenvectors are assumed to be mass-weighted.",
  68         "The position along each eigenvector is randomly taken from a Gaussian",
  69         "distribution with variance kT/eigenvalue.[PAR]",
  70         "By default the starting eigenvector is set to 7, since the first six",
  71         "normal modes are the translational and rotational degrees of freedom."
  72     };
  73     static int  nstruct = 100, first = 7, last = -1, seed = 0;
  74     static real temp = 300.0;
  75     t_pargs     pa[] = {
  76         { "-temp", FALSE, etREAL, { &temp }, "Temperature in Kelvin" },
  77         { "-seed", FALSE, etINT, { &seed }, "Random seed (0 means generate)" },
  78         { "-num", FALSE, etINT, { &nstruct }, "Number of structures to generate" },
  79         { "-first", FALSE, etINT, { &first }, "First eigenvector to use (-1 is select)" },
  80         { "-last", FALSE, etINT, { &last }, "Last eigenvector to use (-1 is till the last)" }
  81     };
  82 #define NPA asize(pa)
  83
  84     t_trxstatus*        out;
  85     t_topology          top;
  86     int                 ePBC;
  87     t_atoms*            atoms;
  88     rvec *              xtop, *xref, *xav, *xout1, *xout2;
  89     gmx_bool            bDMR, bDMA, bFit;
  90     int                 nvec, *eignr = nullptr;
  91     rvec**              eigvec = nullptr;
  92     matrix              box;
  93     real *              eigval, *invsqrtm, t, disp;
  94     int                 natoms;
  95     char*               grpname;
  96     const char*         indexfile;
  97     int                 i, j, d, s, v;
  98     int                 nout, *iout, noutvec, *outvec;
  99     int*                index;
 100     real                rfac, rhalf, jr;
 101     gmx_output_env_t*   oenv;
 102     int                 jran;
 103     const unsigned long im = 0xffff;
 104     const unsigned long ia = 1093;
 105     const unsigned long ic = 18257;
 106
 107
 108     t_filenm fnm[] = { { efTRN, "-v", "eigenvec", ffREAD },
 109                        { efXVG, "-e", "eigenval", ffREAD },
 110                        { efTPS, nullptr, nullptr, ffREAD },
 111                        { efNDX, nullptr, nullptr, ffOPTRD },
 112                        { efTRO, "-o", "ensemble", ffWRITE } };
 113 #define NFILE asize(fnm)
 114
 115     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
 116     {
 117         return 0;
 118     }
 119
 120     indexfile = ftp2fn_null(efNDX, NFILE, fnm);
 121
 122     read_eigenvectors(opt2fn("-v", NFILE, fnm), &natoms, &bFit, &xref, &bDMR, &xav, &bDMA, &nvec,
 123                       &eignr, &eigvec, &eigval);
 124
 125     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xtop, nullptr, box, bDMA);
 126     atoms = &top.atoms;
 127
 128     printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
 129     get_index(atoms, indexfile, 1, &i, &index, &grpname);
 130     if (i != natoms)
 131     {
 132         gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
 133     }
 134     printf("\n");
 135
 136     snew(invsqrtm, natoms);
 137     if (bDMA)
 138     {
 139         for (i = 0; (i < natoms); i++)
 140         {
 141             invsqrtm[i] = gmx::invsqrt(atoms->atom[index[i]].m);
 142         }
 143     }
 144     else
 145     {
 146         for (i = 0; (i < natoms); i++)
 147         {
 148             invsqrtm[i] = 1.0;
 149         }
 150     }
 151
 152     if (last == -1)
 153     {
 154         last = natoms * DIM;
 155     }
 156     if (first > -1)
 157     {
 158         /* make an index from first to last */
 159         nout = last - first + 1;
 160         snew(iout, nout);
 161         for (i = 0; i < nout; i++)
 162         {
 163             iout[i] = first - 1 + i;
 164         }
 165     }
 166     else
 167     {
 168         printf("Select eigenvectors for output, end your selection with 0\n");
 169         nout = -1;
 170         iout = nullptr;
 171         do
 172         {
 173             nout++;
 174             srenew(iout, nout + 1);
 175             if (1 != scanf("%d", &iout[nout]))
 176             {
 177                 gmx_fatal(FARGS, "Error reading user input");
 178             }
 179             iout[nout]--;
 180         } while (iout[nout] >= 0);
 181         printf("\n");
 182     }
 183
 184     /* make an index of the eigenvectors which are present */
 185     snew(outvec, nout);
 186     noutvec = 0;
 187     for (i = 0; i < nout; i++)
 188     {
 189         j = 0;
 190         while ((j < nvec) && (eignr[j] != iout[i]))
 191         {
 192             j++;
 193         }
 194         if ((j < nvec) && (eignr[j] == iout[i]))
 195         {
 196             outvec[noutvec] = j;
 197             iout[noutvec]   = iout[i];
 198             noutvec++;
 199         }
 200     }
 201
 202     fprintf(stderr, "%d eigenvectors selected for output\n", noutvec);
 203
 204
 205     if (seed == 0)
 206     {
 207         // Make do with 32 bits for now to avoid changing user input to hex
 208         seed = static_cast<int>(gmx::makeRandomSeed());
 209     }
 210
 211     gmx::DefaultRandomEngine rng(seed);
 212
 213     fprintf(stderr, "Using random seed %d and a temperature of %g K.\n", seed, temp);
 214
 215     gmx::UniformIntDistribution<int> dist(0, im - 1);
 216     jran = dist(rng);
 217
 218     snew(xout1, natoms);
 219     snew(xout2, atoms->nr);
 220     out = open_trx(ftp2fn(efTRO, NFILE, fnm), "w");
 221
 222     for (s = 0; s < nstruct; s++)
 223     {
 224         for (i = 0; i < natoms; i++)
 225         {
 226             copy_rvec(xav[i], xout1[i]);
 227         }
 228         for (j = 0; j < noutvec; j++)
 229         {
 230             v = outvec[j];
 231             /* (r-0.5) n times:  var_n = n * var_1 = n/12
 232                n=4:  var_n = 1/3, so multiply with 3 */
 233
 234             rfac  = std::sqrt(3.0 * BOLTZ * temp / eigval[iout[j]]);
 235             rhalf = 2.0 * rfac;
 236             rfac  = rfac / im;
 237
 238             jran = (jran * ia + ic) & im;
 239             jr   = jran;
 240             jran = (jran * ia + ic) & im;
 241             jr += jran;
 242             jran = (jran * ia + ic) & im;
 243             jr += jran;
 244             jran = (jran * ia + ic) & im;
 245             jr += jran;
 246             disp = rfac * jr - rhalf;
 247
 248             for (i = 0; i < natoms; i++)
 249             {
 250                 for (d = 0; d < DIM; d++)
 251                 {
 252                     xout1[i][d] += disp * eigvec[v][i][d] * invsqrtm[i];
 253                 }
 254             }
 255         }
 256         for (i = 0; i < natoms; i++)
 257         {
 258             copy_rvec(xout1[i], xout2[index[i]]);
 259         }
 260         t = s + 1;
 261         write_trx(out, natoms, index, atoms, 0, t, box, xout2, nullptr, nullptr);
 262         fprintf(stderr, "\rGenerated %d structures", s + 1);
 263         fflush(stderr);
 264     }
 265     fprintf(stderr, "\n");
 266     close_trx(out);
 267
 268     return 0;
 269 }