| File: | gromacs/gmxana/gmx_nmtraj.c |
| Location: | line 180, column 9 |
| Description: | Value stored to 'dum' is never read |
| 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 | #ifdef HAVE_CONFIG_H1 |
| 38 | #include <config.h> |
| 39 | #endif |
| 40 | |
| 41 | #include <math.h> |
| 42 | #include <stdlib.h> |
| 43 | #include <string.h> |
| 44 | |
| 45 | #include "gromacs/commandline/pargs.h" |
| 46 | #include "typedefs.h" |
| 47 | #include "gromacs/utility/smalloc.h" |
| 48 | #include "macros.h" |
| 49 | #include "gromacs/utility/fatalerror.h" |
| 50 | #include "gromacs/math/vec.h" |
| 51 | #include "gromacs/utility/futil.h" |
| 52 | #include "index.h" |
| 53 | #include "gromacs/fileio/pdbio.h" |
| 54 | #include "gromacs/fileio/tpxio.h" |
| 55 | #include "gromacs/fileio/trxio.h" |
| 56 | #include "txtdump.h" |
| 57 | #include "physics.h" |
| 58 | #include "eigio.h" |
| 59 | #include "gmx_ana.h" |
| 60 | |
| 61 | |
| 62 | int gmx_nmtraj(int argc, char *argv[]) |
| 63 | { |
| 64 | const char *desc[] = |
| 65 | { |
| 66 | "[THISMODULE] generates an virtual trajectory from an eigenvector, ", |
| 67 | "corresponding to a harmonic Cartesian oscillation around the average ", |
| 68 | "structure. The eigenvectors should normally be mass-weighted, but you can ", |
| 69 | "use non-weighted eigenvectors to generate orthogonal motions. ", |
| 70 | "The output frames are written as a trajectory file covering an entire period, and ", |
| 71 | "the first frame is the average structure. If you write the trajectory in (or convert to) ", |
| 72 | "PDB format you can view it directly in PyMol and also render a photorealistic movie. ", |
| 73 | "Motion amplitudes are calculated from the eigenvalues and a preset temperature, ", |
| 74 | "assuming equipartition of the energy over all modes. To make the motion clearly visible ", |
| 75 | "in PyMol you might want to amplify it by setting an unrealistically high temperature. ", |
| 76 | "However, be aware that both the linear Cartesian displacements and mass weighting will ", |
| 77 | "lead to serious structure deformation for high amplitudes - this is is simply a limitation ", |
| 78 | "of the Cartesian normal mode model. By default the selected eigenvector is set to 7, since ", |
| 79 | " the first six normal modes are the translational and rotational degrees of freedom." |
| 80 | }; |
| 81 | |
| 82 | static real refamplitude = 0.25; |
| 83 | static int nframes = 30; |
| 84 | static real temp = 300.0; |
| 85 | static const char *eignrvec = "7"; |
| 86 | static const char *phasevec = "0.0"; |
| 87 | |
| 88 | t_pargs pa[] = |
| 89 | { |
| 90 | { "-eignr", FALSE0, etSTR, {&eignrvec}, "String of eigenvectors to use (first is 1)" }, |
| 91 | { "-phases", FALSE0, etSTR, {&phasevec}, "String of phases (default is 0.0)" }, |
| 92 | { "-temp", FALSE0, etREAL, {&temp}, "Temperature (K)" }, |
| 93 | { "-amplitude", FALSE0, etREAL, {&refamplitude}, "Amplitude for modes with eigenvalue<=0" }, |
| 94 | { "-nframes", FALSE0, etINT, {&nframes}, "Number of frames to generate" } |
| 95 | }; |
| 96 | |
| 97 | #define NPA((int)(sizeof(pa)/sizeof((pa)[0]))) asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))) |
| 98 | |
| 99 | t_trxstatus *out; |
| 100 | t_topology top; |
| 101 | int ePBC; |
| 102 | t_atoms *atoms; |
| 103 | rvec *xtop, *xref, *xav, *xout; |
| 104 | int nvec, *eignr = NULL((void*)0); |
| 105 | int *eigvalnr; |
| 106 | rvec **eigvec = NULL((void*)0); |
| 107 | matrix box; |
| 108 | int natoms; |
| 109 | int i, j, k, kmode, d, s, v; |
| 110 | gmx_bool bDMR, bDMA, bFit; |
| 111 | char * indexfile; |
| 112 | |
| 113 | char * grpname; |
| 114 | real * eigval; |
| 115 | int neigval; |
| 116 | int * dummy; |
| 117 | real * invsqrtm; |
| 118 | char title[STRLEN4096]; |
| 119 | real fraction; |
| 120 | int *out_eigidx; |
| 121 | real *out_eigval; |
| 122 | rvec * this_eigvec; |
| 123 | real omega, Ekin, sum, m, vel; |
| 124 | gmx_bool found; |
| 125 | int nmodes, nphases; |
| 126 | int *imodes; |
| 127 | real *amplitude; |
| 128 | real *phases; |
| 129 | real dum; |
| 130 | const char *p; |
| 131 | char *pe; |
| 132 | output_env_t oenv; |
| 133 | |
| 134 | t_filenm fnm[] = |
| 135 | { |
| 136 | { efTPS, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, |
| 137 | { efTRN, "-v", "eigenvec", ffREAD1<<1 }, |
| 138 | { efTRO, "-o", "nmtraj", ffWRITE1<<2 } |
| 139 | }; |
| 140 | |
| 141 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) |
| 142 | |
| 143 | if (!parse_common_args(&argc, argv, PCA_BE_NICE(1<<13), |
| 144 | NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv)) |
| 145 | { |
| 146 | return 0; |
| 147 | } |
| 148 | |
| 149 | read_eigenvectors(opt2fn("-v", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &natoms, &bFit, |
| 150 | &xref, &bDMR, &xav, &bDMA, &nvec, &eignr, &eigvec, &eigval); |
| 151 | |
| 152 | read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), title, &top, &ePBC, &xtop, NULL((void*)0), box, bDMA); |
| 153 | |
| 154 | /* Find vectors and phases */ |
| 155 | |
| 156 | /* first find number of args in string */ |
| 157 | nmodes = 0; |
| 158 | p = eignrvec; |
| 159 | while (*p != 0) |
| 160 | { |
| 161 | dum = strtod(p, &pe); |
| 162 | p = pe; |
| 163 | nmodes++; |
| 164 | } |
| 165 | |
| 166 | snew(imodes, nmodes)(imodes) = save_calloc("imodes", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmtraj.c" , 166, (nmodes), sizeof(*(imodes))); |
| 167 | p = eignrvec; |
| 168 | for (i = 0; i < nmodes; i++) |
| 169 | { |
| 170 | /* C indices start on 0 */ |
| 171 | imodes[i] = strtol(p, &pe, 10)-1; |
| 172 | p = pe; |
| 173 | } |
| 174 | |
| 175 | /* Now read phases */ |
| 176 | nphases = 0; |
| 177 | p = phasevec; |
| 178 | while (*p != 0) |
| 179 | { |
| 180 | dum = strtod(p, &pe); |
Value stored to 'dum' is never read | |
| 181 | p = pe; |
| 182 | nphases++; |
| 183 | } |
| 184 | if (nphases > nmodes) |
| 185 | { |
| 186 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmtraj.c" , 186, "More phases than eigenvector indices specified.\n"); |
| 187 | } |
| 188 | |
| 189 | snew(phases, nmodes)(phases) = save_calloc("phases", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmtraj.c" , 189, (nmodes), sizeof(*(phases))); |
| 190 | p = phasevec; |
| 191 | |
| 192 | for (i = 0; i < nphases; i++) |
| 193 | { |
| 194 | phases[i] = strtod(p, &pe); |
| 195 | p = pe; |
| 196 | } |
| 197 | |
| 198 | if (nmodes > nphases) |
| 199 | { |
| 200 | printf("Warning: Setting phase of last %d modes to zero...\n", nmodes-nphases); |
| 201 | } |
| 202 | |
| 203 | for (i = nphases; i < nmodes; i++) |
| 204 | { |
| 205 | phases[i] = 0; |
| 206 | } |
| 207 | |
| 208 | atoms = &top.atoms; |
| 209 | |
| 210 | if (atoms->nr != natoms) |
| 211 | { |
| 212 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmtraj.c" , 212, "Different number of atoms in topology and eigenvectors.\n"); |
| 213 | } |
| 214 | |
| 215 | snew(dummy, natoms)(dummy) = save_calloc("dummy", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmtraj.c" , 215, (natoms), sizeof(*(dummy))); |
| 216 | for (i = 0; i < natoms; i++) |
| 217 | { |
| 218 | dummy[i] = i; |
| 219 | } |
| 220 | |
| 221 | /* Find the eigenvalue/vector to match our select one */ |
| 222 | snew(out_eigidx, nmodes)(out_eigidx) = save_calloc("out_eigidx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmtraj.c" , 222, (nmodes), sizeof(*(out_eigidx))); |
| 223 | for (i = 0; i < nmodes; i++) |
| 224 | { |
| 225 | out_eigidx[i] = -1; |
| 226 | } |
| 227 | |
| 228 | for (i = 0; i < nvec; i++) |
| 229 | { |
| 230 | for (j = 0; j < nmodes; j++) |
| 231 | { |
| 232 | if (imodes[j] == eignr[i]) |
| 233 | { |
| 234 | out_eigidx[j] = i; |
| 235 | } |
| 236 | } |
| 237 | } |
| 238 | for (i = 0; i < nmodes; i++) |
| 239 | { |
| 240 | if (out_eigidx[i] == -1) |
| 241 | { |
| 242 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmtraj.c" , 242, "Could not find mode %d in eigenvector file.\n", imodes[i]); |
| 243 | } |
| 244 | } |
| 245 | |
| 246 | |
| 247 | snew(invsqrtm, natoms)(invsqrtm) = save_calloc("invsqrtm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmtraj.c" , 247, (natoms), sizeof(*(invsqrtm))); |
| 248 | |
| 249 | if (bDMA) |
| 250 | { |
| 251 | for (i = 0; (i < natoms); i++) |
| 252 | { |
| 253 | invsqrtm[i] = gmx_invsqrt(atoms->atom[i].m)gmx_software_invsqrt(atoms->atom[i].m); |
| 254 | } |
| 255 | } |
| 256 | else |
| 257 | { |
| 258 | for (i = 0; (i < natoms); i++) |
| 259 | { |
| 260 | invsqrtm[i] = 1.0; |
| 261 | } |
| 262 | } |
| 263 | |
| 264 | snew(xout, natoms)(xout) = save_calloc("xout", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmtraj.c" , 264, (natoms), sizeof(*(xout))); |
| 265 | snew(amplitude, nmodes)(amplitude) = save_calloc("amplitude", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmtraj.c" , 265, (nmodes), sizeof(*(amplitude))); |
| 266 | |
| 267 | printf("mode phases: %g %g\n", phases[0], phases[1]); |
| 268 | |
| 269 | for (i = 0; i < nmodes; i++) |
| 270 | { |
| 271 | kmode = out_eigidx[i]; |
| 272 | this_eigvec = eigvec[kmode]; |
| 273 | |
| 274 | if ( (kmode >= 6) && (eigval[kmode] > 0)) |
| 275 | { |
| 276 | /* Derive amplitude from temperature and eigenvalue if we can */ |
| 277 | |
| 278 | /* Convert eigenvalue to angular frequency, in units s^(-1) */ |
| 279 | omega = sqrt(eigval[kmode]*1.0E21/(AVOGADRO(6.0221367e23)*AMU(1.6605402e-27))); |
| 280 | /* Harmonic motion will be x=x0 + A*sin(omega*t)*eigenvec. |
| 281 | * The velocity is thus: |
| 282 | * |
| 283 | * v = A*omega*cos(omega*t)*eigenvec. |
| 284 | * |
| 285 | * And the average kinetic energy the integral of mass*v*v/2 over a |
| 286 | * period: |
| 287 | * |
| 288 | * (1/4)*mass*A*omega*eigenvec |
| 289 | * |
| 290 | * For t =2*pi*n, all energy will be kinetic, and v=A*omega*eigenvec. |
| 291 | * The kinetic energy will be sum(0.5*mass*v*v) if we temporarily set A to 1, |
| 292 | * and the average over a period half of this. |
| 293 | */ |
| 294 | |
| 295 | Ekin = 0; |
| 296 | for (k = 0; k < natoms; k++) |
| 297 | { |
| 298 | m = atoms->atom[k].m; |
| 299 | for (d = 0; d < DIM3; d++) |
| 300 | { |
| 301 | vel = omega*this_eigvec[k][d]; |
| 302 | Ekin += 0.5*0.5*m*vel*vel; |
| 303 | } |
| 304 | } |
| 305 | |
| 306 | /* Convert Ekin from amu*(nm/s)^2 to J, i.e., kg*(m/s)^2 |
| 307 | * This will also be proportional to A^2 |
| 308 | */ |
| 309 | Ekin *= AMU(1.6605402e-27)*1E-18; |
| 310 | |
| 311 | /* Set the amplitude so the energy is kT/2 */ |
| 312 | amplitude[i] = sqrt(0.5*BOLTZMANN(1.380658e-23)*temp/Ekin); |
| 313 | } |
| 314 | else |
| 315 | { |
| 316 | amplitude[i] = refamplitude; |
| 317 | } |
| 318 | } |
| 319 | |
| 320 | out = open_trx(ftp2fn(efTRO, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "w"); |
| 321 | |
| 322 | /* Write a sine oscillation around the average structure, |
| 323 | * modulated by the eigenvector with selected amplitude. |
| 324 | */ |
| 325 | |
| 326 | for (i = 0; i < nframes; i++) |
| 327 | { |
| 328 | fraction = (real)i/(real)nframes; |
| 329 | for (j = 0; j < natoms; j++) |
| 330 | { |
| 331 | copy_rvec(xav[j], xout[j]); |
| 332 | } |
| 333 | |
| 334 | for (k = 0; k < nmodes; k++) |
| 335 | { |
| 336 | kmode = out_eigidx[k]; |
| 337 | this_eigvec = eigvec[kmode]; |
| 338 | |
| 339 | for (j = 0; j < natoms; j++) |
| 340 | { |
| 341 | for (d = 0; d < DIM3; d++) |
| 342 | { |
| 343 | xout[j][d] += amplitude[k]*sin(2*M_PI3.14159265358979323846*(fraction+phases[k]/360.0))*this_eigvec[j][d]; |
| 344 | } |
| 345 | } |
| 346 | } |
| 347 | write_trx(out, natoms, dummy, atoms, i, (real)i/(real)nframes, box, xout, NULL((void*)0), NULL((void*)0)); |
| 348 | } |
| 349 | |
| 350 | fprintf(stderrstderr, "\n"); |
| 351 | close_trx(out); |
| 352 | |
| 353 | return 0; |
| 354 | } |