| File: | gromacs/gmxana/gmx_make_edi.c |
| Location: | line 871, column 5 |
| Description: | Value stored to 'bTop' is never read |
| 1 | /* |
| 2 | * This file is part of the GROMACS molecular simulation package. |
| 3 | * |
| 4 | * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by |
| 5 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, |
| 6 | * and including many others, as listed in the AUTHORS file in the |
| 7 | * top-level source directory and at http://www.gromacs.org. |
| 8 | * |
| 9 | * GROMACS is free software; you can redistribute it and/or |
| 10 | * modify it under the terms of the GNU Lesser General Public License |
| 11 | * as published by the Free Software Foundation; either version 2.1 |
| 12 | * of the License, or (at your option) any later version. |
| 13 | * |
| 14 | * GROMACS is distributed in the hope that it will be useful, |
| 15 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 16 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| 17 | * Lesser General Public License for more details. |
| 18 | * |
| 19 | * You should have received a copy of the GNU Lesser General Public |
| 20 | * License along with GROMACS; if not, see |
| 21 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, |
| 22 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |
| 23 | * |
| 24 | * If you want to redistribute modifications to GROMACS, please |
| 25 | * consider that scientific software is very special. Version |
| 26 | * control is crucial - bugs must be traceable. We will be happy to |
| 27 | * consider code for inclusion in the official distribution, but |
| 28 | * derived work must not be called official GROMACS. Details are found |
| 29 | * in the README & COPYING files - if they are missing, get the |
| 30 | * official version at http://www.gromacs.org. |
| 31 | * |
| 32 | * To help us fund GROMACS development, we humbly ask that you cite |
| 33 | * the research papers on the package. Check out http://www.gromacs.org. |
| 34 | */ |
| 35 | /* The make_edi program was generously contributed by Oliver Lange, based |
| 36 | * on the code from g_anaeig. You can reach him as olange@gwdg.de. He |
| 37 | * probably also holds copyright to the following code. |
| 38 | */ |
| 39 | #ifdef HAVE_CONFIG_H1 |
| 40 | #include <config.h> |
| 41 | #endif |
| 42 | |
| 43 | #include <math.h> |
| 44 | #include <stdlib.h> |
| 45 | #include <ctype.h> |
| 46 | #include <string.h> |
| 47 | #include "readinp.h" |
| 48 | #include "gromacs/commandline/pargs.h" |
| 49 | #include "typedefs.h" |
| 50 | #include "gromacs/utility/smalloc.h" |
| 51 | #include "macros.h" |
| 52 | #include "gromacs/utility/fatalerror.h" |
| 53 | #include "gromacs/math/vec.h" |
| 54 | #include "pbc.h" |
| 55 | #include "gromacs/utility/futil.h" |
| 56 | #include "gromacs/fileio/pdbio.h" |
| 57 | #include "gromacs/fileio/confio.h" |
| 58 | #include "gromacs/fileio/tpxio.h" |
| 59 | #include "gromacs/fileio/matio.h" |
| 60 | #include "mshift.h" |
| 61 | #include "gromacs/fileio/xvgr.h" |
| 62 | #include "rmpbc.h" |
| 63 | #include "txtdump.h" |
| 64 | #include "eigio.h" |
| 65 | #include "index.h" |
| 66 | #include "gromacs/utility/cstringutil.h" |
| 67 | |
| 68 | typedef struct |
| 69 | { |
| 70 | real deltaF0; |
| 71 | gmx_bool bHarmonic; |
| 72 | gmx_bool bConstForce; /* Do constant force flooding instead of |
| 73 | evaluating a flooding potential */ |
| 74 | real tau; |
| 75 | real deltaF; |
| 76 | real kT; |
| 77 | real constEfl; |
| 78 | real alpha2; |
| 79 | } t_edflood; |
| 80 | |
| 81 | |
| 82 | /* This type is for the average, reference, target, and origin structure */ |
| 83 | typedef struct edix |
| 84 | { |
| 85 | int nr; /* number of atoms this structure contains */ |
| 86 | int *anrs; /* atom index numbers */ |
| 87 | rvec *x; /* positions */ |
| 88 | real *sqrtm; /* sqrt of the masses used for mass- |
| 89 | * weighting of analysis */ |
| 90 | } t_edix; |
| 91 | |
| 92 | |
| 93 | typedef struct edipar |
| 94 | { |
| 95 | int nini; /* total Nr of atoms */ |
| 96 | gmx_bool fitmas; /* true if trans fit with cm */ |
| 97 | gmx_bool pcamas; /* true if mass-weighted PCA */ |
| 98 | int presteps; /* number of steps to run without any |
| 99 | * perturbations ... just monitoring */ |
| 100 | int outfrq; /* freq (in steps) of writing to edo */ |
| 101 | int maxedsteps; /* max nr of steps per cycle */ |
| 102 | struct edix sref; /* reference positions, to these fitting |
| 103 | * will be done */ |
| 104 | struct edix sav; /* average positions */ |
| 105 | struct edix star; /* target positions */ |
| 106 | struct edix sori; /* origin positions */ |
| 107 | real slope; /* minimal slope in acceptance radexp */ |
| 108 | int ned; /* Nr of atoms in essdyn buffer */ |
| 109 | t_edflood flood; /* parameters especially for flooding */ |
| 110 | } t_edipar; |
| 111 | |
| 112 | |
| 113 | |
| 114 | void make_t_edx(struct edix *edx, int natoms, rvec *pos, atom_id index[]) |
| 115 | { |
| 116 | edx->nr = natoms; |
| 117 | edx->anrs = index; |
| 118 | edx->x = pos; |
| 119 | } |
| 120 | |
| 121 | void write_t_edx(FILE *fp, struct edix edx, const char *comment) |
| 122 | { |
| 123 | /*here we copy only the pointers into the t_edx struct |
| 124 | no data is copied and edx.box is ignored */ |
| 125 | int i; |
| 126 | fprintf(fp, "#%s \n %d \n", comment, edx.nr); |
| 127 | for (i = 0; i < edx.nr; i++) |
| 128 | { |
| 129 | fprintf(fp, "%d %f %f %f\n", (edx.anrs)[i]+1, (edx.x)[i][XX0], (edx.x)[i][YY1], (edx.x)[i][ZZ2]); |
| 130 | } |
| 131 | } |
| 132 | |
| 133 | int sscan_list(int *list[], const char *str, const char *listname) |
| 134 | { |
| 135 | /*this routine scans a string of the form 1,3-6,9 and returns the |
| 136 | selected numbers (in this case 1 3 4 5 6 9) in NULL-terminated array of integers. |
| 137 | memory for this list will be allocated in this routine -- sscan_list expects *list to |
| 138 | be a NULL-Pointer |
| 139 | |
| 140 | listname is a string used in the errormessage*/ |
| 141 | |
| 142 | |
| 143 | int i, istep; |
| 144 | char c; |
| 145 | char *pos, *startpos, *step; |
| 146 | int n = strlen(str); |
| 147 | |
| 148 | /*enums to define the different lexical stati */ |
| 149 | enum { |
| 150 | sBefore, sNumber, sMinus, sRange, sZero, sSmaller, sError, sSteppedRange |
| 151 | }; |
| 152 | |
| 153 | int status = sBefore; /*status of the deterministic automat to scan str */ |
| 154 | int number = 0; |
| 155 | int end_number = 0; |
| 156 | |
| 157 | char *start = NULL((void*)0); /*holds the string of the number behind a ','*/ |
| 158 | char *end = NULL((void*)0); /*holds the string of the number behind a '-' */ |
| 159 | |
| 160 | int nvecs = 0; /* counts the number of vectors in the list*/ |
| 161 | |
| 162 | step = NULL((void*)0); |
| 163 | snew(pos, n+4)(pos) = save_calloc("pos", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 163, (n+4), sizeof(*(pos))); |
| 164 | startpos = pos; |
| 165 | strcpy(pos, str); |
| 166 | pos[n] = ','; |
| 167 | pos[n+1] = '1'; |
| 168 | pos[n+2] = '\0'; |
| 169 | |
| 170 | *list = NULL((void*)0); |
| 171 | |
| 172 | while ((c = *pos) != 0) |
| 173 | { |
| 174 | switch (status) |
| 175 | { |
| 176 | /* expect a number */ |
| 177 | case sBefore: if (isdigit(c)((*__ctype_b_loc ())[(int) ((c))] & (unsigned short int) _ISdigit )) |
| 178 | { |
| 179 | start = pos; |
| 180 | status = sNumber; |
| 181 | break; |
| 182 | } |
| 183 | else |
| 184 | { |
| 185 | status = sError; |
| 186 | } break; |
| 187 | |
| 188 | /* have read a number, expect ',' or '-' */ |
| 189 | case sNumber: if (c == ',') |
| 190 | { |
| 191 | /*store number*/ |
| 192 | srenew(*list, nvecs+1)(*list) = save_realloc("*list", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 192, (*list), (nvecs+1), sizeof(*(*list))); |
| 193 | (*list)[nvecs++] = number = strtol(start, NULL((void*)0), 10); |
| 194 | status = sBefore; |
| 195 | if (number == 0) |
| 196 | { |
| 197 | status = sZero; |
| 198 | } |
| 199 | break; |
| 200 | } |
| 201 | else if (c == '-') |
| 202 | { |
| 203 | status = sMinus; break; |
| 204 | } |
| 205 | else if (isdigit(c)((*__ctype_b_loc ())[(int) ((c))] & (unsigned short int) _ISdigit )) |
| 206 | { |
| 207 | break; |
| 208 | } |
| 209 | else |
| 210 | { |
| 211 | status = sError; |
| 212 | } break; |
| 213 | |
| 214 | /* have read a '-' -> expect a number */ |
| 215 | case sMinus: |
| 216 | if (isdigit(c)((*__ctype_b_loc ())[(int) ((c))] & (unsigned short int) _ISdigit )) |
| 217 | { |
| 218 | end = pos; |
| 219 | status = sRange; break; |
| 220 | } |
| 221 | else |
| 222 | { |
| 223 | status = sError; |
| 224 | } break; |
| 225 | |
| 226 | case sSteppedRange: |
| 227 | if (isdigit(c)((*__ctype_b_loc ())[(int) ((c))] & (unsigned short int) _ISdigit )) |
| 228 | { |
| 229 | if (step) |
| 230 | { |
| 231 | status = sError; break; |
| 232 | } |
| 233 | else |
| 234 | { |
| 235 | step = pos; |
| 236 | } |
| 237 | status = sRange; |
| 238 | break; |
| 239 | } |
| 240 | else |
| 241 | { |
| 242 | status = sError; |
| 243 | } break; |
| 244 | |
| 245 | /* have read the number after a minus, expect ',' or ':' */ |
| 246 | case sRange: |
| 247 | if (c == ',') |
| 248 | { |
| 249 | /*store numbers*/ |
| 250 | end_number = strtol(end, NULL((void*)0), 10); |
| 251 | number = strtol(start, NULL((void*)0), 10); |
| 252 | status = sBefore; |
| 253 | if (number == 0) |
| 254 | { |
| 255 | status = sZero; break; |
| 256 | } |
| 257 | if (end_number <= number) |
| 258 | { |
| 259 | status = sSmaller; break; |
| 260 | } |
| 261 | srenew(*list, nvecs+end_number-number+1)(*list) = save_realloc("*list", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 261, (*list), (nvecs+end_number-number+1), sizeof(*(*list)) ); |
| 262 | if (step) |
| 263 | { |
| 264 | istep = strtol(step, NULL((void*)0), 10); |
| 265 | step = NULL((void*)0); |
| 266 | } |
| 267 | else |
| 268 | { |
| 269 | istep = 1; |
| 270 | } |
| 271 | for (i = number; i <= end_number; i += istep) |
| 272 | { |
| 273 | (*list)[nvecs++] = i; |
| 274 | } |
| 275 | break; |
| 276 | } |
| 277 | else if (c == ':') |
| 278 | { |
| 279 | status = sSteppedRange; |
| 280 | break; |
| 281 | } |
| 282 | else if (isdigit(c)((*__ctype_b_loc ())[(int) ((c))] & (unsigned short int) _ISdigit )) |
| 283 | { |
| 284 | break; |
| 285 | } |
| 286 | else |
| 287 | { |
| 288 | status = sError; |
| 289 | } break; |
| 290 | |
| 291 | /* format error occured */ |
| 292 | case sError: |
| 293 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 293, "Error in the list of eigenvectors for %s at pos %d with char %c", listname, pos-startpos, *(pos-1)); |
| 294 | break; |
| 295 | /* logical error occured */ |
| 296 | case sZero: |
| 297 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 297, "Error in the list of eigenvectors for %s at pos %d: eigenvector 0 is not valid", listname, pos-startpos); |
| 298 | break; |
| 299 | case sSmaller: |
| 300 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 300, "Error in the list of eigenvectors for %s at pos %d: second index %d is not bigger than %d", listname, pos-startpos, end_number, number); |
| 301 | break; |
| 302 | } |
| 303 | ++pos; /* read next character */ |
| 304 | } /*scanner has finished */ |
| 305 | |
| 306 | /* append zero to list of eigenvectors */ |
| 307 | srenew(*list, nvecs+1)(*list) = save_realloc("*list", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 307, (*list), (nvecs+1), sizeof(*(*list))); |
| 308 | (*list)[nvecs] = 0; |
| 309 | sfree(startpos)save_free("startpos", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 309, (startpos)); |
| 310 | return nvecs; |
| 311 | } /*sscan_list*/ |
| 312 | |
| 313 | void write_eigvec(FILE* fp, int natoms, int eig_list[], rvec** eigvecs, int nvec, const char *grouptitle, real steps[]) |
| 314 | { |
| 315 | /* eig_list is a zero-terminated list of indices into the eigvecs array. |
| 316 | eigvecs are coordinates of eigenvectors |
| 317 | grouptitle to write in the comment line |
| 318 | steps -- array with stepsizes for evLINFIX, evLINACC and evRADACC |
| 319 | */ |
| 320 | |
| 321 | int n = 0, i; rvec x; |
| 322 | real sum; |
| 323 | while (eig_list[n++]) |
| 324 | { |
| 325 | ; /*count selected eigenvecs*/ |
| 326 | |
| 327 | } |
| 328 | fprintf(fp, "# NUMBER OF EIGENVECTORS + %s\n %d\n", grouptitle, n-1); |
| 329 | |
| 330 | /* write list of eigenvector indicess */ |
| 331 | for (n = 0; eig_list[n]; n++) |
| 332 | { |
| 333 | if (steps) |
| 334 | { |
| 335 | fprintf(fp, "%8d %g\n", eig_list[n], steps[n]); |
| 336 | } |
| 337 | else |
| 338 | { |
| 339 | fprintf(fp, "%8d %g\n", eig_list[n], 1.0); |
| 340 | } |
| 341 | } |
| 342 | n = 0; |
| 343 | |
| 344 | /* dump coordinates of the selected eigenvectors */ |
| 345 | while (eig_list[n]) |
| 346 | { |
| 347 | sum = 0; |
| 348 | for (i = 0; i < natoms; i++) |
| 349 | { |
| 350 | if (eig_list[n] > nvec) |
| 351 | { |
| 352 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 352, "Selected eigenvector %d is higher than maximum number %d of available eigenvectors", eig_list[n], nvec); |
| 353 | } |
| 354 | copy_rvec(eigvecs[eig_list[n]-1][i], x); |
| 355 | sum += norm2(x); |
| 356 | fprintf(fp, "%8.5f %8.5f %8.5f\n", x[XX0], x[YY1], x[ZZ2]); |
| 357 | } |
| 358 | n++; |
| 359 | } |
| 360 | } |
| 361 | |
| 362 | |
| 363 | /*enum referring to the different lists of eigenvectors*/ |
| 364 | enum { |
| 365 | evLINFIX, evLINACC, evFLOOD, evRADFIX, evRADACC, evRADCON, evMON, evNr |
| 366 | }; |
| 367 | #define oldMAGIC666 666 |
| 368 | #define MAGIC670 670 |
| 369 | |
| 370 | |
| 371 | void write_the_whole_thing(FILE* fp, t_edipar *edpars, rvec** eigvecs, |
| 372 | int nvec, int *eig_listen[], real* evStepList[]) |
| 373 | { |
| 374 | /* write edi-file */ |
| 375 | |
| 376 | /*Header*/ |
| 377 | fprintf(fp, "#MAGIC\n %d \n#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n", |
| 378 | MAGIC670, edpars->nini, edpars->fitmas, edpars->pcamas); |
| 379 | fprintf(fp, "#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n", |
| 380 | edpars->outfrq, edpars->maxedsteps, edpars->slope); |
| 381 | fprintf(fp, "#PRESTEPS\n %d\n#DELTA_F0\n %f\n#INIT_DELTA_F\n %f\n#TAU\n %f\n#EFL_NULL\n %f\n#ALPHA2\n %f\n#KT\n %f\n#HARMONIC\n %d\n#CONST_FORCE_FLOODING\n %d\n", |
| 382 | edpars->presteps, edpars->flood.deltaF0, edpars->flood.deltaF, edpars->flood.tau, edpars->flood.constEfl, |
| 383 | edpars->flood.alpha2, edpars->flood.kT, edpars->flood.bHarmonic, edpars->flood.bConstForce); |
| 384 | |
| 385 | /* Average and reference positions */ |
| 386 | write_t_edx(fp, edpars->sref, "NREF, XREF"); |
| 387 | write_t_edx(fp, edpars->sav, "NAV, XAV"); |
| 388 | |
| 389 | /*Eigenvectors */ |
| 390 | |
| 391 | write_eigvec(fp, edpars->ned, eig_listen[evMON], eigvecs, nvec, "COMPONENTS GROUP 1", NULL((void*)0)); |
| 392 | write_eigvec(fp, edpars->ned, eig_listen[evLINFIX], eigvecs, nvec, "COMPONENTS GROUP 2", evStepList[evLINFIX]); |
| 393 | write_eigvec(fp, edpars->ned, eig_listen[evLINACC], eigvecs, nvec, "COMPONENTS GROUP 3", evStepList[evLINACC]); |
| 394 | write_eigvec(fp, edpars->ned, eig_listen[evRADFIX], eigvecs, nvec, "COMPONENTS GROUP 4", evStepList[evRADFIX]); |
| 395 | write_eigvec(fp, edpars->ned, eig_listen[evRADACC], eigvecs, nvec, "COMPONENTS GROUP 5", NULL((void*)0)); |
| 396 | write_eigvec(fp, edpars->ned, eig_listen[evRADCON], eigvecs, nvec, "COMPONENTS GROUP 6", NULL((void*)0)); |
| 397 | write_eigvec(fp, edpars->ned, eig_listen[evFLOOD], eigvecs, nvec, "COMPONENTS GROUP 7", evStepList[evFLOOD]); |
| 398 | |
| 399 | |
| 400 | /*Target and Origin positions */ |
| 401 | write_t_edx(fp, edpars->star, "NTARGET, XTARGET"); |
| 402 | write_t_edx(fp, edpars->sori, "NORIGIN, XORIGIN"); |
| 403 | } |
| 404 | |
| 405 | int read_conffile(const char *confin, char *title, rvec *x[]) |
| 406 | { |
| 407 | /* read coordinates out of STX file */ |
| 408 | int natoms; |
| 409 | t_atoms confat; |
| 410 | matrix box; |
| 411 | printf("read coordnumber from file %s\n", confin); |
| 412 | get_stx_coordnum(confin, &natoms); |
| 413 | printf("number of coordinates in file %d\n", natoms); |
| 414 | /* if (natoms != ncoords) |
| 415 | gmx_fatal(FARGS,"number of coordinates in coordinate file (%s, %d)\n" |
| 416 | " does not match topology (= %d)", |
| 417 | confin,natoms,ncoords); |
| 418 | else {*/ |
| 419 | /* make space for coordinates and velocities */ |
| 420 | init_t_atoms(&confat, natoms, FALSE0); |
| 421 | snew(*x, natoms)(*x) = save_calloc("*x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 421, (natoms), sizeof(*(*x))); |
| 422 | read_stx_conf(confin, title, &confat, *x, NULL((void*)0), NULL((void*)0), box); |
| 423 | return natoms; |
| 424 | } |
| 425 | |
| 426 | |
| 427 | void read_eigenvalues(int vecs[], const char *eigfile, real values[], |
| 428 | gmx_bool bHesse, real kT) |
| 429 | { |
| 430 | int neig, nrow, i; |
| 431 | double **eigval; |
| 432 | |
| 433 | neig = read_xvg(eigfile, &eigval, &nrow); |
| 434 | |
| 435 | fprintf(stderrstderr, "Read %d eigenvalues\n", neig); |
| 436 | for (i = bHesse ? 6 : 0; i < neig; i++) |
| 437 | { |
| 438 | if (eigval[1][i] < -0.001 && bHesse) |
| 439 | { |
| 440 | fprintf(stderrstderr, |
| 441 | "WARNING: The Hessian Matrix has negative eigenvalue %f, we set it to zero (no flooding in this direction)\n\n", eigval[1][i]); |
| 442 | } |
| 443 | |
| 444 | if (eigval[1][i] < 0) |
| 445 | { |
| 446 | eigval[1][i] = 0; |
| 447 | } |
| 448 | } |
| 449 | if (bHesse) |
| 450 | { |
| 451 | for (i = 0; vecs[i]; i++) |
| 452 | { |
| 453 | if (vecs[i] < 7) |
| 454 | { |
| 455 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 455, "ERROR: You have chosen one of the first 6 eigenvectors of the HESSE Matrix. That does not make sense, since they correspond to the 6 rotational and translational degrees of freedom.\n\n"); |
| 456 | } |
| 457 | values[i] = eigval[1][vecs[i]-1]/kT; |
| 458 | } |
| 459 | } |
| 460 | else |
| 461 | { |
| 462 | for (i = 0; vecs[i]; i++) |
| 463 | { |
| 464 | if (vecs[i] > (neig-6)) |
| 465 | { |
| 466 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 466, "ERROR: You have chosen one of the last 6 eigenvectors of the COVARIANCE Matrix. That does not make sense, since they correspond to the 6 rotational and translational degrees of freedom.\n\n"); |
| 467 | } |
| 468 | values[i] = 1/eigval[1][vecs[i]-1]; |
| 469 | } |
| 470 | } |
| 471 | /* free memory */ |
| 472 | for (i = 0; i < nrow; i++) |
| 473 | { |
| 474 | sfree(eigval[i])save_free("eigval[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 474, (eigval[i])); |
| 475 | } |
| 476 | sfree(eigval)save_free("eigval", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 476, (eigval)); |
| 477 | } |
| 478 | |
| 479 | |
| 480 | static real *scan_vecparams(const char *str, const char * par, int nvecs) |
| 481 | { |
| 482 | char f0[256], f1[256]; /*format strings adapted every pass of the loop*/ |
| 483 | double d, tcap = 0; |
| 484 | int i; |
| 485 | real *vec_params; |
| 486 | |
| 487 | snew(vec_params, nvecs)(vec_params) = save_calloc("vec_params", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 487, (nvecs), sizeof(*(vec_params))); |
| 488 | if (str) |
| 489 | { |
| 490 | f0[0] = '\0'; |
| 491 | for (i = 0; (i < nvecs); i++) |
| 492 | { |
| 493 | strcpy(f1, f0); /*f0 is the format string for the "to-be-ignored" numbers*/ |
| 494 | strcat(f1, "%lf"); /*and f1 to read the actual number in this pass of the loop*/ |
| 495 | if (sscanf(str, f1, &d) != 1) |
| 496 | { |
| 497 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 497, "Not enough elements for %s parameter (I need %d)", par, nvecs); |
| 498 | } |
| 499 | vec_params[i] = d; |
| 500 | tcap += d; |
| 501 | strcat(f0, "%*s"); |
| 502 | } |
| 503 | } |
| 504 | return vec_params; |
| 505 | } |
| 506 | |
| 507 | |
| 508 | void init_edx(struct edix *edx) |
| 509 | { |
| 510 | edx->nr = 0; |
| 511 | snew(edx->x, 1)(edx->x) = save_calloc("edx->x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 511, (1), sizeof(*(edx->x))); |
| 512 | snew(edx->anrs, 1)(edx->anrs) = save_calloc("edx->anrs", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 512, (1), sizeof(*(edx->anrs))); |
| 513 | } |
| 514 | |
| 515 | void filter2edx(struct edix *edx, int nindex, atom_id index[], int ngro, |
| 516 | atom_id igro[], rvec *x, const char* structure) |
| 517 | { |
| 518 | /* filter2edx copies coordinates from x to edx which are given in index |
| 519 | */ |
| 520 | |
| 521 | int pos, i; |
| 522 | int ix = edx->nr; |
| 523 | edx->nr += nindex; |
| 524 | srenew(edx->x, edx->nr)(edx->x) = save_realloc("edx->x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 524, (edx->x), (edx->nr), sizeof(*(edx->x))); |
| 525 | srenew(edx->anrs, edx->nr)(edx->anrs) = save_realloc("edx->anrs", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 525, (edx->anrs), (edx->nr), sizeof(*(edx->anrs))); |
| 526 | for (i = 0; i < nindex; i++, ix++) |
| 527 | { |
| 528 | for (pos = 0; pos < ngro-1 && igro[pos] != index[i]; ++pos) |
| 529 | { |
| 530 | } |
| 531 | ; /*search element in igro*/ |
| 532 | if (igro[pos] != index[i]) |
| 533 | { |
| 534 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 534, "Couldn't find atom with index %d in structure %s", index[i], structure); |
| 535 | } |
| 536 | edx->anrs[ix] = index[i]; |
| 537 | copy_rvec(x[pos], edx->x[ix]); |
| 538 | } |
| 539 | } |
| 540 | |
| 541 | void get_structure(t_atoms *atoms, const char *IndexFile, |
| 542 | const char *StructureFile, struct edix *edx, int nfit, |
| 543 | atom_id ifit[], int nav, atom_id index[]) |
| 544 | { |
| 545 | atom_id *igro; /*index corresponding to target or origin structure*/ |
| 546 | int ngro; |
| 547 | int ntar; |
| 548 | rvec *xtar; |
| 549 | char title[STRLEN4096]; |
| 550 | char * grpname; |
| 551 | |
| 552 | |
| 553 | ntar = read_conffile(StructureFile, title, &xtar); |
| 554 | printf("Select an index group of %d elements that corresponds to the atoms in the structure file %s\n", |
| 555 | ntar, StructureFile); |
| 556 | get_index(atoms, IndexFile, 1, &ngro, &igro, &grpname); |
| 557 | if (ngro != ntar) |
| 558 | { |
| 559 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 559, "You selected an index group with %d elements instead of %d", ngro, ntar); |
| 560 | } |
| 561 | init_edx(edx); |
| 562 | filter2edx(edx, nfit, ifit, ngro, igro, xtar, StructureFile); |
| 563 | |
| 564 | /* If average and reference/fitting structure differ, append the average structure as well */ |
| 565 | if (ifit != index) /*if fit structure is different append these coordinates, too -- don't mind duplicates*/ |
| 566 | { |
| 567 | filter2edx(edx, nav, index, ngro, igro, xtar, StructureFile); |
| 568 | } |
| 569 | } |
| 570 | |
| 571 | int gmx_make_edi(int argc, char *argv[]) |
| 572 | { |
| 573 | |
| 574 | static const char *desc[] = { |
| 575 | "[THISMODULE] generates an essential dynamics (ED) sampling input file to be used with [TT]mdrun[tt]", |
| 576 | "based on eigenvectors of a covariance matrix ([gmx-covar]) or from a", |
| 577 | "normal modes analysis ([gmx-nmeig]).", |
| 578 | "ED sampling can be used to manipulate the position along collective coordinates", |
| 579 | "(eigenvectors) of (biological) macromolecules during a simulation. Particularly,", |
| 580 | "it may be used to enhance the sampling efficiency of MD simulations by stimulating", |
| 581 | "the system to explore new regions along these collective coordinates. A number", |
| 582 | "of different algorithms are implemented to drive the system along the eigenvectors", |
| 583 | "([TT]-linfix[tt], [TT]-linacc[tt], [TT]-radfix[tt], [TT]-radacc[tt], [TT]-radcon[tt]),", |
| 584 | "to keep the position along a certain (set of) coordinate(s) fixed ([TT]-linfix[tt]),", |
| 585 | "or to only monitor the projections of the positions onto", |
| 586 | "these coordinates ([TT]-mon[tt]).[PAR]", |
| 587 | "References:[BR]", |
| 588 | "A. Amadei, A.B.M. Linssen, B.L. de Groot, D.M.F. van Aalten and ", |
| 589 | "H.J.C. Berendsen; An efficient method for sampling the essential subspace ", |
| 590 | "of proteins., J. Biomol. Struct. Dyn. 13:615-626 (1996)[BR]", |
| 591 | "B.L. de Groot, A. Amadei, D.M.F. van Aalten and H.J.C. Berendsen; ", |
| 592 | "Towards an exhaustive sampling of the configurational spaces of the ", |
| 593 | "two forms of the peptide hormone guanylin,", |
| 594 | "J. Biomol. Struct. Dyn. 13 : 741-751 (1996)[BR]", |
| 595 | "B.L. de Groot, A.Amadei, R.M. Scheek, N.A.J. van Nuland and H.J.C. Berendsen; ", |
| 596 | "An extended sampling of the configurational space of HPr from E. coli", |
| 597 | "Proteins: Struct. Funct. Gen. 26: 314-322 (1996)", |
| 598 | "[PAR]You will be prompted for one or more index groups that correspond to the eigenvectors,", |
| 599 | "reference structure, target positions, etc.[PAR]", |
| 600 | |
| 601 | "[TT]-mon[tt]: monitor projections of the coordinates onto selected eigenvectors.[PAR]", |
| 602 | "[TT]-linfix[tt]: perform fixed-step linear expansion along selected eigenvectors.[PAR]", |
| 603 | "[TT]-linacc[tt]: perform acceptance linear expansion along selected eigenvectors.", |
| 604 | "(steps in the desired directions will be accepted, others will be rejected).[PAR]", |
| 605 | "[TT]-radfix[tt]: perform fixed-step radius expansion along selected eigenvectors.[PAR]", |
| 606 | "[TT]-radacc[tt]: perform acceptance radius expansion along selected eigenvectors.", |
| 607 | "(steps in the desired direction will be accepted, others will be rejected).", |
| 608 | "[BB]Note:[bb] by default the starting MD structure will be taken as origin of the first", |
| 609 | "expansion cycle for radius expansion. If [TT]-ori[tt] is specified, you will be able", |
| 610 | "to read in a structure file that defines an external origin.[PAR]", |
| 611 | "[TT]-radcon[tt]: perform acceptance radius contraction along selected eigenvectors", |
| 612 | "towards a target structure specified with [TT]-tar[tt].[PAR]", |
| 613 | "NOTE: each eigenvector can be selected only once. [PAR]", |
| 614 | "[TT]-outfrq[tt]: frequency (in steps) of writing out projections etc. to [TT].xvg[tt] file[PAR]", |
| 615 | "[TT]-slope[tt]: minimal slope in acceptance radius expansion. A new expansion", |
| 616 | "cycle will be started if the spontaneous increase of the radius (in nm/step)", |
| 617 | "is less than the value specified.[PAR]", |
| 618 | "[TT]-maxedsteps[tt]: maximum number of steps per cycle in radius expansion", |
| 619 | "before a new cycle is started.[PAR]", |
| 620 | "Note on the parallel implementation: since ED sampling is a 'global' thing", |
| 621 | "(collective coordinates etc.), at least on the 'protein' side, ED sampling", |
| 622 | "is not very parallel-friendly from an implementation point of view. Because", |
| 623 | "parallel ED requires some extra communication, expect the performance to be", |
| 624 | "lower as in a free MD simulation, especially on a large number of nodes and/or", |
| 625 | "when the ED group contains a lot of atoms. [PAR]", |
| 626 | "Please also note that if your ED group contains more than a single protein,", |
| 627 | "then the [TT].tpr[tt] file must contain the correct PBC representation of the ED group.", |
| 628 | "Take a look on the initial RMSD from the reference structure, which is printed", |
| 629 | "out at the start of the simulation; if this is much higher than expected, one", |
| 630 | "of the ED molecules might be shifted by a box vector. [PAR]", |
| 631 | "All ED-related output of [TT]mdrun[tt] (specify with [TT]-eo[tt]) is written to a [TT].xvg[tt] file", |
| 632 | "as a function of time in intervals of OUTFRQ steps.[PAR]", |
| 633 | "[BB]Note[bb] that you can impose multiple ED constraints and flooding potentials in", |
| 634 | "a single simulation (on different molecules) if several [TT].edi[tt] files were concatenated", |
| 635 | "first. The constraints are applied in the order they appear in the [TT].edi[tt] file. ", |
| 636 | "Depending on what was specified in the [TT].edi[tt] input file, the output file contains for each ED dataset[PAR]", |
| 637 | "[TT]*[tt] the RMSD of the fitted molecule to the reference structure (for atoms involved in fitting prior to calculating the ED constraints)[BR]", |
| 638 | "[TT]*[tt] projections of the positions onto selected eigenvectors[BR]", |
| 639 | "[PAR][PAR]", |
| 640 | "FLOODING:[PAR]", |
| 641 | "with [TT]-flood[tt], you can specify which eigenvectors are used to compute a flooding potential,", |
| 642 | "which will lead to extra forces expelling the structure out of the region described", |
| 643 | "by the covariance matrix. If you switch -restrain the potential is inverted and the structure", |
| 644 | "is kept in that region.", |
| 645 | "[PAR]", |
| 646 | "The origin is normally the average structure stored in the [TT]eigvec.trr[tt] file.", |
| 647 | "It can be changed with [TT]-ori[tt] to an arbitrary position in configuration space.", |
| 648 | "With [TT]-tau[tt], [TT]-deltaF0[tt], and [TT]-Eflnull[tt] you control the flooding behaviour.", |
| 649 | "Efl is the flooding strength, it is updated according to the rule of adaptive flooding.", |
| 650 | "Tau is the time constant of adaptive flooding, high [GRK]tau[grk] means slow adaption (i.e. growth). ", |
| 651 | "DeltaF0 is the flooding strength you want to reach after tau ps of simulation.", |
| 652 | "To use constant Efl set [TT]-tau[tt] to zero.", |
| 653 | "[PAR]", |
| 654 | "[TT]-alpha[tt] is a fudge parameter to control the width of the flooding potential. A value of 2 has been found", |
| 655 | "to give good results for most standard cases in flooding of proteins.", |
| 656 | "[GRK]alpha[grk] basically accounts for incomplete sampling, if you sampled further the width of the ensemble would", |
| 657 | "increase, this is mimicked by [GRK]alpha[grk] > 1.", |
| 658 | "For restraining, [GRK]alpha[grk] < 1 can give you smaller width in the restraining potential.", |
| 659 | "[PAR]", |
| 660 | "RESTART and FLOODING:", |
| 661 | "If you want to restart a crashed flooding simulation please find the values deltaF and Efl in", |
| 662 | "the output file and manually put them into the [TT].edi[tt] file under DELTA_F0 and EFL_NULL." |
| 663 | }; |
| 664 | |
| 665 | /* Save all the params in this struct and then save it in an edi file. |
| 666 | * ignoring fields nmass,massnrs,mass,tmass,nfit,fitnrs,edo |
| 667 | */ |
| 668 | static t_edipar edi_params; |
| 669 | |
| 670 | enum { |
| 671 | evStepNr = evRADFIX + 1 |
| 672 | }; |
| 673 | static const char* evSelections[evNr] = {NULL((void*)0), NULL((void*)0), NULL((void*)0), NULL((void*)0), NULL((void*)0), NULL((void*)0)}; |
| 674 | static const char* evOptions[evNr] = {"-linfix", "-linacc", "-flood", "-radfix", "-radacc", "-radcon", "-mon"}; |
| 675 | static const char* evParams[evStepNr] = {NULL((void*)0), NULL((void*)0)}; |
| 676 | static const char* evStepOptions[evStepNr] = {"-linstep", "-accdir", "-not_used", "-radstep"}; |
| 677 | static const char* ConstForceStr; |
| 678 | static real * evStepList[evStepNr]; |
| 679 | static real radstep = 0.0; |
| 680 | static real deltaF0 = 150; |
| 681 | static real deltaF = 0; |
| 682 | static real tau = .1; |
| 683 | static real constEfl = 0.0; |
| 684 | static real alpha = 1; |
| 685 | static int eqSteps = 0; |
| 686 | static int * listen[evNr]; |
| 687 | static real T = 300.0; |
| 688 | const real kB = 2.5 / 300.0; /* k_boltzmann in MD units */ |
| 689 | static gmx_bool bRestrain = FALSE0; |
| 690 | static gmx_bool bHesse = FALSE0; |
| 691 | static gmx_bool bHarmonic = FALSE0; |
| 692 | t_pargs pa[] = { |
| 693 | { "-mon", FALSE0, etSTR, {&evSelections[evMON]}, |
| 694 | "Indices of eigenvectors for projections of x (e.g. 1,2-5,9) or 1-100:10 means 1 11 21 31 ... 91" }, |
| 695 | { "-linfix", FALSE0, etSTR, {&evSelections[0]}, |
| 696 | "Indices of eigenvectors for fixed increment linear sampling" }, |
| 697 | { "-linacc", FALSE0, etSTR, {&evSelections[1]}, |
| 698 | "Indices of eigenvectors for acceptance linear sampling" }, |
| 699 | { "-radfix", FALSE0, etSTR, {&evSelections[3]}, |
| 700 | "Indices of eigenvectors for fixed increment radius expansion" }, |
| 701 | { "-radacc", FALSE0, etSTR, {&evSelections[4]}, |
| 702 | "Indices of eigenvectors for acceptance radius expansion" }, |
| 703 | { "-radcon", FALSE0, etSTR, {&evSelections[5]}, |
| 704 | "Indices of eigenvectors for acceptance radius contraction" }, |
| 705 | { "-flood", FALSE0, etSTR, {&evSelections[2]}, |
| 706 | "Indices of eigenvectors for flooding"}, |
| 707 | { "-outfrq", FALSE0, etINT, {&edi_params.outfrq}, |
| 708 | "Freqency (in steps) of writing output in [TT].xvg[tt] file" }, |
| 709 | { "-slope", FALSE0, etREAL, { &edi_params.slope}, |
| 710 | "Minimal slope in acceptance radius expansion"}, |
| 711 | { "-linstep", FALSE0, etSTR, {&evParams[0]}, |
| 712 | "Stepsizes (nm/step) for fixed increment linear sampling (put in quotes! \"1.0 2.3 5.1 -3.1\")"}, |
| 713 | { "-accdir", FALSE0, etSTR, {&evParams[1]}, |
| 714 | "Directions for acceptance linear sampling - only sign counts! (put in quotes! \"-1 +1 -1.1\")"}, |
| 715 | { "-radstep", FALSE0, etREAL, {&radstep}, |
| 716 | "Stepsize (nm/step) for fixed increment radius expansion"}, |
| 717 | { "-maxedsteps", FALSE0, etINT, {&edi_params.maxedsteps}, |
| 718 | "Maximum number of steps per cycle" }, |
| 719 | { "-eqsteps", FALSE0, etINT, {&eqSteps}, |
| 720 | "Number of steps to run without any perturbations "}, |
| 721 | { "-deltaF0", FALSE0, etREAL, {&deltaF0}, |
| 722 | "Target destabilization energy for flooding"}, |
| 723 | { "-deltaF", FALSE0, etREAL, {&deltaF}, |
| 724 | "Start deltaF with this parameter - default 0, nonzero values only needed for restart"}, |
| 725 | { "-tau", FALSE0, etREAL, {&tau}, |
| 726 | "Coupling constant for adaption of flooding strength according to deltaF0, 0 = infinity i.e. constant flooding strength"}, |
| 727 | { "-Eflnull", FALSE0, etREAL, {&constEfl}, |
| 728 | "The starting value of the flooding strength. The flooding strength is updated " |
| 729 | "according to the adaptive flooding scheme. For a constant flooding strength use [TT]-tau[tt] 0. "}, |
| 730 | { "-T", FALSE0, etREAL, {&T}, |
| 731 | "T is temperature, the value is needed if you want to do flooding "}, |
| 732 | { "-alpha", FALSE0, etREAL, {&alpha}, |
| 733 | "Scale width of gaussian flooding potential with alpha^2 "}, |
| 734 | { "-restrain", FALSE0, etBOOL, {&bRestrain}, |
| 735 | "Use the flooding potential with inverted sign -> effects as quasiharmonic restraining potential"}, |
| 736 | { "-hessian", FALSE0, etBOOL, {&bHesse}, |
| 737 | "The eigenvectors and eigenvalues are from a Hessian matrix"}, |
| 738 | { "-harmonic", FALSE0, etBOOL, {&bHarmonic}, |
| 739 | "The eigenvalues are interpreted as spring constant"}, |
| 740 | { "-constF", FALSE0, etSTR, {&ConstForceStr}, |
| 741 | "Constant force flooding: manually set the forces for the eigenvectors selected with -flood " |
| 742 | "(put in quotes! \"1.0 2.3 5.1 -3.1\"). No other flooding parameters are needed when specifying the forces directly."} |
| 743 | }; |
| 744 | #define NPA((int)(sizeof(pa)/sizeof((pa)[0]))) asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))) |
| 745 | |
| 746 | rvec *xref1; |
| 747 | int nvec1, *eignr1 = NULL((void*)0); |
| 748 | rvec *xav1, **eigvec1 = NULL((void*)0); |
| 749 | t_atoms *atoms = NULL((void*)0); |
| 750 | int nav; /* Number of atoms in the average structure */ |
| 751 | char *grpname; |
| 752 | const char *indexfile; |
| 753 | int i; |
| 754 | atom_id *index, *ifit; |
| 755 | int nfit; /* Number of atoms in the reference/fit structure */ |
| 756 | int ev_class; /* parameter _class i.e. evMON, evRADFIX etc. */ |
| 757 | int nvecs; |
| 758 | real *eigval1 = NULL((void*)0); /* in V3.3 this is parameter of read_eigenvectors */ |
| 759 | |
| 760 | const char *EdiFile; |
| 761 | const char *TargetFile; |
| 762 | const char *OriginFile; |
| 763 | const char *EigvecFile; |
| 764 | |
| 765 | output_env_t oenv; |
| 766 | |
| 767 | /*to read topology file*/ |
| 768 | t_topology top; |
| 769 | int ePBC; |
| 770 | char title[STRLEN4096]; |
| 771 | matrix topbox; |
| 772 | rvec *xtop; |
| 773 | gmx_bool bTop, bFit1; |
| 774 | |
| 775 | t_filenm fnm[] = { |
| 776 | { efTRN, "-f", "eigenvec", ffREAD1<<1 }, |
| 777 | { efXVG, "-eig", "eigenval", ffOPTRD(1<<1 | 1<<3) }, |
| 778 | { efTPS, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, |
| 779 | { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
| 780 | { efSTX, "-tar", "target", ffOPTRD(1<<1 | 1<<3)}, |
| 781 | { efSTX, "-ori", "origin", ffOPTRD(1<<1 | 1<<3)}, |
| 782 | { efEDI, "-o", "sam", ffWRITE1<<2 } |
| 783 | }; |
| 784 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) |
| 785 | edi_params.outfrq = 100; edi_params.slope = 0.0; edi_params.maxedsteps = 0; |
| 786 | if (!parse_common_args(&argc, argv, 0, |
| 787 | 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)) |
| 788 | { |
| 789 | return 0; |
| 790 | } |
| 791 | |
| 792 | indexfile = ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
| 793 | EdiFile = ftp2fn(efEDI, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
| 794 | TargetFile = opt2fn_null("-tar", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
| 795 | OriginFile = opt2fn_null("-ori", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
| 796 | |
| 797 | |
| 798 | for (ev_class = 0; ev_class < evNr; ++ev_class) |
| 799 | { |
| 800 | if (opt2parg_bSet(evOptions[ev_class], NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa)) |
| 801 | { |
| 802 | /*get list of eigenvectors*/ |
| 803 | nvecs = sscan_list(&(listen[ev_class]), opt2parg_str(evOptions[ev_class], NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa), evOptions[ev_class]); |
| 804 | if (ev_class < evStepNr-2) |
| 805 | { |
| 806 | /*if apropriate get list of stepsizes for these eigenvectors*/ |
| 807 | if (opt2parg_bSet(evStepOptions[ev_class], NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa)) |
| 808 | { |
| 809 | evStepList[ev_class] = |
| 810 | scan_vecparams(opt2parg_str(evStepOptions[ev_class], NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa), evStepOptions[ev_class], nvecs); |
| 811 | } |
| 812 | else /*if list is not given fill with zeros */ |
| 813 | { |
| 814 | snew(evStepList[ev_class], nvecs)(evStepList[ev_class]) = save_calloc("evStepList[ev_class]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 814, (nvecs), sizeof(*(evStepList[ev_class]))); |
| 815 | for (i = 0; i < nvecs; i++) |
| 816 | { |
| 817 | evStepList[ev_class][i] = 0.0; |
| 818 | } |
| 819 | } |
| 820 | } |
| 821 | else if (ev_class == evRADFIX) |
| 822 | { |
| 823 | snew(evStepList[ev_class], nvecs)(evStepList[ev_class]) = save_calloc("evStepList[ev_class]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 823, (nvecs), sizeof(*(evStepList[ev_class]))); |
| 824 | for (i = 0; i < nvecs; i++) |
| 825 | { |
| 826 | evStepList[ev_class][i] = radstep; |
| 827 | } |
| 828 | } |
| 829 | else if (ev_class == evFLOOD) |
| 830 | { |
| 831 | snew(evStepList[ev_class], nvecs)(evStepList[ev_class]) = save_calloc("evStepList[ev_class]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 831, (nvecs), sizeof(*(evStepList[ev_class]))); |
| 832 | |
| 833 | /* Are we doing constant force flooding? In that case, we read in |
| 834 | * the fproj values from the command line */ |
| 835 | if (opt2parg_bSet("-constF", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa)) |
| 836 | { |
| 837 | evStepList[ev_class] = scan_vecparams(opt2parg_str("-constF", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa), "-constF", nvecs); |
| 838 | } |
| 839 | } |
| 840 | else |
| 841 | { |
| 842 | }; /*to avoid ambiguity */ |
| 843 | } |
| 844 | else /* if there are no eigenvectors for this option set list to zero */ |
| 845 | { |
| 846 | listen[ev_class] = NULL((void*)0); |
| 847 | snew(listen[ev_class], 1)(listen[ev_class]) = save_calloc("listen[ev_class]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 847, (1), sizeof(*(listen[ev_class]))); |
| 848 | listen[ev_class][0] = 0; |
| 849 | } |
| 850 | } |
| 851 | |
| 852 | /* print the interpreted list of eigenvectors - to give some feedback*/ |
| 853 | for (ev_class = 0; ev_class < evNr; ++ev_class) |
| 854 | { |
| 855 | printf("Eigenvector list %7s consists of the indices: ", evOptions[ev_class]); |
| 856 | i = 0; |
| 857 | while (listen[ev_class][i]) |
| 858 | { |
| 859 | printf("%d ", listen[ev_class][i++]); |
| 860 | } |
| 861 | printf("\n"); |
| 862 | } |
| 863 | |
| 864 | EigvecFile = NULL((void*)0); |
| 865 | EigvecFile = opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
| 866 | |
| 867 | /*read eigenvectors from eigvec.trr*/ |
| 868 | read_eigenvectors(EigvecFile, &nav, &bFit1, |
| 869 | &xref1, &edi_params.fitmas, &xav1, &edi_params.pcamas, &nvec1, &eignr1, &eigvec1, &eigval1); |
| 870 | |
| 871 | bTop = read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
Value stored to 'bTop' is never read | |
| 872 | title, &top, &ePBC, &xtop, NULL((void*)0), topbox, 0); |
| 873 | atoms = &top.atoms; |
| 874 | |
| 875 | |
| 876 | printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", nav); |
| 877 | get_index(atoms, indexfile, 1, &i, &index, &grpname); /*if indexfile != NULL parameter 'atoms' is ignored */ |
| 878 | if (i != nav) |
| 879 | { |
| 880 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 880, "you selected a group with %d elements instead of %d", |
| 881 | i, nav); |
| 882 | } |
| 883 | printf("\n"); |
| 884 | |
| 885 | |
| 886 | if (xref1 == NULL((void*)0)) |
| 887 | { |
| 888 | if (bFit1) |
| 889 | { |
| 890 | /* if g_covar used different coordinate groups to fit and to do the PCA */ |
| 891 | printf("\nNote: the structure in %s should be the same\n" |
| 892 | " as the one used for the fit in g_covar\n", ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)); |
| 893 | printf("\nSelect the index group that was used for the least squares fit in g_covar\n"); |
| 894 | } |
| 895 | else |
| 896 | { |
| 897 | printf("\nNote: Apparently no fitting was done in g_covar.\n" |
| 898 | " However, you need to select a reference group for fitting in mdrun\n"); |
| 899 | } |
| 900 | get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname); |
| 901 | snew(xref1, nfit)(xref1) = save_calloc("xref1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_make_edi.c" , 901, (nfit), sizeof(*(xref1))); |
| 902 | for (i = 0; i < nfit; i++) |
| 903 | { |
| 904 | copy_rvec(xtop[ifit[i]], xref1[i]); |
| 905 | } |
| 906 | } |
| 907 | else |
| 908 | { |
| 909 | nfit = nav; |
| 910 | ifit = index; |
| 911 | } |
| 912 | |
| 913 | if (opt2parg_bSet("-constF", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa)) |
| 914 | { |
| 915 | /* Constant force flooding is special: Most of the normal flooding |
| 916 | * options are not needed. */ |
| 917 | edi_params.flood.bConstForce = TRUE1; |
| 918 | } |
| 919 | else |
| 920 | { |
| 921 | /* For normal flooding read eigenvalues and store them in evSteplist[evFLOOD] */ |
| 922 | |
| 923 | if (listen[evFLOOD][0] != 0) |
| 924 | { |
| 925 | read_eigenvalues(listen[evFLOOD], opt2fn("-eig", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), evStepList[evFLOOD], bHesse, kB*T); |
| 926 | } |
| 927 | |
| 928 | edi_params.flood.tau = tau; |
| 929 | edi_params.flood.deltaF0 = deltaF0; |
| 930 | edi_params.flood.deltaF = deltaF; |
| 931 | edi_params.presteps = eqSteps; |
| 932 | edi_params.flood.kT = kB*T; |
| 933 | edi_params.flood.bHarmonic = bHarmonic; |
| 934 | if (bRestrain) |
| 935 | { |
| 936 | /* Trick: invert sign of Efl and alpha2 then this will give the same sign in the exponential and inverted sign outside */ |
| 937 | edi_params.flood.constEfl = -constEfl; |
| 938 | edi_params.flood.alpha2 = -sqr(alpha); |
| 939 | } |
| 940 | else |
| 941 | { |
| 942 | edi_params.flood.constEfl = constEfl; |
| 943 | edi_params.flood.alpha2 = sqr(alpha); |
| 944 | } |
| 945 | } |
| 946 | |
| 947 | edi_params.ned = nav; |
| 948 | |
| 949 | /*number of system atoms */ |
| 950 | edi_params.nini = atoms->nr; |
| 951 | |
| 952 | |
| 953 | /*store reference and average structure in edi_params*/ |
| 954 | make_t_edx(&edi_params.sref, nfit, xref1, ifit ); |
| 955 | make_t_edx(&edi_params.sav, nav, xav1, index); |
| 956 | |
| 957 | |
| 958 | /* Store target positions in edi_params */ |
| 959 | if (opt2bSet("-tar", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
| 960 | { |
| 961 | if (0 != listen[evFLOOD][0]) |
| 962 | { |
| 963 | fprintf(stderrstderr, "\nNote: Providing a TARGET structure has no effect when using flooding.\n" |
| 964 | " You may want to use -ori to define the flooding potential center.\n\n"); |
| 965 | } |
| 966 | get_structure(atoms, indexfile, TargetFile, &edi_params.star, nfit, ifit, nav, index); |
| 967 | } |
| 968 | else |
| 969 | { |
| 970 | make_t_edx(&edi_params.star, 0, NULL((void*)0), index); |
| 971 | } |
| 972 | |
| 973 | /* Store origin positions */ |
| 974 | if (opt2bSet("-ori", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
| 975 | { |
| 976 | get_structure(atoms, indexfile, OriginFile, &edi_params.sori, nfit, ifit, nav, index); |
| 977 | } |
| 978 | else |
| 979 | { |
| 980 | make_t_edx(&edi_params.sori, 0, NULL((void*)0), index); |
| 981 | } |
| 982 | |
| 983 | /* Write edi-file */ |
| 984 | write_the_whole_thing(gmx_ffopen(EdiFile, "w"), &edi_params, eigvec1, nvec1, listen, evStepList); |
| 985 | |
| 986 | return 0; |
| 987 | } |