********>Bugfix 29: Author: H. Gohlke, A. Koller Date: 02/05/07 Programs: ptraj Description: Eigenvectors for iRED analysis are output in the wrong order. Fix: apply the following patch to amber9/src/ptraj/actions.c ------------------------------------------------------------------------------ *** src/ptraj/actions.c 2007-01-14 11:13:12.000000000 +0100 --- src/ptraj/actions.c 2007-01-14 12:21:34.000000000 +0100 *************** *** 13619,13625 **** * push the vector info onto the vector stack and store it in action object */ action->carg1 = (void *) vectorInfo; ! pushStack(&vectorStack, (void *) vectorInfo); return 0; } --- 13619,13632 ---- * push the vector info onto the vector stack and store it in action object */ action->carg1 = (void *) vectorInfo; ! if(vectorInfo->mode == VECTOR_IRED){ ! /* Invers vector storage necessary for IRED; ! otherwise IRED matrix will be (N,N)->(1,1) instead of (1,1)->(N,N) */ ! pushBottomStack(&vectorStack, (void *) vectorInfo); ! } ! else{ ! pushStack(&vectorStack, (void *) vectorInfo); ! } return 0; } ------------------------------------------------------------------------------ Temporary workarounds: none ********>Bugfix 29: Author: H. Gohlke, A. Koller Date: 02/05/07 Programs: ptraj Description: - transformAnalyzeMatrix only calculates N-1 eigenvectors - Calculation of Cm(t->T) is missing for iRED analysis - Precision of output for iRED analysis is insufficient Fix: apply the following patch to amber9/src/ptraj/analyze.c ------------------------------------------------------------------------------ *** src/ptraj/analyze.c 2007-01-14 11:13:03.000000000 +0100 --- src/ptraj/analyze.c 2007-02-05 16:23:27.000000000 +0100 *************** *** 452,470 **** * Find eigenvalues and eigenvectors */ ! if(nevec > 0){ ! if(nevec >= nelem){ ! nevec = nelem - 1; ! fprintf(stderr, ! "FYI: NEVEC >= NELEM: Number of calculated evecs were reduced to %i\n", ! nevec); } ! if(nevec >= minfo->snap){ ! nevec = minfo->snap - 1; ! fprintf(stderr, ! "FYI: NEVEC >= SNAP: Number of calculated evecs were reduced to %i\n", ! nevec); } neval = nevec; if(2*nevec <= nelem){ --- 452,496 ---- * Find eigenvalues and eigenvectors */ ! if(nevec > nelem){ ! nevec = nelem; ! fprintf(stderr, ! "FYI: NEVEC > NELEM: Number of calculated evecs were reduced to %i\n", ! nevec); ! } ! if(nevec >= minfo->snap){ ! nevec = minfo->snap; ! fprintf(stderr, ! "FYI: NEVEC > SNAP: Number of calculated evecs were reduced to %i\n", ! nevec); ! } ! ! if(nevec == 0 || nelem == nevec){ ! neval = nelem; ! eigval = (double *) safe_malloc(nelem * sizeof(double)); // for the eigenvalues ! workd = (double *) safe_malloc(3 * nelem * sizeof(double)); // for the function dspev to work with ! uplo[0] = 'L'; // lower triangle of matrix is expected as input for dspev ! ! if (nevec == nelem){ // get all eigenvectors ! vout = (double *) safe_malloc(nelem * nelem * sizeof(double)); // for the eigenvectors ! jobz[0] = 'V'; // calculate eigenvalues and eigenvectors ! ldz = nelem; // dimension of "vout" ! } ! else if(nevec == 0){ // get only eigenvalues; only possible if thermo flag is set, otherwise nevec is set to 1 per default ! vout = (double *) safe_malloc(nelem * sizeof(double)); // for the eigenvectors ! jobz[0] = 'N'; // only calculate eigenvalues ! ldz = 1; // dimension of "vout" } ! ! dspev_(jobz, uplo, &nelem, mat, eigval, vout, &ldz, workd, &info); ! if(info != 0){ ! fprintf(stderr," Warning in analyzeMatrix: dspev returned info = %i\n", info); ! return 0; } + + safe_free(workd); + } + else{ // get up to n-1 eigenvectors neval = nevec; if(2*nevec <= nelem){ *************** *** 525,544 **** safe_free(resid); safe_free(select); } - else{ - neval = nelem; - - eigval = (double *) safe_malloc(nelem * sizeof(double)); - vout = (double *) safe_malloc(nelem * sizeof(double)); - workd = (double *) safe_malloc(3 * nelem * sizeof(double)); - - jobz[0] = 'N'; - uplo[0] = 'L'; - ldz = 1; - dspev_(jobz, uplo, &nelem, mat, eigval, vout, &ldz, workd, &info); - - safe_free(workd); - } if(type == MATRIX_MWCOVAR){ /* --- 551,556 ---- *************** *** 2467,2475 **** if(tcr->table != NULL) safe_free(tcr->table); if(tcr->data1 != NULL) safe_free(tcr->data1); if(tcr->data2 != NULL) safe_free(tcr->data2); ! if(tcr->cf != NULL) safe_free(tcr->cf); ! if(tcr->p2cf != NULL) safe_free(tcr->p2cf); ! if(tcr->rcf != NULL) safe_free(tcr->rcf); INITIALIZE_timecorrResults(tcr); safe_free(tcr); } --- 2479,2488 ---- if(tcr->table != NULL) safe_free(tcr->table); if(tcr->data1 != NULL) safe_free(tcr->data1); if(tcr->data2 != NULL) safe_free(tcr->data2); ! if(tcr->cf != NULL) safe_free(tcr->cf); ! if(tcr->cfinf != NULL) safe_free(tcr->cfinf); ! if(tcr->p2cf != NULL) safe_free(tcr->p2cf); ! if(tcr->rcf != NULL) safe_free(tcr->rcf); INITIALIZE_timecorrResults(tcr); safe_free(tcr); } *************** *** 2497,2503 **** double tstep, tcorr; double rave1, r3iave1, r6iave1, rave2, r3iave2, r6iave2; double dnorm; ! double *table, *data1, *data2, *cf, *p2cf, *rcf; double *avgcrd1, *cftmp1, *p2cftmp1, *rcftmp1, *avgcrd2, *cftmp2, *p2cftmp2, *rcftmp2; double *dpt1, *dpt2; --- 2510,2517 ---- double tstep, tcorr; double rave1, r3iave1, r6iave1, rave2, r3iave2, r6iave2; double dnorm; ! double cfinfavgreal, cfinfavgimg; ! double *table, *data1, *data2, *cf, *cfinf, *p2cf, *rcf; double *avgcrd1, *cftmp1, *p2cftmp1, *rcftmp1, *avgcrd2, *cftmp2, *p2cftmp2, *rcftmp2; double *dpt1, *dpt2; *************** *** 2753,2758 **** --- 2767,2776 ---- timecorr->cf = cf = (double *) safe_malloc(sizeof(double) * nvect * nsteps); for(i = 0; i < nvect * nsteps; i++) cf[i] = 0.0; + + timecorr->cfinf = cfinf = (double *) safe_malloc(sizeof(double) * nvect); + for(i = 0; i < nvect; i++) + cfinf[i] = 0.0; } else if(type == ATCT_NORMAL){ timecorr->p2cf = p2cf = (double *) safe_malloc(sizeof(double) * nsteps); *************** *** 2798,2808 **** /* * Loop over all snapshots */ for(k = 0; k < frame; k++){ ind3 = 2 * (nvect * mtot * k + ind2 + i); ! data1[2*k ] = cftmp1[ind3 ]; ! data1[2*k+1] = cftmp1[ind3+1]; } if(drct){ /* --- 2816,2838 ---- /* * Loop over all snapshots */ + cfinfavgreal = 0.0; + cfinfavgimg = 0.0; for(k = 0; k < frame; k++){ ind3 = 2 * (nvect * mtot * k + ind2 + i); ! data1[2*k ] = cftmp1[ind3 ]; ! cfinfavgreal += cftmp1[ind3 ]; ! data1[2*k+1] = cftmp1[ind3+1]; ! cfinfavgimg += cftmp1[ind3+1]; } + cfinfavgreal /= frame; + cfinfavgimg /= frame; + + /* + * Calc plateau value of correlation function (= C(m,t->T) in Bruschweiler paper (A20)) + */ + + cfinf[i] += (cfinfavgreal * cfinfavgreal) + (cfinfavgimg * cfinfavgimg); if(drct){ /* *************** *** 2950,2977 **** if(type == ATCT_IRED){ ! fprintf(fp, "%10s", "Time"); for(i = 1; i <= nvect; i++){ if(i < 10) ! fprintf(fp, " Mode%i", i); else if(i < 100) ! fprintf(fp, " Mode%i", i); else if(i < 1000) ! fprintf(fp, " Mode%i", i); else if(i < 10000) ! fprintf(fp, " Mode%i", i); else if(i < 100000) ! fprintf(fp, " Mode%i", i); } fprintf(fp, "\n"); for(i = 0; i < nsteps; i++){ ! fprintf(fp, "%10.3f", (double) i * tstep); for(j = 0; j < nvect; j++){ if(norm) ! fprintf(fp, "%10.3f", cf[nsteps * j + i] * frame / (cf[nsteps * j] * (frame - i))); else ! fprintf(fp, "%10.3f", factor * cf[nsteps * j + i] / (frame - i)); } fprintf(fp, "\n"); } --- 2980,3019 ---- if(type == ATCT_IRED){ ! /* Print header */ ! fprintf(fp, "%12s", "XXX"); for(i = 1; i <= nvect; i++){ if(i < 10) ! fprintf(fp, " Mode%i", i); else if(i < 100) ! fprintf(fp, " Mode%i", i); else if(i < 1000) ! fprintf(fp, " Mode%i", i); else if(i < 10000) ! fprintf(fp, " Mode%i", i); else if(i < 100000) ! fprintf(fp, " Mode%i", i); ! } ! fprintf(fp, "\n"); ! ! /* Print cfinf */ ! fprintf(fp, "%12s", "C(m,t->T)"); ! for(i = 0; i < nvect; i++){ ! if(norm) ! fprintf(fp, "%12.8f", cfinf[i] / cf[nsteps * i] * frame); ! else ! fprintf(fp, "%12.8f", cfinf[i]); } fprintf(fp, "\n"); + /* Print cf */ for(i = 0; i < nsteps; i++){ ! fprintf(fp, "%12.8f", (double) i * tstep); for(j = 0; j < nvect; j++){ if(norm) ! fprintf(fp, "%12.8f", cf[nsteps * j + i] * frame / (cf[nsteps * j] * (frame - i))); else ! fprintf(fp, "%12.8f", factor * cf[nsteps * j + i] / (frame - i)); } fprintf(fp, "\n"); } ------------------------------------------------------------------------------ Temporary workarounds: none ********>Bugfix 29: Author: H. Gohlke, A. Koller Date: 02/05/07 Programs: ptraj Description: - Calculation of Cm(t->T) is missing for iRED analysis Fix: apply the following patch to amber9/src/ptraj/analyze.h ------------------------------------------------------------------------------ *** src/ptraj/analyze.h 2007-01-14 14:14:27.000000000 +0100 --- src/ptraj/analyze.h 2007-01-14 13:43:28.000000000 +0100 *************** *** 272,277 **** --- 272,278 ---- double *data1; double *data2; double *cf; + double *cfinf; double *p2cf; double *rcf; } timecorrResults; *************** *** 290,295 **** --- 291,297 ---- _p_->data1 = NULL; \ _p_->data2 = NULL; \ _p_->cf = NULL; \ + _p_->cfinf = NULL; \ _p_->p2cf = NULL; \ _p_->rcf = NULL; ------------------------------------------------------------------------------ Temporary workarounds: none