| File: | gromacs/gmxana/gmx_anaeig.c |
| Location: | line 750, column 16 |
| Description: | Potential leak of memory pointed to by 'resnm' |
| 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 "pbc.h" | |||
| 52 | #include "gromacs/utility/futil.h" | |||
| 53 | #include "index.h" | |||
| 54 | #include "gromacs/fileio/pdbio.h" | |||
| 55 | #include "gromacs/fileio/confio.h" | |||
| 56 | #include "gromacs/fileio/tpxio.h" | |||
| 57 | #include "gromacs/fileio/trxio.h" | |||
| 58 | #include "gromacs/fileio/matio.h" | |||
| 59 | #include "mshift.h" | |||
| 60 | #include "gromacs/fileio/xvgr.h" | |||
| 61 | #include "viewit.h" | |||
| 62 | #include "rmpbc.h" | |||
| 63 | #include "txtdump.h" | |||
| 64 | #include "eigio.h" | |||
| 65 | #include "physics.h" | |||
| 66 | #include "gmx_ana.h" | |||
| 67 | ||||
| 68 | #include "gromacs/math/do_fit.h" | |||
| 69 | ||||
| 70 | static void calc_entropy_qh(FILE *fp, int n, real eigval[], real temp, int nskip) | |||
| 71 | { | |||
| 72 | int i; | |||
| 73 | double hwkT, w, dS, S = 0; | |||
| 74 | double hbar, lambda; | |||
| 75 | ||||
| 76 | hbar = PLANCK1(6.6262e-34)/(2*M_PI3.14159265358979323846); | |||
| 77 | for (i = 0; (i < n-nskip); i++) | |||
| 78 | { | |||
| 79 | if (eigval[i] > 0) | |||
| 80 | { | |||
| 81 | lambda = eigval[i]*AMU(1.6605402e-27); | |||
| 82 | w = sqrt(BOLTZMANN(1.380658e-23)*temp/lambda)/NANO(1e-9); | |||
| 83 | hwkT = (hbar*w)/(BOLTZMANN(1.380658e-23)*temp); | |||
| 84 | dS = (hwkT/(exp(hwkT) - 1) - log(1-exp(-hwkT))); | |||
| 85 | S += dS; | |||
| 86 | if (debug) | |||
| 87 | { | |||
| 88 | fprintf(debug, "i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n", | |||
| 89 | i, w, lambda, hwkT, dS); | |||
| 90 | } | |||
| 91 | } | |||
| 92 | else | |||
| 93 | { | |||
| 94 | fprintf(stderrstderr, "eigval[%d] = %g\n", i, eigval[i]); | |||
| 95 | w = 0; | |||
| 96 | } | |||
| 97 | } | |||
| 98 | fprintf(fp, "The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n", | |||
| 99 | S*RGAS((1.380658e-23)*(6.0221367e23))); | |||
| 100 | } | |||
| 101 | ||||
| 102 | static void calc_entropy_schlitter(FILE *fp, int n, int nskip, | |||
| 103 | real *eigval, real temp) | |||
| 104 | { | |||
| 105 | double dd, deter; | |||
| 106 | int *indx; | |||
| 107 | int i, j, k, m; | |||
| 108 | char buf[256]; | |||
| 109 | double hbar, kt, kteh, S; | |||
| 110 | ||||
| 111 | hbar = PLANCK1(6.6262e-34)/(2*M_PI3.14159265358979323846); | |||
| 112 | kt = BOLTZMANN(1.380658e-23)*temp; | |||
| 113 | kteh = kt*exp(2.0)/(hbar*hbar)*AMU(1.6605402e-27)*(NANO(1e-9)*NANO(1e-9)); | |||
| 114 | if (debug) | |||
| 115 | { | |||
| 116 | fprintf(debug, "n = %d, nskip = %d kteh = %g\n", n, nskip, kteh); | |||
| 117 | } | |||
| 118 | ||||
| 119 | deter = 0; | |||
| 120 | for (i = 0; (i < n-nskip); i++) | |||
| 121 | { | |||
| 122 | dd = 1+kteh*eigval[i]; | |||
| 123 | deter += log(dd); | |||
| 124 | } | |||
| 125 | S = 0.5*RGAS((1.380658e-23)*(6.0221367e23))*deter; | |||
| 126 | ||||
| 127 | fprintf(fp, "The Entropy due to the Schlitter formula is %g J/mol K\n", S); | |||
| 128 | } | |||
| 129 | ||||
| 130 | const char *proj_unit; | |||
| 131 | ||||
| 132 | static real tick_spacing(real range, int minticks) | |||
| 133 | { | |||
| 134 | real sp; | |||
| 135 | ||||
| 136 | if (range <= 0) | |||
| 137 | { | |||
| 138 | return 1.0; | |||
| 139 | } | |||
| 140 | ||||
| 141 | sp = 0.2*exp(log(10)*ceil(log(range)/log(10))); | |||
| 142 | while (range/sp < minticks-1) | |||
| 143 | { | |||
| 144 | sp = sp/2; | |||
| 145 | } | |||
| 146 | ||||
| 147 | return sp; | |||
| 148 | } | |||
| 149 | ||||
| 150 | static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph, | |||
| 151 | const char *title, const char *subtitle, | |||
| 152 | const char *xlabel, const char **ylabel, | |||
| 153 | int n, real *x, real **y, real ***sy, | |||
| 154 | real scale_x, gmx_bool bZero, gmx_bool bSplit, | |||
| 155 | const output_env_t oenv) | |||
| 156 | { | |||
| 157 | FILE *out; | |||
| 158 | int g, s, i; | |||
| 159 | real min, max, xsp, ysp; | |||
| 160 | ||||
| 161 | out = gmx_ffopen(file, "w"); | |||
| 162 | if (output_env_get_xvg_format(oenv) == exvgXMGRACE) | |||
| 163 | { | |||
| 164 | fprintf(out, "@ autoscale onread none\n"); | |||
| 165 | } | |||
| 166 | for (g = 0; g < ngraphs; g++) | |||
| 167 | { | |||
| 168 | if (y) | |||
| 169 | { | |||
| 170 | min = y[g][0]; | |||
| 171 | max = y[g][0]; | |||
| 172 | for (i = 0; i < n; i++) | |||
| 173 | { | |||
| 174 | if (y[g][i] < min) | |||
| 175 | { | |||
| 176 | min = y[g][i]; | |||
| 177 | } | |||
| 178 | if (y[g][i] > max) | |||
| 179 | { | |||
| 180 | max = y[g][i]; | |||
| 181 | } | |||
| 182 | } | |||
| 183 | } | |||
| 184 | else | |||
| 185 | { | |||
| 186 | min = sy[g][0][0]; | |||
| 187 | max = sy[g][0][0]; | |||
| 188 | for (s = 0; s < nsetspergraph; s++) | |||
| 189 | { | |||
| 190 | for (i = 0; i < n; i++) | |||
| 191 | { | |||
| 192 | if (sy[g][s][i] < min) | |||
| 193 | { | |||
| 194 | min = sy[g][s][i]; | |||
| 195 | } | |||
| 196 | if (sy[g][s][i] > max) | |||
| 197 | { | |||
| 198 | max = sy[g][s][i]; | |||
| 199 | } | |||
| 200 | } | |||
| 201 | } | |||
| 202 | } | |||
| 203 | if (bZero) | |||
| 204 | { | |||
| 205 | min = 0; | |||
| 206 | } | |||
| 207 | else | |||
| 208 | { | |||
| 209 | min = min-0.1*(max-min); | |||
| 210 | } | |||
| 211 | max = max+0.1*(max-min); | |||
| 212 | xsp = tick_spacing((x[n-1]-x[0])*scale_x, 4); | |||
| 213 | ysp = tick_spacing(max-min, 3); | |||
| 214 | if (output_env_get_print_xvgr_codes(oenv)) | |||
| 215 | { | |||
| 216 | fprintf(out, "@ with g%d\n@ g%d on\n", g, g); | |||
| 217 | if (g == 0) | |||
| 218 | { | |||
| 219 | fprintf(out, "@ title \"%s\"\n", title); | |||
| 220 | if (subtitle) | |||
| 221 | { | |||
| 222 | fprintf(out, "@ subtitle \"%s\"\n", subtitle); | |||
| 223 | } | |||
| 224 | } | |||
| 225 | if (g == ngraphs-1) | |||
| 226 | { | |||
| 227 | fprintf(out, "@ xaxis label \"%s\"\n", xlabel); | |||
| 228 | } | |||
| 229 | else | |||
| 230 | { | |||
| 231 | fprintf(out, "@ xaxis ticklabel off\n"); | |||
| 232 | } | |||
| 233 | if (n > 1) | |||
| 234 | { | |||
| 235 | fprintf(out, "@ world xmin %g\n", x[0]*scale_x); | |||
| 236 | fprintf(out, "@ world xmax %g\n", x[n-1]*scale_x); | |||
| 237 | fprintf(out, "@ world ymin %g\n", min); | |||
| 238 | fprintf(out, "@ world ymax %g\n", max); | |||
| 239 | } | |||
| 240 | fprintf(out, "@ view xmin 0.15\n"); | |||
| 241 | fprintf(out, "@ view xmax 0.85\n"); | |||
| 242 | fprintf(out, "@ view ymin %g\n", 0.15+(ngraphs-1-g)*0.7/ngraphs); | |||
| 243 | fprintf(out, "@ view ymax %g\n", 0.15+(ngraphs-g)*0.7/ngraphs); | |||
| 244 | fprintf(out, "@ yaxis label \"%s\"\n", ylabel[g]); | |||
| 245 | fprintf(out, "@ xaxis tick major %g\n", xsp); | |||
| 246 | fprintf(out, "@ xaxis tick minor %g\n", xsp/2); | |||
| 247 | fprintf(out, "@ xaxis ticklabel start type spec\n"); | |||
| 248 | fprintf(out, "@ xaxis ticklabel start %g\n", ceil(min/xsp)*xsp); | |||
| 249 | fprintf(out, "@ yaxis tick major %g\n", ysp); | |||
| 250 | fprintf(out, "@ yaxis tick minor %g\n", ysp/2); | |||
| 251 | fprintf(out, "@ yaxis ticklabel start type spec\n"); | |||
| 252 | fprintf(out, "@ yaxis ticklabel start %g\n", ceil(min/ysp)*ysp); | |||
| 253 | if ((min < 0) && (max > 0)) | |||
| 254 | { | |||
| 255 | fprintf(out, "@ zeroxaxis bar on\n"); | |||
| 256 | fprintf(out, "@ zeroxaxis bar linestyle 3\n"); | |||
| 257 | } | |||
| 258 | } | |||
| 259 | for (s = 0; s < nsetspergraph; s++) | |||
| 260 | { | |||
| 261 | for (i = 0; i < n; i++) | |||
| 262 | { | |||
| 263 | if (bSplit && i > 0 && abs(x[i]) < 1e-5) | |||
| 264 | { | |||
| 265 | if (output_env_get_print_xvgr_codes(oenv)) | |||
| 266 | { | |||
| 267 | fprintf(out, "&\n"); | |||
| 268 | } | |||
| 269 | } | |||
| 270 | fprintf(out, "%10.4f %10.5f\n", | |||
| 271 | x[i]*scale_x, y ? y[g][i] : sy[g][s][i]); | |||
| 272 | } | |||
| 273 | if (output_env_get_print_xvgr_codes(oenv)) | |||
| 274 | { | |||
| 275 | fprintf(out, "&\n"); | |||
| 276 | } | |||
| 277 | } | |||
| 278 | } | |||
| 279 | gmx_ffclose(out); | |||
| 280 | } | |||
| 281 | ||||
| 282 | static void | |||
| 283 | compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2, | |||
| 284 | real *eigval1, int neig1, real *eigval2, int neig2) | |||
| 285 | { | |||
| 286 | int n, nrow; | |||
| 287 | int i, j, k; | |||
| 288 | double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip; | |||
| 289 | ||||
| 290 | n = min(n1, n2)(((n1) < (n2)) ? (n1) : (n2) ); | |||
| 291 | ||||
| 292 | n = min(n, min(neig1, neig2))(((n) < ((((neig1) < (neig2)) ? (neig1) : (neig2) ))) ? (n) : ((((neig1) < (neig2)) ? (neig1) : (neig2) )) ); | |||
| 293 | fprintf(stdoutstdout, "Will compare the covariance matrices using %d dimensions\n", n); | |||
| 294 | ||||
| 295 | sum1 = 0; | |||
| 296 | for (i = 0; i < n; i++) | |||
| 297 | { | |||
| 298 | if (eigval1[i] < 0) | |||
| 299 | { | |||
| 300 | eigval1[i] = 0; | |||
| 301 | } | |||
| 302 | sum1 += eigval1[i]; | |||
| 303 | eigval1[i] = sqrt(eigval1[i]); | |||
| 304 | } | |||
| 305 | trace1 = sum1; | |||
| 306 | for (i = n; i < neig1; i++) | |||
| 307 | { | |||
| 308 | trace1 += eigval1[i]; | |||
| 309 | } | |||
| 310 | sum2 = 0; | |||
| 311 | for (i = 0; i < n; i++) | |||
| 312 | { | |||
| 313 | if (eigval2[i] < 0) | |||
| 314 | { | |||
| 315 | eigval2[i] = 0; | |||
| 316 | } | |||
| 317 | sum2 += eigval2[i]; | |||
| 318 | eigval2[i] = sqrt(eigval2[i]); | |||
| 319 | } | |||
| 320 | trace2 = sum2; | |||
| 321 | for (i = n; i < neig2; i++) | |||
| 322 | { | |||
| 323 | trace2 += eigval2[i]; | |||
| 324 | } | |||
| 325 | ||||
| 326 | fprintf(stdoutstdout, "Trace of the two matrices: %g and %g\n", sum1, sum2); | |||
| 327 | if (neig1 != n || neig2 != n) | |||
| 328 | { | |||
| 329 | fprintf(stdoutstdout, "this is %d%% and %d%% of the total trace\n", | |||
| 330 | (int)(100*sum1/trace1+0.5), (int)(100*sum2/trace2+0.5)); | |||
| 331 | } | |||
| 332 | fprintf(stdoutstdout, "Square root of the traces: %g and %g\n", | |||
| 333 | sqrt(sum1), sqrt(sum2)); | |||
| 334 | ||||
| 335 | sab = 0; | |||
| 336 | for (i = 0; i < n; i++) | |||
| 337 | { | |||
| 338 | tmp = 0; | |||
| 339 | for (j = 0; j < n; j++) | |||
| 340 | { | |||
| 341 | ip = 0; | |||
| 342 | for (k = 0; k < natoms; k++) | |||
| 343 | { | |||
| 344 | ip += iprod(eigvec1[i][k], eigvec2[j][k]); | |||
| 345 | } | |||
| 346 | tmp += eigval2[j]*ip*ip; | |||
| 347 | } | |||
| 348 | sab += eigval1[i]*tmp; | |||
| 349 | } | |||
| 350 | ||||
| 351 | samsb2 = sum1+sum2-2*sab; | |||
| 352 | if (samsb2 < 0) | |||
| 353 | { | |||
| 354 | samsb2 = 0; | |||
| 355 | } | |||
| 356 | ||||
| 357 | fprintf(stdoutstdout, "The overlap of the covariance matrices:\n"); | |||
| 358 | fprintf(stdoutstdout, " normalized: %.3f\n", 1-sqrt(samsb2/(sum1+sum2))); | |||
| 359 | tmp = 1-sab/sqrt(sum1*sum2); | |||
| 360 | if (tmp < 0) | |||
| 361 | { | |||
| 362 | tmp = 0; | |||
| 363 | } | |||
| 364 | fprintf(stdoutstdout, " shape: %.3f\n", 1-sqrt(tmp)); | |||
| 365 | } | |||
| 366 | ||||
| 367 | ||||
| 368 | static void inprod_matrix(const char *matfile, int natoms, | |||
| 369 | int nvec1, int *eignr1, rvec **eigvec1, | |||
| 370 | int nvec2, int *eignr2, rvec **eigvec2, | |||
| 371 | gmx_bool bSelect, int noutvec, int *outvec) | |||
| 372 | { | |||
| 373 | FILE *out; | |||
| 374 | real **mat; | |||
| 375 | int i, x1, y1, x, y, nlevels; | |||
| 376 | int nx, ny; | |||
| 377 | real inp, *t_x, *t_y, max; | |||
| 378 | t_rgb rlo, rhi; | |||
| 379 | ||||
| 380 | snew(t_y, nvec2)(t_y) = save_calloc("t_y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 380, (nvec2), sizeof(*(t_y))); | |||
| 381 | if (bSelect) | |||
| 382 | { | |||
| 383 | nx = noutvec; | |||
| 384 | ny = 0; | |||
| 385 | for (y1 = 0; y1 < nx; y1++) | |||
| 386 | { | |||
| 387 | if (outvec[y1] < nvec2) | |||
| 388 | { | |||
| 389 | t_y[ny] = eignr2[outvec[y1]]+1; | |||
| 390 | ny++; | |||
| 391 | } | |||
| 392 | } | |||
| 393 | } | |||
| 394 | else | |||
| 395 | { | |||
| 396 | nx = nvec1; | |||
| 397 | ny = nvec2; | |||
| 398 | for (y = 0; y < ny; y++) | |||
| 399 | { | |||
| 400 | t_y[y] = eignr2[y]+1; | |||
| 401 | } | |||
| 402 | } | |||
| 403 | ||||
| 404 | fprintf(stderrstderr, "Calculating inner-product matrix of %dx%d eigenvectors\n", | |||
| 405 | nx, nvec2); | |||
| 406 | ||||
| 407 | snew(mat, nx)(mat) = save_calloc("mat", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 407, (nx), sizeof(*(mat))); | |||
| 408 | snew(t_x, nx)(t_x) = save_calloc("t_x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 408, (nx), sizeof(*(t_x))); | |||
| 409 | max = 0; | |||
| 410 | for (x1 = 0; x1 < nx; x1++) | |||
| 411 | { | |||
| 412 | snew(mat[x1], ny)(mat[x1]) = save_calloc("mat[x1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 412, (ny), sizeof(*(mat[x1]))); | |||
| 413 | if (bSelect) | |||
| 414 | { | |||
| 415 | x = outvec[x1]; | |||
| 416 | } | |||
| 417 | else | |||
| 418 | { | |||
| 419 | x = x1; | |||
| 420 | } | |||
| 421 | t_x[x1] = eignr1[x]+1; | |||
| 422 | fprintf(stderrstderr, " %d", eignr1[x]+1); | |||
| 423 | for (y1 = 0; y1 < ny; y1++) | |||
| 424 | { | |||
| 425 | inp = 0; | |||
| 426 | if (bSelect) | |||
| 427 | { | |||
| 428 | while (outvec[y1] >= nvec2) | |||
| 429 | { | |||
| 430 | y1++; | |||
| 431 | } | |||
| 432 | y = outvec[y1]; | |||
| 433 | } | |||
| 434 | else | |||
| 435 | { | |||
| 436 | y = y1; | |||
| 437 | } | |||
| 438 | for (i = 0; i < natoms; i++) | |||
| 439 | { | |||
| 440 | inp += iprod(eigvec1[x][i], eigvec2[y][i]); | |||
| 441 | } | |||
| 442 | mat[x1][y1] = fabs(inp); | |||
| 443 | if (mat[x1][y1] > max) | |||
| 444 | { | |||
| 445 | max = mat[x1][y1]; | |||
| 446 | } | |||
| 447 | } | |||
| 448 | } | |||
| 449 | fprintf(stderrstderr, "\n"); | |||
| 450 | rlo.r = 1; rlo.g = 1; rlo.b = 1; | |||
| 451 | rhi.r = 0; rhi.g = 0; rhi.b = 0; | |||
| 452 | nlevels = 41; | |||
| 453 | out = gmx_ffopen(matfile, "w"); | |||
| 454 | write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2", | |||
| 455 | nx, ny, t_x, t_y, mat, 0.0, max, rlo, rhi, &nlevels); | |||
| 456 | gmx_ffclose(out); | |||
| 457 | } | |||
| 458 | ||||
| 459 | static void overlap(const char *outfile, int natoms, | |||
| 460 | rvec **eigvec1, | |||
| 461 | int nvec2, int *eignr2, rvec **eigvec2, | |||
| 462 | int noutvec, int *outvec, | |||
| 463 | const output_env_t oenv) | |||
| 464 | { | |||
| 465 | FILE *out; | |||
| 466 | int i, v, vec, x; | |||
| 467 | real overlap, inp; | |||
| 468 | ||||
| 469 | fprintf(stderrstderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n"); | |||
| 470 | for (i = 0; i < noutvec; i++) | |||
| 471 | { | |||
| 472 | fprintf(stderrstderr, "%d ", outvec[i]+1); | |||
| 473 | } | |||
| 474 | fprintf(stderrstderr, "\n"); | |||
| 475 | ||||
| 476 | out = xvgropen(outfile, "Subspace overlap", | |||
| 477 | "Eigenvectors of trajectory 2", "Overlap", oenv); | |||
| 478 | fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n", | |||
| 479 | noutvec); | |||
| 480 | overlap = 0; | |||
| 481 | for (x = 0; x < nvec2; x++) | |||
| 482 | { | |||
| 483 | for (v = 0; v < noutvec; v++) | |||
| 484 | { | |||
| 485 | vec = outvec[v]; | |||
| 486 | inp = 0; | |||
| 487 | for (i = 0; i < natoms; i++) | |||
| 488 | { | |||
| 489 | inp += iprod(eigvec1[vec][i], eigvec2[x][i]); | |||
| 490 | } | |||
| 491 | overlap += sqr(inp); | |||
| 492 | } | |||
| 493 | fprintf(out, "%5d %5.3f\n", eignr2[x]+1, overlap/noutvec); | |||
| 494 | } | |||
| 495 | ||||
| 496 | gmx_ffclose(out); | |||
| 497 | } | |||
| 498 | ||||
| 499 | static void project(const char *trajfile, t_topology *top, int ePBC, matrix topbox, | |||
| 500 | const char *projfile, const char *twodplotfile, | |||
| 501 | const char *threedplotfile, const char *filterfile, int skip, | |||
| 502 | const char *extremefile, gmx_bool bExtrAll, real extreme, | |||
| 503 | int nextr, t_atoms *atoms, int natoms, atom_id *index, | |||
| 504 | gmx_bool bFit, rvec *xref, int nfit, atom_id *ifit, real *w_rls, | |||
| 505 | real *sqrtm, rvec *xav, | |||
| 506 | int *eignr, rvec **eigvec, | |||
| 507 | int noutvec, int *outvec, gmx_bool bSplit, | |||
| 508 | const output_env_t oenv) | |||
| 509 | { | |||
| 510 | FILE *xvgrout = NULL((void*)0); | |||
| 511 | int nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame; | |||
| 512 | t_trxstatus *out = NULL((void*)0); | |||
| 513 | t_trxstatus *status; | |||
| 514 | int noutvec_extr, imin, imax; | |||
| 515 | real *pmin, *pmax; | |||
| 516 | atom_id *all_at; | |||
| 517 | matrix box; | |||
| 518 | rvec *xread, *x; | |||
| 519 | real t, inp, **inprod = NULL((void*)0), min = 0, max = 0; | |||
| 520 | char str[STRLEN4096], str2[STRLEN4096], **ylabel, *c; | |||
| 521 | real fact; | |||
| 522 | gmx_rmpbc_t gpbc = NULL((void*)0); | |||
| 523 | ||||
| 524 | snew(x, natoms)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 524, (natoms), sizeof(*(x))); | |||
| 525 | ||||
| 526 | if (bExtrAll) | |||
| ||||
| 527 | { | |||
| 528 | noutvec_extr = noutvec; | |||
| 529 | } | |||
| 530 | else | |||
| 531 | { | |||
| 532 | noutvec_extr = 1; | |||
| 533 | } | |||
| 534 | ||||
| 535 | ||||
| 536 | if (trajfile) | |||
| 537 | { | |||
| 538 | snew(inprod, noutvec+1)(inprod) = save_calloc("inprod", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 538, (noutvec+1), sizeof(*(inprod))); | |||
| 539 | ||||
| 540 | if (filterfile) | |||
| 541 | { | |||
| 542 | fprintf(stderrstderr, "Writing a filtered trajectory to %s using eigenvectors\n", | |||
| 543 | filterfile); | |||
| 544 | for (i = 0; i < noutvec; i++) | |||
| 545 | { | |||
| 546 | fprintf(stderrstderr, "%d ", outvec[i]+1); | |||
| 547 | } | |||
| 548 | fprintf(stderrstderr, "\n"); | |||
| 549 | out = open_trx(filterfile, "w"); | |||
| 550 | } | |||
| 551 | snew_size = 0; | |||
| 552 | nfr = 0; | |||
| 553 | nframes = 0; | |||
| 554 | nat = read_first_x(oenv, &status, trajfile, &t, &xread, box); | |||
| 555 | if (nat > atoms->nr) | |||
| 556 | { | |||
| 557 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 557, "the number of atoms in your trajectory (%d) is larger than the number of atoms in your structure file (%d)", nat, atoms->nr); | |||
| 558 | } | |||
| 559 | snew(all_at, nat)(all_at) = save_calloc("all_at", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 559, (nat), sizeof(*(all_at))); | |||
| 560 | ||||
| 561 | if (top) | |||
| 562 | { | |||
| 563 | gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat); | |||
| 564 | } | |||
| 565 | ||||
| 566 | for (i = 0; i < nat; i++) | |||
| 567 | { | |||
| 568 | all_at[i] = i; | |||
| 569 | } | |||
| 570 | do | |||
| 571 | { | |||
| 572 | if (nfr % skip == 0) | |||
| 573 | { | |||
| 574 | if (top) | |||
| 575 | { | |||
| 576 | gmx_rmpbc(gpbc, nat, box, xread); | |||
| 577 | } | |||
| 578 | if (nframes >= snew_size) | |||
| 579 | { | |||
| 580 | snew_size += 100; | |||
| 581 | for (i = 0; i < noutvec+1; i++) | |||
| 582 | { | |||
| 583 | srenew(inprod[i], snew_size)(inprod[i]) = save_realloc("inprod[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 583, (inprod[i]), (snew_size), sizeof(*(inprod[i]))); | |||
| 584 | } | |||
| 585 | } | |||
| 586 | inprod[noutvec][nframes] = t; | |||
| 587 | /* calculate x: a fitted struture of the selected atoms */ | |||
| 588 | if (bFit) | |||
| 589 | { | |||
| 590 | reset_x(nfit, ifit, nat, NULL((void*)0), xread, w_rls); | |||
| 591 | do_fit(nat, w_rls, xref, xread); | |||
| 592 | } | |||
| 593 | for (i = 0; i < natoms; i++) | |||
| 594 | { | |||
| 595 | copy_rvec(xread[index[i]], x[i]); | |||
| 596 | } | |||
| 597 | ||||
| 598 | for (v = 0; v < noutvec; v++) | |||
| 599 | { | |||
| 600 | vec = outvec[v]; | |||
| 601 | /* calculate (mass-weighted) projection */ | |||
| 602 | inp = 0; | |||
| 603 | for (i = 0; i < natoms; i++) | |||
| 604 | { | |||
| 605 | inp += (eigvec[vec][i][0]*(x[i][0]-xav[i][0])+ | |||
| 606 | eigvec[vec][i][1]*(x[i][1]-xav[i][1])+ | |||
| 607 | eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i]; | |||
| 608 | } | |||
| 609 | inprod[v][nframes] = inp; | |||
| 610 | } | |||
| 611 | if (filterfile) | |||
| 612 | { | |||
| 613 | for (i = 0; i < natoms; i++) | |||
| 614 | { | |||
| 615 | for (d = 0; d < DIM3; d++) | |||
| 616 | { | |||
| 617 | /* misuse xread for output */ | |||
| 618 | xread[index[i]][d] = xav[i][d]; | |||
| 619 | for (v = 0; v < noutvec; v++) | |||
| 620 | { | |||
| 621 | xread[index[i]][d] += | |||
| 622 | inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i]; | |||
| 623 | } | |||
| 624 | } | |||
| 625 | } | |||
| 626 | write_trx(out, natoms, index, atoms, 0, t, box, xread, NULL((void*)0), NULL((void*)0)); | |||
| 627 | } | |||
| 628 | nframes++; | |||
| 629 | } | |||
| 630 | nfr++; | |||
| 631 | } | |||
| 632 | while (read_next_x(oenv, status, &t, xread, box)); | |||
| 633 | close_trx(status); | |||
| 634 | sfree(x)save_free("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 634, (x)); | |||
| 635 | if (filterfile) | |||
| 636 | { | |||
| 637 | close_trx(out); | |||
| 638 | } | |||
| 639 | } | |||
| 640 | else | |||
| 641 | { | |||
| 642 | snew(xread, atoms->nr)(xread) = save_calloc("xread", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 642, (atoms->nr), sizeof(*(xread))); | |||
| 643 | } | |||
| 644 | ||||
| 645 | if (top) | |||
| 646 | { | |||
| 647 | gmx_rmpbc_done(gpbc); | |||
| 648 | } | |||
| 649 | ||||
| 650 | ||||
| 651 | if (projfile) | |||
| 652 | { | |||
| 653 | snew(ylabel, noutvec)(ylabel) = save_calloc("ylabel", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 653, (noutvec), sizeof(*(ylabel))); | |||
| 654 | for (v = 0; v < noutvec; v++) | |||
| 655 | { | |||
| 656 | sprintf(str, "vec %d", eignr[outvec[v]]+1); | |||
| 657 | ylabel[v] = strdup(str)(__extension__ (__builtin_constant_p (str) && ((size_t )(const void *)((str) + 1) - (size_t)(const void *)(str) == 1 ) ? (((const char *) (str))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (str) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, str, __len); __retval ; })) : __strdup (str))); | |||
| 658 | } | |||
| 659 | sprintf(str, "projection on eigenvectors (%s)", proj_unit); | |||
| 660 | write_xvgr_graphs(projfile, noutvec, 1, str, NULL((void*)0), output_env_get_xvgr_tlabel(oenv), | |||
| 661 | (const char **)ylabel, | |||
| 662 | nframes, inprod[noutvec], inprod, NULL((void*)0), | |||
| 663 | output_env_get_time_factor(oenv), FALSE0, bSplit, oenv); | |||
| 664 | } | |||
| 665 | ||||
| 666 | if (twodplotfile) | |||
| 667 | { | |||
| 668 | sprintf(str, "projection on eigenvector %d (%s)", | |||
| 669 | eignr[outvec[0]]+1, proj_unit); | |||
| 670 | sprintf(str2, "projection on eigenvector %d (%s)", | |||
| 671 | eignr[outvec[noutvec-1]]+1, proj_unit); | |||
| 672 | xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2, | |||
| 673 | oenv); | |||
| 674 | for (i = 0; i < nframes; i++) | |||
| 675 | { | |||
| 676 | if (bSplit && i > 0 && abs(inprod[noutvec][i]) < 1e-5) | |||
| 677 | { | |||
| 678 | fprintf(xvgrout, "&\n"); | |||
| 679 | } | |||
| 680 | fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]); | |||
| 681 | } | |||
| 682 | gmx_ffclose(xvgrout); | |||
| 683 | } | |||
| 684 | ||||
| 685 | if (threedplotfile) | |||
| 686 | { | |||
| 687 | t_atoms atoms; | |||
| 688 | rvec *x; | |||
| 689 | real *b = NULL((void*)0); | |||
| 690 | matrix box; | |||
| 691 | char *resnm, *atnm, pdbform[STRLEN4096]; | |||
| 692 | gmx_bool bPDB, b4D; | |||
| 693 | FILE *out; | |||
| 694 | ||||
| 695 | if (noutvec < 3) | |||
| 696 | { | |||
| 697 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 697, "You have selected less than 3 eigenvectors"); | |||
| 698 | } | |||
| 699 | ||||
| 700 | /* initialize */ | |||
| 701 | bPDB = fn2ftp(threedplotfile) == efPDB; | |||
| 702 | clear_mat(box); | |||
| 703 | box[XX0][XX0] = box[YY1][YY1] = box[ZZ2][ZZ2] = 1; | |||
| 704 | ||||
| 705 | b4D = bPDB && (noutvec >= 4); | |||
| 706 | if (b4D) | |||
| 707 | { | |||
| 708 | fprintf(stderrstderr, "You have selected four or more eigenvectors:\n" | |||
| 709 | "fourth eigenvector will be plotted " | |||
| 710 | "in bfactor field of pdb file\n"); | |||
| 711 | sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d", | |||
| 712 | eignr[outvec[0]]+1, eignr[outvec[1]]+1, | |||
| 713 | eignr[outvec[2]]+1, eignr[outvec[3]]+1); | |||
| 714 | } | |||
| 715 | else | |||
| 716 | { | |||
| 717 | sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d", | |||
| 718 | eignr[outvec[0]]+1, eignr[outvec[1]]+1, eignr[outvec[2]]+1); | |||
| 719 | } | |||
| 720 | init_t_atoms(&atoms, nframes, FALSE0); | |||
| 721 | snew(x, nframes)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 721, (nframes), sizeof(*(x))); | |||
| 722 | snew(b, nframes)(b) = save_calloc("b", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 722, (nframes), sizeof(*(b))); | |||
| 723 | atnm = strdup("C")(__extension__ (__builtin_constant_p ("C") && ((size_t )(const void *)(("C") + 1) - (size_t)(const void *)("C") == 1 ) ? (((const char *) ("C"))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen ("C") + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, "C", __len); __retval ; })) : __strdup ("C"))); | |||
| 724 | resnm = strdup("PRJ")(__extension__ (__builtin_constant_p ("PRJ") && ((size_t )(const void *)(("PRJ") + 1) - (size_t)(const void *)("PRJ") == 1) ? (((const char *) ("PRJ"))[0] == '\0' ? (char *) calloc ( (size_t) 1, (size_t) 1) : ({ size_t __len = strlen ("PRJ") + 1 ; char *__retval = (char *) malloc (__len); if (__retval != ( (void*)0)) __retval = (char *) memcpy (__retval, "PRJ", __len ); __retval; })) : __strdup ("PRJ"))); | |||
| 725 | ||||
| 726 | if (nframes > 10000) | |||
| 727 | { | |||
| 728 | fact = 10000.0/nframes; | |||
| 729 | } | |||
| 730 | else | |||
| 731 | { | |||
| 732 | fact = 1.0; | |||
| 733 | } | |||
| 734 | ||||
| 735 | for (i = 0; i < nframes; i++) | |||
| 736 | { | |||
| 737 | atoms.atomname[i] = &atnm; | |||
| 738 | atoms.atom[i].resind = i; | |||
| 739 | atoms.resinfo[i].name = &resnm; | |||
| 740 | atoms.resinfo[i].nr = ceil(i*fact); | |||
| 741 | atoms.resinfo[i].ic = ' '; | |||
| 742 | x[i][XX0] = inprod[0][i]; | |||
| 743 | x[i][YY1] = inprod[1][i]; | |||
| 744 | x[i][ZZ2] = inprod[2][i]; | |||
| 745 | if (b4D) | |||
| 746 | { | |||
| 747 | b[i] = inprod[3][i]; | |||
| 748 | } | |||
| 749 | } | |||
| 750 | if ( ( b4D || bSplit ) && bPDB) | |||
| ||||
| 751 | { | |||
| 752 | strcpy(pdbform, get_pdbformat()); | |||
| 753 | strcat(pdbform, "%8.4f%8.4f\n"); | |||
| 754 | ||||
| 755 | out = gmx_ffopen(threedplotfile, "w"); | |||
| 756 | fprintf(out, "HEADER %s\n", str); | |||
| 757 | if (b4D) | |||
| 758 | { | |||
| 759 | fprintf(out, "REMARK %s\n", "fourth dimension plotted as B-factor"); | |||
| 760 | } | |||
| 761 | j = 0; | |||
| 762 | for (i = 0; i < atoms.nr; i++) | |||
| 763 | { | |||
| 764 | if (j > 0 && bSplit && abs(inprod[noutvec][i]) < 1e-5) | |||
| 765 | { | |||
| 766 | fprintf(out, "TER\n"); | |||
| 767 | j = 0; | |||
| 768 | } | |||
| 769 | fprintf(out, pdbform, "ATOM", i+1, "C", "PRJ", ' ', j+1, | |||
| 770 | 10*x[i][XX0], 10*x[i][YY1], 10*x[i][ZZ2], 1.0, 10*b[i]); | |||
| 771 | if (j > 0) | |||
| 772 | { | |||
| 773 | fprintf(out, "CONECT%5d%5d\n", i, i+1); | |||
| 774 | } | |||
| 775 | j++; | |||
| 776 | } | |||
| 777 | fprintf(out, "TER\n"); | |||
| 778 | gmx_ffclose(out); | |||
| 779 | } | |||
| 780 | else | |||
| 781 | { | |||
| 782 | write_sto_conf(threedplotfile, str, &atoms, x, NULL((void*)0), ePBC, box); | |||
| 783 | } | |||
| 784 | free_t_atoms(&atoms, FALSE0); | |||
| 785 | } | |||
| 786 | ||||
| 787 | if (extremefile) | |||
| 788 | { | |||
| 789 | snew(pmin, noutvec_extr)(pmin) = save_calloc("pmin", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 789, (noutvec_extr), sizeof(*(pmin))); | |||
| 790 | snew(pmax, noutvec_extr)(pmax) = save_calloc("pmax", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 790, (noutvec_extr), sizeof(*(pmax))); | |||
| 791 | if (extreme == 0) | |||
| 792 | { | |||
| 793 | fprintf(stderrstderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum"); | |||
| 794 | fprintf(stderrstderr, | |||
| 795 | "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame"); | |||
| 796 | imin = 0; | |||
| 797 | imax = 0; | |||
| 798 | for (v = 0; v < noutvec_extr; v++) | |||
| 799 | { | |||
| 800 | for (i = 0; i < nframes; i++) | |||
| 801 | { | |||
| 802 | if (inprod[v][i] < inprod[v][imin]) | |||
| 803 | { | |||
| 804 | imin = i; | |||
| 805 | } | |||
| 806 | if (inprod[v][i] > inprod[v][imax]) | |||
| 807 | { | |||
| 808 | imax = i; | |||
| 809 | } | |||
| 810 | } | |||
| 811 | pmin[v] = inprod[v][imin]; | |||
| 812 | pmax[v] = inprod[v][imax]; | |||
| 813 | fprintf(stderrstderr, "%7d %10.6f %10d %10.6f %10d\n", | |||
| 814 | eignr[outvec[v]]+1, | |||
| 815 | pmin[v], imin, pmax[v], imax); | |||
| 816 | } | |||
| 817 | } | |||
| 818 | else | |||
| 819 | { | |||
| 820 | pmin[0] = -extreme; | |||
| 821 | pmax[0] = extreme; | |||
| 822 | } | |||
| 823 | /* build format string for filename: */ | |||
| 824 | strcpy(str, extremefile); /* copy filename */ | |||
| 825 | c = strrchr(str, '.'); /* find where extention begins */ | |||
| 826 | strcpy(str2, c); /* get extention */ | |||
| 827 | sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */ | |||
| 828 | for (v = 0; v < noutvec_extr; v++) | |||
| 829 | { | |||
| 830 | /* make filename using format string */ | |||
| 831 | if (noutvec_extr == 1) | |||
| 832 | { | |||
| 833 | strcpy(str2, extremefile); | |||
| 834 | } | |||
| 835 | else | |||
| 836 | { | |||
| 837 | sprintf(str2, str, eignr[outvec[v]]+1); | |||
| 838 | } | |||
| 839 | fprintf(stderrstderr, "Writing %d frames along eigenvector %d to %s\n", | |||
| 840 | nextr, outvec[v]+1, str2); | |||
| 841 | out = open_trx(str2, "w"); | |||
| 842 | for (frame = 0; frame < nextr; frame++) | |||
| 843 | { | |||
| 844 | if ((extreme == 0) && (nextr <= 3)) | |||
| 845 | { | |||
| 846 | for (i = 0; i < natoms; i++) | |||
| 847 | { | |||
| 848 | atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame; | |||
| 849 | } | |||
| 850 | } | |||
| 851 | for (i = 0; i < natoms; i++) | |||
| 852 | { | |||
| 853 | for (d = 0; d < DIM3; d++) | |||
| 854 | { | |||
| 855 | xread[index[i]][d] = | |||
| 856 | (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1) | |||
| 857 | *eigvec[outvec[v]][i][d]/sqrtm[i]); | |||
| 858 | } | |||
| 859 | } | |||
| 860 | write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, NULL((void*)0), NULL((void*)0)); | |||
| 861 | } | |||
| 862 | close_trx(out); | |||
| 863 | } | |||
| 864 | sfree(pmin)save_free("pmin", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 864, (pmin)); | |||
| 865 | sfree(pmax)save_free("pmax", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 865, (pmax)); | |||
| 866 | } | |||
| 867 | fprintf(stderrstderr, "\n"); | |||
| 868 | } | |||
| 869 | ||||
| 870 | static void components(const char *outfile, int natoms, | |||
| 871 | int *eignr, rvec **eigvec, | |||
| 872 | int noutvec, int *outvec, | |||
| 873 | const output_env_t oenv) | |||
| 874 | { | |||
| 875 | int g, s, v, i; | |||
| 876 | real *x, ***y; | |||
| 877 | char str[STRLEN4096], **ylabel; | |||
| 878 | ||||
| 879 | fprintf(stderrstderr, "Writing eigenvector components to %s\n", outfile); | |||
| 880 | ||||
| 881 | snew(ylabel, noutvec)(ylabel) = save_calloc("ylabel", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 881, (noutvec), sizeof(*(ylabel))); | |||
| 882 | snew(y, noutvec)(y) = save_calloc("y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 882, (noutvec), sizeof(*(y))); | |||
| 883 | snew(x, natoms)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 883, (natoms), sizeof(*(x))); | |||
| 884 | for (i = 0; i < natoms; i++) | |||
| 885 | { | |||
| 886 | x[i] = i+1; | |||
| 887 | } | |||
| 888 | for (g = 0; g < noutvec; g++) | |||
| 889 | { | |||
| 890 | v = outvec[g]; | |||
| 891 | sprintf(str, "vec %d", eignr[v]+1); | |||
| 892 | ylabel[g] = strdup(str)(__extension__ (__builtin_constant_p (str) && ((size_t )(const void *)((str) + 1) - (size_t)(const void *)(str) == 1 ) ? (((const char *) (str))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (str) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, str, __len); __retval ; })) : __strdup (str))); | |||
| 893 | snew(y[g], 4)(y[g]) = save_calloc("y[g]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 893, (4), sizeof(*(y[g]))); | |||
| 894 | for (s = 0; s < 4; s++) | |||
| 895 | { | |||
| 896 | snew(y[g][s], natoms)(y[g][s]) = save_calloc("y[g][s]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 896, (natoms), sizeof(*(y[g][s]))); | |||
| 897 | } | |||
| 898 | for (i = 0; i < natoms; i++) | |||
| 899 | { | |||
| 900 | y[g][0][i] = norm(eigvec[v][i]); | |||
| 901 | for (s = 0; s < 3; s++) | |||
| 902 | { | |||
| 903 | y[g][s+1][i] = eigvec[v][i][s]; | |||
| 904 | } | |||
| 905 | } | |||
| 906 | } | |||
| 907 | write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components", | |||
| 908 | "black: total, red: x, green: y, blue: z", | |||
| 909 | "Atom number", (const char **)ylabel, | |||
| 910 | natoms, x, NULL((void*)0), y, 1, FALSE0, FALSE0, oenv); | |||
| 911 | fprintf(stderrstderr, "\n"); | |||
| 912 | } | |||
| 913 | ||||
| 914 | static void rmsf(const char *outfile, int natoms, real *sqrtm, | |||
| 915 | int *eignr, rvec **eigvec, | |||
| 916 | int noutvec, int *outvec, | |||
| 917 | real *eigval, int neig, | |||
| 918 | const output_env_t oenv) | |||
| 919 | { | |||
| 920 | int nrow, g, v, i; | |||
| 921 | real *x, **y; | |||
| 922 | char str[STRLEN4096], **ylabel; | |||
| 923 | ||||
| 924 | for (i = 0; i < neig; i++) | |||
| 925 | { | |||
| 926 | if (eigval[i] < 0) | |||
| 927 | { | |||
| 928 | eigval[i] = 0; | |||
| 929 | } | |||
| 930 | } | |||
| 931 | ||||
| 932 | fprintf(stderrstderr, "Writing rmsf to %s\n", outfile); | |||
| 933 | ||||
| 934 | snew(ylabel, noutvec)(ylabel) = save_calloc("ylabel", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 934, (noutvec), sizeof(*(ylabel))); | |||
| 935 | snew(y, noutvec)(y) = save_calloc("y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 935, (noutvec), sizeof(*(y))); | |||
| 936 | snew(x, natoms)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 936, (natoms), sizeof(*(x))); | |||
| 937 | for (i = 0; i < natoms; i++) | |||
| 938 | { | |||
| 939 | x[i] = i+1; | |||
| 940 | } | |||
| 941 | for (g = 0; g < noutvec; g++) | |||
| 942 | { | |||
| 943 | v = outvec[g]; | |||
| 944 | if (eignr[v] >= neig) | |||
| 945 | { | |||
| 946 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 946, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig); | |||
| 947 | } | |||
| 948 | sprintf(str, "vec %d", eignr[v]+1); | |||
| 949 | ylabel[g] = strdup(str)(__extension__ (__builtin_constant_p (str) && ((size_t )(const void *)((str) + 1) - (size_t)(const void *)(str) == 1 ) ? (((const char *) (str))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (str) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, str, __len); __retval ; })) : __strdup (str))); | |||
| 950 | snew(y[g], natoms)(y[g]) = save_calloc("y[g]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 950, (natoms), sizeof(*(y[g]))); | |||
| 951 | for (i = 0; i < natoms; i++) | |||
| 952 | { | |||
| 953 | y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i]; | |||
| 954 | } | |||
| 955 | } | |||
| 956 | write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL((void*)0), | |||
| 957 | "Atom number", (const char **)ylabel, | |||
| 958 | natoms, x, y, NULL((void*)0), 1, TRUE1, FALSE0, oenv); | |||
| 959 | fprintf(stderrstderr, "\n"); | |||
| 960 | } | |||
| 961 | ||||
| 962 | int gmx_anaeig(int argc, char *argv[]) | |||
| 963 | { | |||
| 964 | static const char *desc[] = { | |||
| 965 | "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a", | |||
| 966 | "covariance matrix ([gmx-covar]) or of a Normal Modes analysis", | |||
| 967 | "([gmx-nmeig]).[PAR]", | |||
| 968 | ||||
| 969 | "When a trajectory is projected on eigenvectors, all structures are", | |||
| 970 | "fitted to the structure in the eigenvector file, if present, otherwise", | |||
| 971 | "to the structure in the structure file. When no run input file is", | |||
| 972 | "supplied, periodicity will not be taken into account. Most analyses", | |||
| 973 | "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when", | |||
| 974 | "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]", | |||
| 975 | ||||
| 976 | "[TT]-comp[tt]: plot the vector components per atom of eigenvectors", | |||
| 977 | "[TT]-first[tt] to [TT]-last[tt].[PAR]", | |||
| 978 | ||||
| 979 | "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors", | |||
| 980 | "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]", | |||
| 981 | ||||
| 982 | "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors", | |||
| 983 | "[TT]-first[tt] to [TT]-last[tt].", | |||
| 984 | "The projections of a trajectory on the eigenvectors of its", | |||
| 985 | "covariance matrix are called principal components (pc's).", | |||
| 986 | "It is often useful to check the cosine content of the pc's,", | |||
| 987 | "since the pc's of random diffusion are cosines with the number", | |||
| 988 | "of periods equal to half the pc index.", | |||
| 989 | "The cosine content of the pc's can be calculated with the program", | |||
| 990 | "[gmx-analyze].[PAR]", | |||
| 991 | ||||
| 992 | "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors", | |||
| 993 | "[TT]-first[tt] and [TT]-last[tt].[PAR]", | |||
| 994 | ||||
| 995 | "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first", | |||
| 996 | "three selected eigenvectors.[PAR]", | |||
| 997 | ||||
| 998 | "[TT]-filt[tt]: filter the trajectory to show only the motion along", | |||
| 999 | "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]", | |||
| 1000 | ||||
| 1001 | "[TT]-extr[tt]: calculate the two extreme projections along a trajectory", | |||
| 1002 | "on the average structure and interpolate [TT]-nframes[tt] frames", | |||
| 1003 | "between them, or set your own extremes with [TT]-max[tt]. The", | |||
| 1004 | "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and", | |||
| 1005 | "[TT]-last[tt] have been set explicitly, in which case all eigenvectors", | |||
| 1006 | "will be written to separate files. Chain identifiers will be added", | |||
| 1007 | "when writing a [TT].pdb[tt] file with two or three structures (you", | |||
| 1008 | "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]", | |||
| 1009 | ||||
| 1010 | " Overlap calculations between covariance analysis:[BR]", | |||
| 1011 | " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]", | |||
| 1012 | ||||
| 1013 | "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in", | |||
| 1014 | "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]", | |||
| 1015 | "in file [TT]-v[tt].[PAR]", | |||
| 1016 | ||||
| 1017 | "[TT]-inpr[tt]: calculate a matrix of inner-products between", | |||
| 1018 | "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors", | |||
| 1019 | "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]", | |||
| 1020 | "have been set explicitly.[PAR]", | |||
| 1021 | ||||
| 1022 | "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,", | |||
| 1023 | "a single number for the overlap between the covariance matrices is", | |||
| 1024 | "generated. The formulas are:[BR]", | |||
| 1025 | " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]", | |||
| 1026 | "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]", | |||
| 1027 | " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]", | |||
| 1028 | "where M1 and M2 are the two covariance matrices and tr is the trace", | |||
| 1029 | "of a matrix. The numbers are proportional to the overlap of the square", | |||
| 1030 | "root of the fluctuations. The normalized overlap is the most useful", | |||
| 1031 | "number, it is 1 for identical matrices and 0 when the sampled", | |||
| 1032 | "subspaces are orthogonal.[PAR]", | |||
| 1033 | "When the [TT]-entropy[tt] flag is given an entropy estimate will be", | |||
| 1034 | "computed based on the Quasiharmonic approach and based on", | |||
| 1035 | "Schlitter's formula." | |||
| 1036 | }; | |||
| 1037 | static int first = 1, last = -1, skip = 1, nextr = 2, nskip = 6; | |||
| 1038 | static real max = 0.0, temp = 298.15; | |||
| 1039 | static gmx_bool bSplit = FALSE0, bEntropy = FALSE0; | |||
| 1040 | t_pargs pa[] = { | |||
| 1041 | { "-first", FALSE0, etINT, {&first}, | |||
| 1042 | "First eigenvector for analysis (-1 is select)" }, | |||
| 1043 | { "-last", FALSE0, etINT, {&last}, | |||
| 1044 | "Last eigenvector for analysis (-1 is till the last)" }, | |||
| 1045 | { "-skip", FALSE0, etINT, {&skip}, | |||
| 1046 | "Only analyse every nr-th frame" }, | |||
| 1047 | { "-max", FALSE0, etREAL, {&max}, | |||
| 1048 | "Maximum for projection of the eigenvector on the average structure, " | |||
| 1049 | "max=0 gives the extremes" }, | |||
| 1050 | { "-nframes", FALSE0, etINT, {&nextr}, | |||
| 1051 | "Number of frames for the extremes output" }, | |||
| 1052 | { "-split", FALSE0, etBOOL, {&bSplit}, | |||
| 1053 | "Split eigenvector projections where time is zero" }, | |||
| 1054 | { "-entropy", FALSE0, etBOOL, {&bEntropy}, | |||
| 1055 | "Compute entropy according to the Quasiharmonic formula or Schlitter's method." }, | |||
| 1056 | { "-temp", FALSE0, etREAL, {&temp}, | |||
| 1057 | "Temperature for entropy calculations" }, | |||
| 1058 | { "-nevskip", FALSE0, etINT, {&nskip}, | |||
| 1059 | "Number of eigenvalues to skip when computing the entropy due to the quasi harmonic approximation. When you do a rotational and/or translational fit prior to the covariance analysis, you get 3 or 6 eigenvalues that are very close to zero, and which should not be taken into account when computing the entropy." } | |||
| 1060 | }; | |||
| 1061 | #define NPA((int)(sizeof(pa)/sizeof((pa)[0]))) asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))) | |||
| 1062 | ||||
| 1063 | FILE *out; | |||
| 1064 | int status, trjout; | |||
| 1065 | t_topology top; | |||
| 1066 | int ePBC = -1; | |||
| 1067 | t_atoms *atoms = NULL((void*)0); | |||
| 1068 | rvec *xtop, *xref1, *xref2, *xrefp = NULL((void*)0); | |||
| 1069 | gmx_bool bDMR1, bDMA1, bDMR2, bDMA2; | |||
| 1070 | int nvec1, nvec2, *eignr1 = NULL((void*)0), *eignr2 = NULL((void*)0); | |||
| 1071 | rvec *x, *xread, *xav1, *xav2, **eigvec1 = NULL((void*)0), **eigvec2 = NULL((void*)0); | |||
| 1072 | matrix topbox; | |||
| 1073 | real xid, totmass, *sqrtm, *w_rls, t, lambda; | |||
| 1074 | int natoms, step; | |||
| 1075 | char *grpname; | |||
| 1076 | const char *indexfile; | |||
| 1077 | char title[STRLEN4096]; | |||
| 1078 | int i, j, d; | |||
| 1079 | int nout, *iout, noutvec, *outvec, nfit; | |||
| 1080 | atom_id *index, *ifit; | |||
| 1081 | const char *VecFile, *Vec2File, *topfile; | |||
| 1082 | const char *EigFile, *Eig2File; | |||
| 1083 | const char *CompFile, *RmsfFile, *ProjOnVecFile; | |||
| 1084 | const char *TwoDPlotFile, *ThreeDPlotFile; | |||
| 1085 | const char *FilterFile, *ExtremeFile; | |||
| 1086 | const char *OverlapFile, *InpMatFile; | |||
| 1087 | gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj; | |||
| 1088 | gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D; | |||
| 1089 | real *eigval1 = NULL((void*)0), *eigval2 = NULL((void*)0); | |||
| 1090 | int neig1, neig2; | |||
| 1091 | double **xvgdata; | |||
| 1092 | output_env_t oenv; | |||
| 1093 | gmx_rmpbc_t gpbc; | |||
| 1094 | ||||
| 1095 | t_filenm fnm[] = { | |||
| 1096 | { efTRN, "-v", "eigenvec", ffREAD1<<1 }, | |||
| 1097 | { efTRN, "-v2", "eigenvec2", ffOPTRD(1<<1 | 1<<3) }, | |||
| 1098 | { efTRX, "-f", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, | |||
| 1099 | { efTPS, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, | |||
| 1100 | { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, | |||
| 1101 | { efXVG, "-eig", "eigenval", ffOPTRD(1<<1 | 1<<3) }, | |||
| 1102 | { efXVG, "-eig2", "eigenval2", ffOPTRD(1<<1 | 1<<3) }, | |||
| 1103 | { efXVG, "-comp", "eigcomp", ffOPTWR(1<<2| 1<<3) }, | |||
| 1104 | { efXVG, "-rmsf", "eigrmsf", ffOPTWR(1<<2| 1<<3) }, | |||
| 1105 | { efXVG, "-proj", "proj", ffOPTWR(1<<2| 1<<3) }, | |||
| 1106 | { efXVG, "-2d", "2dproj", ffOPTWR(1<<2| 1<<3) }, | |||
| 1107 | { efSTO, "-3d", "3dproj.pdb", ffOPTWR(1<<2| 1<<3) }, | |||
| 1108 | { efTRX, "-filt", "filtered", ffOPTWR(1<<2| 1<<3) }, | |||
| 1109 | { efTRX, "-extr", "extreme.pdb", ffOPTWR(1<<2| 1<<3) }, | |||
| 1110 | { efXVG, "-over", "overlap", ffOPTWR(1<<2| 1<<3) }, | |||
| 1111 | { efXPM, "-inpr", "inprod", ffOPTWR(1<<2| 1<<3) } | |||
| 1112 | }; | |||
| 1113 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) | |||
| 1114 | ||||
| 1115 | if (!parse_common_args(&argc, argv, | |||
| 1116 | PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_TIME_UNIT(1<<15) | PCA_CAN_VIEW(1<<5) | PCA_BE_NICE(1<<13), | |||
| 1117 | 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)) | |||
| 1118 | { | |||
| 1119 | return 0; | |||
| 1120 | } | |||
| 1121 | ||||
| 1122 | indexfile = ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1123 | ||||
| 1124 | VecFile = opt2fn("-v", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1125 | Vec2File = opt2fn_null("-v2", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1126 | topfile = ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1127 | EigFile = opt2fn_null("-eig", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1128 | Eig2File = opt2fn_null("-eig2", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1129 | CompFile = opt2fn_null("-comp", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1130 | RmsfFile = opt2fn_null("-rmsf", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1131 | ProjOnVecFile = opt2fn_null("-proj", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1132 | TwoDPlotFile = opt2fn_null("-2d", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1133 | ThreeDPlotFile = opt2fn_null("-3d", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1134 | FilterFile = opt2fn_null("-filt", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1135 | ExtremeFile = opt2fn_null("-extr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1136 | OverlapFile = opt2fn_null("-over", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1137 | InpMatFile = ftp2fn_null(efXPM, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1138 | ||||
| 1139 | bTop = fn2bTPX(topfile); | |||
| 1140 | bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile | |||
| 1141 | || FilterFile || ExtremeFile; | |||
| 1142 | bFirstLastSet = | |||
| 1143 | opt2parg_bSet("-first", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa) && opt2parg_bSet("-last", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa); | |||
| 1144 | bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile || | |||
| 1145 | OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet); | |||
| 1146 | bVec2 = Vec2File || OverlapFile || InpMatFile; | |||
| 1147 | bM = RmsfFile || bProj; | |||
| 1148 | bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0)) | |||
| 1149 | || TwoDPlotFile || ThreeDPlotFile; | |||
| 1150 | bIndex = bM || bProj; | |||
| 1151 | bTPS = ftp2bSet(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || bM || bTraj || | |||
| 1152 | FilterFile || (bIndex && indexfile); | |||
| 1153 | bCompare = Vec2File || Eig2File; | |||
| 1154 | bPDB3D = fn2ftp(ThreeDPlotFile) == efPDB; | |||
| 1155 | ||||
| 1156 | read_eigenvectors(VecFile, &natoms, &bFit1, | |||
| 1157 | &xref1, &bDMR1, &xav1, &bDMA1, | |||
| 1158 | &nvec1, &eignr1, &eigvec1, &eigval1); | |||
| 1159 | neig1 = DIM3*natoms; | |||
| 1160 | ||||
| 1161 | /* Overwrite eigenvalues from separate files if the user provides them */ | |||
| 1162 | if (EigFile != NULL((void*)0)) | |||
| 1163 | { | |||
| 1164 | int neig_tmp = read_xvg(EigFile, &xvgdata, &i); | |||
| 1165 | if (neig_tmp != neig1) | |||
| 1166 | { | |||
| 1167 | fprintf(stderrstderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms); | |||
| 1168 | } | |||
| 1169 | neig1 = neig_tmp; | |||
| 1170 | srenew(eigval1, neig1)(eigval1) = save_realloc("eigval1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1170, (eigval1), (neig1), sizeof(*(eigval1))); | |||
| 1171 | for (j = 0; j < neig1; j++) | |||
| 1172 | { | |||
| 1173 | real tmp = eigval1[j]; | |||
| 1174 | eigval1[j] = xvgdata[1][j]; | |||
| 1175 | if (debug && (eigval1[j] != tmp)) | |||
| 1176 | { | |||
| 1177 | fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n", | |||
| 1178 | j, tmp, eigval1[j]); | |||
| 1179 | } | |||
| 1180 | } | |||
| 1181 | for (j = 0; j < i; j++) | |||
| 1182 | { | |||
| 1183 | sfree(xvgdata[j])save_free("xvgdata[j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1183, (xvgdata[j])); | |||
| 1184 | } | |||
| 1185 | sfree(xvgdata)save_free("xvgdata", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1185, (xvgdata)); | |||
| 1186 | fprintf(stderrstderr, "Read %d eigenvalues from %s\n", neig1, EigFile); | |||
| 1187 | } | |||
| 1188 | ||||
| 1189 | if (bEntropy) | |||
| 1190 | { | |||
| 1191 | if (bDMA1) | |||
| 1192 | { | |||
| 1193 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1193, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting"); | |||
| 1194 | } | |||
| 1195 | calc_entropy_qh(stdoutstdout, neig1, eigval1, temp, nskip); | |||
| 1196 | calc_entropy_schlitter(stdoutstdout, neig1, nskip, eigval1, temp); | |||
| 1197 | } | |||
| 1198 | ||||
| 1199 | if (bVec2) | |||
| 1200 | { | |||
| 1201 | if (!Vec2File) | |||
| 1202 | { | |||
| 1203 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1203, "Need a second eigenvector file to do this analysis."); | |||
| 1204 | } | |||
| 1205 | read_eigenvectors(Vec2File, &neig2, &bFit2, | |||
| 1206 | &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2); | |||
| 1207 | ||||
| 1208 | neig2 = DIM3*neig2; | |||
| 1209 | if (neig2 != neig1) | |||
| 1210 | { | |||
| 1211 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1211, "Dimensions in the eigenvector files don't match"); | |||
| 1212 | } | |||
| 1213 | } | |||
| 1214 | ||||
| 1215 | if (Eig2File != NULL((void*)0)) | |||
| 1216 | { | |||
| 1217 | neig2 = read_xvg(Eig2File, &xvgdata, &i); | |||
| 1218 | srenew(eigval2, neig2)(eigval2) = save_realloc("eigval2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1218, (eigval2), (neig2), sizeof(*(eigval2))); | |||
| 1219 | for (j = 0; j < neig2; j++) | |||
| 1220 | { | |||
| 1221 | eigval2[j] = xvgdata[1][j]; | |||
| 1222 | } | |||
| 1223 | for (j = 0; j < i; j++) | |||
| 1224 | { | |||
| 1225 | sfree(xvgdata[j])save_free("xvgdata[j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1225, (xvgdata[j])); | |||
| 1226 | } | |||
| 1227 | sfree(xvgdata)save_free("xvgdata", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1227, (xvgdata)); | |||
| 1228 | fprintf(stderrstderr, "Read %d eigenvalues from %s\n", neig2, Eig2File); | |||
| 1229 | } | |||
| 1230 | ||||
| 1231 | ||||
| 1232 | if ((!bFit1 || xref1) && !bDMR1 && !bDMA1) | |||
| 1233 | { | |||
| 1234 | bM = FALSE0; | |||
| 1235 | } | |||
| 1236 | if ((xref1 == NULL((void*)0)) && (bM || bTraj)) | |||
| 1237 | { | |||
| 1238 | bTPS = TRUE1; | |||
| 1239 | } | |||
| 1240 | ||||
| 1241 | xtop = NULL((void*)0); | |||
| 1242 | nfit = 0; | |||
| 1243 | ifit = NULL((void*)0); | |||
| 1244 | w_rls = NULL((void*)0); | |||
| 1245 | ||||
| 1246 | if (!bTPS) | |||
| 1247 | { | |||
| 1248 | bTop = FALSE0; | |||
| 1249 | } | |||
| 1250 | else | |||
| 1251 | { | |||
| 1252 | bTop = read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
| 1253 | title, &top, &ePBC, &xtop, NULL((void*)0), topbox, bM); | |||
| 1254 | atoms = &top.atoms; | |||
| 1255 | gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr); | |||
| 1256 | gmx_rmpbc(gpbc, atoms->nr, topbox, xtop); | |||
| 1257 | /* Fitting is only required for the projection */ | |||
| 1258 | if (bProj && bFit1) | |||
| 1259 | { | |||
| 1260 | if (xref1 == NULL((void*)0)) | |||
| 1261 | { | |||
| 1262 | printf("\nNote: the structure in %s should be the same\n" | |||
| 1263 | " as the one used for the fit in g_covar\n", topfile); | |||
| 1264 | } | |||
| 1265 | printf("\nSelect the index group that was used for the least squares fit in g_covar\n"); | |||
| 1266 | get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname); | |||
| 1267 | ||||
| 1268 | snew(w_rls, atoms->nr)(w_rls) = save_calloc("w_rls", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1268, (atoms->nr), sizeof(*(w_rls))); | |||
| 1269 | for (i = 0; (i < nfit); i++) | |||
| 1270 | { | |||
| 1271 | if (bDMR1) | |||
| 1272 | { | |||
| 1273 | w_rls[ifit[i]] = atoms->atom[ifit[i]].m; | |||
| 1274 | } | |||
| 1275 | else | |||
| 1276 | { | |||
| 1277 | w_rls[ifit[i]] = 1.0; | |||
| 1278 | } | |||
| 1279 | } | |||
| 1280 | ||||
| 1281 | snew(xrefp, atoms->nr)(xrefp) = save_calloc("xrefp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1281, (atoms->nr), sizeof(*(xrefp))); | |||
| 1282 | if (xref1 != NULL((void*)0)) | |||
| 1283 | { | |||
| 1284 | /* Safety check between selected fit-group and reference structure read from the eigenvector file */ | |||
| 1285 | if (natoms != nfit) | |||
| 1286 | { | |||
| 1287 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1287, "you selected a group with %d elements instead of %d, your selection does not fit the reference structure in the eigenvector file.", nfit, natoms); | |||
| 1288 | } | |||
| 1289 | for (i = 0; (i < nfit); i++) | |||
| 1290 | { | |||
| 1291 | copy_rvec(xref1[i], xrefp[ifit[i]]); | |||
| 1292 | } | |||
| 1293 | } | |||
| 1294 | else | |||
| 1295 | { | |||
| 1296 | /* The top coordinates are the fitting reference */ | |||
| 1297 | for (i = 0; (i < nfit); i++) | |||
| 1298 | { | |||
| 1299 | copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]); | |||
| 1300 | } | |||
| 1301 | reset_x(nfit, ifit, atoms->nr, NULL((void*)0), xrefp, w_rls); | |||
| 1302 | } | |||
| 1303 | } | |||
| 1304 | gmx_rmpbc_done(gpbc); | |||
| 1305 | } | |||
| 1306 | ||||
| 1307 | if (bIndex) | |||
| 1308 | { | |||
| 1309 | printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms); | |||
| 1310 | get_index(atoms, indexfile, 1, &i, &index, &grpname); | |||
| 1311 | if (i != natoms) | |||
| 1312 | { | |||
| 1313 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1313, "you selected a group with %d elements instead of %d", i, natoms); | |||
| 1314 | } | |||
| 1315 | printf("\n"); | |||
| 1316 | } | |||
| 1317 | ||||
| 1318 | snew(sqrtm, natoms)(sqrtm) = save_calloc("sqrtm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1318, (natoms), sizeof(*(sqrtm))); | |||
| 1319 | if (bM && bDMA1) | |||
| 1320 | { | |||
| 1321 | proj_unit = "u\\S1/2\\Nnm"; | |||
| 1322 | for (i = 0; (i < natoms); i++) | |||
| 1323 | { | |||
| 1324 | sqrtm[i] = sqrt(atoms->atom[index[i]].m); | |||
| 1325 | } | |||
| 1326 | } | |||
| 1327 | else | |||
| 1328 | { | |||
| 1329 | proj_unit = "nm"; | |||
| 1330 | for (i = 0; (i < natoms); i++) | |||
| 1331 | { | |||
| 1332 | sqrtm[i] = 1.0; | |||
| 1333 | } | |||
| 1334 | } | |||
| 1335 | ||||
| 1336 | if (bVec2) | |||
| 1337 | { | |||
| 1338 | t = 0; | |||
| 1339 | totmass = 0; | |||
| 1340 | for (i = 0; (i < natoms); i++) | |||
| 1341 | { | |||
| 1342 | for (d = 0; (d < DIM3); d++) | |||
| 1343 | { | |||
| 1344 | t += sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]); | |||
| 1345 | totmass += sqr(sqrtm[i]); | |||
| 1346 | } | |||
| 1347 | } | |||
| 1348 | fprintf(stdoutstdout, "RMSD (without fit) between the two average structures:" | |||
| 1349 | " %.3f (nm)\n\n", sqrt(t/totmass)); | |||
| 1350 | } | |||
| 1351 | ||||
| 1352 | if (last == -1) | |||
| 1353 | { | |||
| 1354 | last = natoms*DIM3; | |||
| 1355 | } | |||
| 1356 | if (first > -1) | |||
| 1357 | { | |||
| 1358 | if (bFirstToLast) | |||
| 1359 | { | |||
| 1360 | /* make an index from first to last */ | |||
| 1361 | nout = last-first+1; | |||
| 1362 | snew(iout, nout)(iout) = save_calloc("iout", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1362, (nout), sizeof(*(iout))); | |||
| 1363 | for (i = 0; i < nout; i++) | |||
| 1364 | { | |||
| 1365 | iout[i] = first-1+i; | |||
| 1366 | } | |||
| 1367 | } | |||
| 1368 | else if (ThreeDPlotFile) | |||
| 1369 | { | |||
| 1370 | /* make an index of first+(0,1,2) and last */ | |||
| 1371 | nout = bPDB3D ? 4 : 3; | |||
| 1372 | nout = min(last-first+1, nout)(((last-first+1) < (nout)) ? (last-first+1) : (nout) ); | |||
| 1373 | snew(iout, nout)(iout) = save_calloc("iout", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1373, (nout), sizeof(*(iout))); | |||
| 1374 | iout[0] = first-1; | |||
| 1375 | iout[1] = first; | |||
| 1376 | if (nout > 3) | |||
| 1377 | { | |||
| 1378 | iout[2] = first+1; | |||
| 1379 | } | |||
| 1380 | iout[nout-1] = last-1; | |||
| 1381 | } | |||
| 1382 | else | |||
| 1383 | { | |||
| 1384 | /* make an index of first and last */ | |||
| 1385 | nout = 2; | |||
| 1386 | snew(iout, nout)(iout) = save_calloc("iout", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1386, (nout), sizeof(*(iout))); | |||
| 1387 | iout[0] = first-1; | |||
| 1388 | iout[1] = last-1; | |||
| 1389 | } | |||
| 1390 | } | |||
| 1391 | else | |||
| 1392 | { | |||
| 1393 | printf("Select eigenvectors for output, end your selection with 0\n"); | |||
| 1394 | nout = -1; | |||
| 1395 | iout = NULL((void*)0); | |||
| 1396 | ||||
| 1397 | do | |||
| 1398 | { | |||
| 1399 | nout++; | |||
| 1400 | srenew(iout, nout+1)(iout) = save_realloc("iout", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1400, (iout), (nout+1), sizeof(*(iout))); | |||
| 1401 | if (1 != scanf("%d", &iout[nout])) | |||
| 1402 | { | |||
| 1403 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1403, "Error reading user input"); | |||
| 1404 | } | |||
| 1405 | iout[nout]--; | |||
| 1406 | } | |||
| 1407 | while (iout[nout] >= 0); | |||
| 1408 | ||||
| 1409 | printf("\n"); | |||
| 1410 | } | |||
| 1411 | /* make an index of the eigenvectors which are present */ | |||
| 1412 | snew(outvec, nout)(outvec) = save_calloc("outvec", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1412, (nout), sizeof(*(outvec))); | |||
| 1413 | noutvec = 0; | |||
| 1414 | for (i = 0; i < nout; i++) | |||
| 1415 | { | |||
| 1416 | j = 0; | |||
| 1417 | while ((j < nvec1) && (eignr1[j] != iout[i])) | |||
| 1418 | { | |||
| 1419 | j++; | |||
| 1420 | } | |||
| 1421 | if ((j < nvec1) && (eignr1[j] == iout[i])) | |||
| 1422 | { | |||
| 1423 | outvec[noutvec] = j; | |||
| 1424 | noutvec++; | |||
| 1425 | } | |||
| 1426 | } | |||
| 1427 | fprintf(stderrstderr, "%d eigenvectors selected for output", noutvec); | |||
| 1428 | if (noutvec <= 100) | |||
| 1429 | { | |||
| 1430 | fprintf(stderrstderr, ":"); | |||
| 1431 | for (j = 0; j < noutvec; j++) | |||
| 1432 | { | |||
| 1433 | fprintf(stderrstderr, " %d", eignr1[outvec[j]]+1); | |||
| 1434 | } | |||
| 1435 | } | |||
| 1436 | fprintf(stderrstderr, "\n"); | |||
| 1437 | ||||
| 1438 | if (CompFile) | |||
| 1439 | { | |||
| 1440 | components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv); | |||
| 1441 | } | |||
| 1442 | ||||
| 1443 | if (RmsfFile) | |||
| 1444 | { | |||
| 1445 | rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1, | |||
| 1446 | neig1, oenv); | |||
| 1447 | } | |||
| 1448 | ||||
| 1449 | if (bProj) | |||
| 1450 | { | |||
| 1451 | project(bTraj ? opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) : NULL((void*)0), | |||
| 1452 | bTop ? &top : NULL((void*)0), ePBC, topbox, | |||
| 1453 | ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip, | |||
| 1454 | ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index, | |||
| 1455 | bFit1, xrefp, nfit, ifit, w_rls, | |||
| 1456 | sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit, | |||
| 1457 | oenv); | |||
| 1458 | } | |||
| 1459 | ||||
| 1460 | if (OverlapFile) | |||
| 1461 | { | |||
| 1462 | overlap(OverlapFile, natoms, | |||
| 1463 | eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv); | |||
| 1464 | } | |||
| 1465 | ||||
| 1466 | if (InpMatFile) | |||
| 1467 | { | |||
| 1468 | inprod_matrix(InpMatFile, natoms, | |||
| 1469 | nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2, | |||
| 1470 | bFirstLastSet, noutvec, outvec); | |||
| 1471 | } | |||
| 1472 | ||||
| 1473 | if (bCompare) | |||
| 1474 | { | |||
| 1475 | compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2); | |||
| 1476 | } | |||
| 1477 | ||||
| 1478 | ||||
| 1479 | if (!CompFile && !bProj && !OverlapFile && !InpMatFile && | |||
| 1480 | !bCompare && !bEntropy) | |||
| 1481 | { | |||
| 1482 | fprintf(stderrstderr, "\nIf you want some output," | |||
| 1483 | " set one (or two or ...) of the output file options\n"); | |||
| 1484 | } | |||
| 1485 | ||||
| 1486 | ||||
| 1487 | view_all(oenv, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 1488 | ||||
| 1489 | return 0; | |||
| 1490 | } |