/* _______________________________________________________________________ * * RDPARM/PTRAJ: APR 2000 * _______________________________________________________________________ * * $Header: /thr/gamow/cvsroot/amber7/src/ptraj/analyze.c,v 1.6 2001/07/13 21:59:19 cheatham Exp $ * * Revision: $Revision: 1.6 $ * Date: $Date: 2001/07/13 21:59:19 $ * Last checked in by $Author: cheatham $ * * This source is now archived under CVS at Scripps in the amber7 tree. * * NOTE: this is a beta, pre-release version, is under constant development, * probably contains some bugs and is under constant revision; therefore * take care. Please report bugs (and fixes!) to cheatham@chpc.utah.edu or * cheatham@cgl.ucsf.edu * * Do not distribute this code without explicit permission. * Do not incorporate into other codes without explicit permission. * * Current contact information: * * Thomas Cheatham, III * 2000 East, 30 South, Skaggs Hall 201 * Department of Medicinal Chemistry * University of Utah * Salt Lake City, UT 84112-5820 * cheatham@chpc.utah.edu * (801) 587-9653 * FAX: (801) 585-5366 * * Other contributors: * * David Case (Scripps) * Michael Crowley (Scripps) * Jed Pitera (UCSF) * Vickie Tsui (Scripps) * * _______________________________________________________________________ * * * The code defined herein is used to analyze and data accumulated by the * various actions during the trajectory processing. */ #include #include #include #include #define ANALYZE_MODULE #include "ptraj.h" int analyzeTest(analyzeInformation *analyze, stackType *sp, int mode) { stackType **argumentStackPointer; if (mode == PTRAJ_SETUP) { /* * ANALYZE: PTRAJ_SETUP * * Parse arguments off the stack and fill in any necessary * information in the analyzeInformation structure. * * This mode is invoked by ptrajSetupAnalyze(). The current * argumentStack is passed in as action->carg1. The * "analyze" structure is already partially setup via * the following code in ptrajSetupAnalyze(): * * statep = ptrajCurrentState(); * state = *statep; * * analyze = (analyzeInformation *) * safe_malloc(sizeof(analyzeInformation)); * INITIALIZE_analyzeInformation(analyze); * * analyze->state = ptrajCopyState(ptrajCurrentState()); * * analyze->type = TRANSFORM_TEST; * analyze->fxn = (analyzeFunction) analyzeTest; */ argumentStackPointer = (stackType **) analyze->carg1; analyze->carg1 = NULL; if (prnlev > 4) { printStack(argumentStackPointer, printString, NULL); } /* * See the comments on argument processing in dispatch.c or in the similar * routine actionTest() in actions.c * * Note that after argument processing, assuming successful return * (i.e. return 0 is performed by this routine), the analyze structure is * placed on the transformAnalyzeStack by ptrajSetupAnalyze() with * the following code: * * pushBottomStack( &transformAnalyzeStack, (void *) action ); * * If an error was encountered during argument processing, return -1 * and this action will be freed and not placed on the * transformActionStack */ printf("In analyzeTest, PTRAJ_SETUP mode\n"); return 0; } if (mode == PTRAJ_STATUS) { /* * Print out a summary of information about this action */ printf("In analyzeTest, PTRAJ_STATUS mode\n"); } /* * Other possible modes: * * PTRAJ_PRINT -- dump information to output files * PTRAJ_CLEANUP -- free any allocated memory * PTRAJ_NOOP -- NULL mode */ if (mode != PTRAJ_ACTION) return 0; /* * Perform action on coordinates, etc. */ printf("In analyzeTest, PTRAJ_ACTION mode\n"); return 1; } int analyzeCorrelationCoefficient(analyzeInformation *analyze, stackType *sp, int mode) { stackType **argumentStackPointer; char *buffer; /* * usage: * * analyze correlationcoefficient {name | ALL} {name | ALL} * * argument usage: * * iarg1: * 0 -- a specific scalarInfo entry for first argument * 1 -- ALL scalarInfo entries for first argument * iarg2: * 0 -- a specific scalarInfo entry for the second argument * 1 -- ALL scalarInfo entries for the second argument * carg1: if iarg1 == 0, the scalarInfo entry * carg2: if iarg2 == 0, the scalarInfo entry */ scalarInfo *scalar1, *scalar2; correlationCoefficientResults Results1, Results2, *results1, *results2; stackType *stack1; int i, total; results1 = &Results1; results2 = &Results2; if (mode == PTRAJ_SETUP) { /* * PTRAJ_SETUP: */ argumentStackPointer = (stackType **) analyze->carg1; analyze->carg1 = NULL; /* * process argument 1 */ buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { fprintf(stdout, "ptraj(), analyzeCorrelationCoefficient: missing the first name on\n"); fprintf(stdout, "scalarStack to process... Returning.\n"); return -1; } if (strcmp(buffer, "ALL") == 0 || strcmp(buffer, "all") == 0) { analyze->iarg1 = 1; } else { scalar1 = scalarStackGetName(&sp, buffer); if (scalar1 == NULL) { /* * try searching for a match without a leading $ just in case the name * was padded */ scalar1 = scalarStackGetName(&sp, buffer+1); } if (scalar1 == NULL) { fprintf(stdout, "ptraj(), analyzeCorrelationCoefficient: cannot find a match in the\n"); fprintf(stdout, "scalarStack for name (%s), returning.\n", buffer); return -1; } analyze->carg1 = (void *) scalar1; } safe_free(buffer); /* * process argument 2 */ buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { fprintf(stdout, "ptraj(), analyzeCorrelationCoefficient: missing the second name on\n"); fprintf(stdout, "scalarStack to process... Returning.\n"); return -1; } if (strcmp(buffer, "ALL") == 0 || strcmp(buffer, "all") == 0) { analyze->iarg2 = 1; } else { scalar2 = scalarStackGetName(&sp, buffer); if (scalar2 == NULL) { /* * try searching for a match without a leading $ just in case the name * was padded */ scalar2 = scalarStackGetName(&sp, buffer+1); } if (scalar2 == NULL) { fprintf(stdout, "ptraj(), analyzeCorrelationCoefficient: cannot find a match in the\n"); fprintf(stdout, "scalarStack for name (%s), returning.\n", buffer); return -1; } analyze->carg2 = (void *) scalar2; } safe_free(buffer); } else if (mode == PTRAJ_STATUS) { /* * PTRAJ_STATUS: */ scalar1 = (scalarInfo *) analyze->carg1; scalar2 = (scalarInfo *) analyze->carg2; fprintf(stdout, " ANALYZE CORRELATIONCOEFFICIENT: comparing %s to %s\n", (analyze->iarg1 ? "ALL" : scalar1->name), (analyze->iarg2 ? "ALL" : scalar2->name)); } /* * Other possible modes: * * PTRAJ_PRINT -- dump information to output files * PTRAJ_CLEANUP -- clean up * PTRAJ_NOOP -- NULL mode */ if (mode != PTRAJ_ACTION) return 0; /* * PTRAJ_ACTION: */ scalar1 = (scalarInfo *) analyze->carg1; scalar2 = (scalarInfo *) analyze->carg2; /* * If ALL is set, count the number of entries on the stack */ total = 0; if (analyze->iarg1 == 1 || analyze->iarg2 == 1) { for (stack1 = sp; stack1 != NULL; stack1 = stack1->next) total++; } /* * Handle the all-to-all case first */ if (analyze->iarg1 == 1 && analyze->iarg2 == 1) { } else if (analyze->iarg1 == 1 || analyze->iarg2 == 1) { if (analyze->iarg2 == 1) scalar1 = scalar2; /* * If ALL is specified (once): no sorting is done */ INITIALIZE_correlationCoefficientResults(results1); INITIALIZE_correlationCoefficientResults(results2); results1->N = scalar1->totalFrames; for (i=0; i < results1->N; i++) { results1->average += scalar1->value[i]; results1->a2 += (scalar1->value[i]*scalar1->value[i]); } results1->average /= results1->N; results1->a2 /= results1->N; results1->stddev = sqrt(results1->a2 - results1->average*results1->average); for (stack1 = sp; stack1 != NULL; stack1 = stack1->next) { scalar2 = (scalarInfo *) stack1->entry; if (scalar1 != scalar2) { results2->N = scalar2->totalFrames; if (results2->N != results1->N) { fprintf(stderr, "Attempting to calculate a correlation coefficient for two data\n"); fprintf(stderr, "sets (%s) and (%s) that have differing numbers of values\n", scalar1->name, scalar2->name); fprintf(stderr, "(%i) and (%i); ignoring this calculation...\n", results1->N, results2->N); } else { for (i=0; i < results2->N; i++) { results2->average += scalar2->value[i]; results2->a2 += (scalar2->value[i]*scalar2->value[i]); } results2->average /= results2->N; results2->a2 /= results2->N; results2->stddev = sqrt(results2->a2 - results2->average*results2->average); for (i=0; i < results1->N; i++) results1->coeff += scalar1->value[i]*scalar2->value[i]; results1->coeff /= results1->N; results1->a = results1->coeff - results1->average - results2->average / (results1->a2 - results2->a2); results1->b = results2->average - results1->a * results1->average; results1->coeff = results1->coeff - (results1->average*results2->average) / (results1->stddev*results2->stddev); if (results1->coeff != 1.0) results1->significance = results1->coeff * sqrt( (results1->N-2)/(1-results1->coeff*results1->coeff)); else results1->significance = 0.0; fprintf(stdout, " CORRELATION COEFFICIENT %6s to %6s IS %10.4f\n", scalar1->name, scalar2->name, results1->coeff); /* fprintf(stdout, " A = %8.4f, B = %8.4f, significance = %8.4f\n", results1->a, results1->b, results1->significance); */ INITIALIZE_correlationCoefficientResults(results2); } } } } else { INITIALIZE_correlationCoefficientResults(results1); INITIALIZE_correlationCoefficientResults(results2); results1->N = scalar1->totalFrames; for (i=0; i < results1->N; i++) { results1->average += scalar1->value[i]; results1->a2 += (scalar1->value[i]*scalar1->value[i]); } results1->average /= results1->N; results1->a2 /= results1->N; results1->stddev = sqrt(results1->a2 - results1->average*results1->average); results2->N = scalar2->totalFrames; if (results2->N != results1->N) { fprintf(stderr, "Attempting to calculate a correlation coefficient for two data\n"); fprintf(stderr, "sets (%s) and (%s) that have differing numbers of values\n", scalar1->name, scalar2->name); fprintf(stderr, "(%i) and (%i); ignoring this calculation...\n", results1->N, results2->N); } else { for (i=0; i < results2->N; i++) { results2->average += scalar2->value[i]; results2->a2 += (scalar2->value[i]*scalar2->value[i]); } results2->average /= results2->N; results2->a2 /= results2->N; results2->stddev = sqrt(results2->a2 - results2->average*results2->average); for (i=0; i < results1->N; i++) results1->coeff += scalar1->value[i]*scalar2->value[i]; results1->coeff /= results1->N; results1->a = (results1->coeff - results1->average - results2->average) / (results1->a2 - results2->a2); results1->b = results2->average - results1->a * results1->average; results1->coeff = (results1->coeff - results1->average*results2->average) / (results1->stddev*results2->stddev); if (results1->coeff != 1.0) results1->significance = results1->coeff * sqrt( (results1->N-2)/(1-results1->coeff*results1->coeff)); else results1->significance = 0.0; fprintf(stdout, " CORRELATION COEFFICIENT %6s to %6s IS %10.4f\n", scalar1->name, scalar2->name, results1->coeff); /* fprintf(stdout, " A = %8.4f, B = %8.4f, significance = %8.4f\n", results1->a, results1->b, results1->significance); */ INITIALIZE_correlationCoefficientResults(results2); } } return 1; } int analyzeHBond(analyzeInformation *analyze, stackType *sp, int mode) { stackType **argumentStackPointer; char *buffer; /* * Usage: * * analyze hbond * [lifetimes [window ]] * [maxoccupied] * [scalar ] * [donor { solvent | resname atomname | mask }] * [acceptor { solvent | res atom res atom | mask } * * argument usage: * * carg1 -- the transformHBondInfo structure pointer * carg2 -- the scalarStack entry * iarg1 -- if set do lifetime analysis, value is window * iarg2 -- if set do maxoccupied analysis, value is window * iarg3 -- if non-zero, a single solvent donor selection * iarg4 -- if non-zero, a single solvent acceptor selection */ transformHBondInfo *hbondInfo, *hbondInfoSelect; scalarInfo *scalar; int *maskD, *maskA, *maskAH1, *maskAH2, *maskAH3; int i, j, k, idonor, iacceptor, iacceptorH, donors, acceptors; int *sortindex, *maxoccupied; float *lifetimes, a, d; int lifetime, lifetimecounter, index, ia, id; if (mode == PTRAJ_SETUP) { argumentStackPointer = (stackType **) analyze->carg1; analyze->carg1 = NULL; /* * name */ buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { fprintf(stdout, "WARNING in ptraj(), analyze: hbond name not specified...\n"); return -1; } hbondInfo = hbondInfoStackGetName(&hbondStack, buffer); if (hbondInfo == NULL) { fprintf(stdout, "WARNING in ptraj(), analyze: hbond name %s not found...\n", buffer); safe_free(buffer); return -1; } safe_free(buffer); analyze->carg1 = (void *) hbondInfo; /* * lifetimes, maxoccupied, window */ analyze->iarg1 = argumentStackContains(argumentStackPointer, "lifetimes"); analyze->iarg2 = argumentStackContains(argumentStackPointer, "maxoccupied"); i = argumentStackKeyToInteger(argumentStackPointer, "window", 0); if (i) { if (analyze->iarg1) analyze->iarg1 = i; if (analyze->iarg2) analyze->iarg2 = i; } /* * scalar * * Place values on the scalar stack using "name" as the key. Note * that this requires selecting an individual h-bond */ scalar = NULL; if (argumentStackContains(argumentStackPointer, "scalar")) { scalar = (scalarInfo *) safe_malloc(sizeof(scalarInfo)); INITIALIZE_scalarInfo(scalar); buffer = argumentStackKeyToString(argumentStackPointer, "distance", NULL); if (buffer != NULL) { scalar->mode = SCALAR_DISTANCE; } else { buffer = argumentStackKeyToString(argumentStackPointer, "angle", NULL); if (buffer == NULL) { fprintf(stdout, "WARNING in ptraj, analyze hbond: scalar command used without\n"); fprintf(stdout, "specifying distance or angle\n"); safe_free(scalar); scalar = NULL; return -1; } scalar->mode = SCALAR_ANGLE; } if ( scalarStackGetName(&scalarStack, buffer) != NULL ) { fprintf(stdout, "WARNING: ptraj(), analyze hbond: The chosen name (%s) has already\n", buffer); fprintf(stdout, "been used. Ignoring scalar command\n"); safe_free(scalar); return -1; } scalar->name = buffer; scalar->totalFrames = -1; } analyze->carg2 = (void *) scalar; /* * SELECTION */ hbondInfoSelect = NULL; maskD = NULL; maskA = NULL; maskAH1 = NULL; maskAH2 = NULL; maskAH3 = NULL; buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) return 0; else if (strcmp(buffer, "donor") != 0) { fprintf(stdout, "WARNING in ptraj, analyze hbond: Error on donor specification in\n"); fprintf(stdout, "select specification... Ignoring command...\n"); return -1; } safe_free(buffer); buffer = (char *) popStack(argumentStackPointer); if ( strstr(buffer, "solvent") != NULL) { if (strcmp(buffer, "solvent") == 0) { analyze->iarg3 = -1; } else if (strncmp(buffer, "solvent", 7) == 0) { sscanf(buffer+7, "%i", &analyze->iarg3); if (analyze->iarg3 > hbondInfo->solventNeighbor) { fprintf(stdout, "WARNING in ptraj, analyze hbond: Number solvent selection out of\n"); fprintf(stdout, "range (%s), max is %i. Ignoring...\n", buffer, hbondInfo->solventNeighbor); return -1; } } safe_free(buffer); } else { pushStack(argumentStackPointer, buffer); maskD = (int *) safe_malloc(sizeof(int) * hbondInfo->state->atoms); for (i=0; i < hbondInfo->state->atoms; i++) maskD[i] = 0; parseHBondDonor(argumentStackPointer, hbondInfo->state, maskD); } buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) return 0; else if (strcmp(buffer, "acceptor") != 0) { fprintf(stdout, "WARNING in ptraj, analyze hbond: Error on acceptor specification in\n"); fprintf(stdout, "select specification... Ignoring command...\n"); return -1; } safe_free(buffer); buffer = (char *) popStack(argumentStackPointer); if ( strstr(buffer, "solvent") != NULL) { if (strcmp(buffer, "solvent") == 0) { analyze->iarg4 = -1; } else if (strncmp(buffer, "solvent", 7) == 0) { sscanf(buffer+7, "%i", &analyze->iarg4); if (analyze->iarg4 > hbondInfo->solventNeighbor) { fprintf(stdout, "WARNING in ptraj, analyze hbond: Number solvent selection out of\n"); fprintf(stdout, "range (%s), max is %i. Ignoring...\n", buffer, hbondInfo->solventNeighbor); return -1; } } safe_free(buffer); } else { pushStack(argumentStackPointer, buffer); maskA = (int *) safe_malloc(sizeof(int) * hbondInfo->state->atoms); maskAH1 = (int *) safe_malloc(sizeof(int) * hbondInfo->state->atoms); maskAH2 = (int *) safe_malloc(sizeof(int) * hbondInfo->state->atoms); maskAH3 = (int *) safe_malloc(sizeof(int) * hbondInfo->state->atoms); for (i=0; i < hbondInfo->state->atoms; i++) { maskA[i] = 0; maskAH1[i] = -1; maskAH2[i] = -1; maskAH3[i] = -1; } parseHBondAcceptor(argumentStackPointer, hbondInfo->state, maskA, maskAH1, maskAH2, maskAH3); } if (maskD || maskA) { hbondInfoSelect = (transformHBondInfo *) safe_malloc(sizeof(transformHBondInfo)); INITIALIZE_transformHBondInfo(hbondInfoSelect); if (maskD) atomMaskIsActive(maskD, hbondInfo->state, &hbondInfoSelect->numdonor, &i); if (maskA) { k = 0; for (i = 0; i < hbondInfo->state->atoms; i++) { if (maskA[i] == 1) { k++; if (maskAH2[i] >= 0) k++; if (maskAH3[i] >= 0) k++; } } hbondInfoSelect->numacceptor = k; } if (maskD) hbondInfoSelect->donor = (int *) safe_malloc(sizeof(int) * hbondInfoSelect->numdonor); if (maskA) { hbondInfoSelect->acceptor = (int *) safe_malloc(sizeof(int) * hbondInfoSelect->numacceptor); hbondInfoSelect->acceptorH = (int *) safe_malloc(sizeof(int) * hbondInfoSelect->numacceptor); } idonor = 0; iacceptor = 0; for (i=0; i < hbondInfo->state->atoms; i++) { if (maskD && maskD[i] == 1) hbondInfoSelect->donor[idonor++] = i; if (maskA && maskA[i] == 1) { hbondInfoSelect->acceptor[iacceptor] = i; hbondInfoSelect->acceptorH[iacceptor] = maskAH1[i]; iacceptor++; if (maskAH2[i] >= 0) { hbondInfoSelect->acceptor[iacceptor] = i; hbondInfoSelect->acceptorH[iacceptor] = maskAH2[i]; iacceptor++; } if (maskAH3[i] >= 0) { hbondInfoSelect->acceptor[iacceptor] = i; hbondInfoSelect->acceptorH[iacceptor] = maskAH3[i]; iacceptor++; } } } } analyze->carg3 = (void *) hbondInfoSelect; if (scalar != NULL) { if (!( (analyze->iarg3 != 0 || (maskD && hbondInfoSelect->numdonor == 1)) && (analyze->iarg4 != 0 || (maskA && hbondInfoSelect->numacceptor == 1)) ) ) { fprintf(stderr, "WARNING in ptraj, analyze hbond: scalar option. More than a single\n"); fprintf(stderr, "hbond selected... This is not supported.\n"); safe_free(scalar); return -1; } } return 0; } hbondInfo = (transformHBondInfo *) analyze->carg1; scalar = (scalarInfo *) analyze->carg2; hbondInfoSelect = (transformHBondInfo *) analyze->carg3; if (mode == PTRAJ_STATUS) { fprintf(stdout, " ANALYZE HBOND: %s -- ", hbondInfo->name); if (analyze->iarg1) fprintf(stdout, "lifetimes, "); if (analyze->iarg2) fprintf(stdout, "max occupied, "); if (analyze->iarg1 > 1 || analyze->iarg2 > 1) fprintf(stdout, "window %i, ", (analyze->iarg1 > 1 ? analyze->iarg1 : analyze->iarg2)); fprintf(stdout, "\n"); if (scalar != NULL) { fprintf(stdout, " Saving hbond %s data to scalar value %s\n", (scalar->mode == SCALAR_DISTANCE ? "distance" : "angle"), scalar->name); } if ( (analyze->iarg3 != 0 || hbondInfoSelect->numdonor > 0) && (analyze->iarg4 != 0 || hbondInfoSelect->numacceptor > 0) ) { fprintf(stdout, " hydrogen bond selection:\n"); if (analyze->iarg3 != 0) { if (analyze->iarg3 < 0) fprintf(stdout, " solvent donor (avg) "); else fprintf(stdout, " solvent donor %i ", analyze->iarg3); } else fprintf(stdout, " donors: %i ", hbondInfoSelect->numdonor); if (analyze->iarg4 != 0) { if (analyze->iarg4 < 0) fprintf(stdout, " solvent acceptor (avg)\n"); else fprintf(stdout, " solvent acceptor %i\n", analyze->iarg4); } else fprintf(stdout, "acceptors %i\n", hbondInfoSelect->numacceptor); if (prnlev > 2 && hbondInfoSelect != NULL) printHBondInfo(hbondInfoSelect->numdonor, hbondInfoSelect->donor, hbondInfoSelect->numacceptor, hbondInfoSelect->acceptor, hbondInfoSelect->acceptorH, hbondInfo->state); } } /* * Other possible modes: * * PTRAJ_PRINT -- dump information to output files * PTRAJ_CLEANUP -- free any allocated memory * PTRAJ_NOOP -- NULL mode */ if (mode != PTRAJ_ACTION) return 0; /* * Perform action on coordinates, etc. */ donors = hbondInfo->numdonor; if (hbondInfo->numSolventAcceptor) donors += hbondInfo->solventNeighbor; acceptors = hbondInfo->numacceptor; if (hbondInfo->numSolventDonor) acceptors += hbondInfo->solventNeighbor; /* * Don't do hbond selection for now... */ if (analyze->iarg1) { lifetimes = (float *) safe_malloc(sizeof(float) * donors * acceptors); for (i=0; i < donors*acceptors; i++) lifetimes[i] = 0.0; } if (analyze->iarg2) { maxoccupied = (int *) safe_malloc(sizeof(int) * donors * acceptors); for (i=0; i < donors*acceptors; i++) maxoccupied[i] = 0; } if (analyze->iarg1) { sortindex = (int *) safe_malloc(sizeof(int) * donors * acceptors); for (k = 0; k < donors*acceptors; k++) sortindex[k] = k; for (i=0; i < hbondInfo->numdonor; i++) { idonor = hbondInfo->donor[i]; for (j = 0; j < hbondInfo->numacceptor; j++) { index = i*acceptors + j; iacceptor = hbondInfo->acceptor[j]; iacceptorH = hbondInfo->acceptorH[j]; lifetime++; lifetimecounter = 0; for (k = 0; k < hbondInfo->visit; k++) { d = hbondInfo->seriesDistance[k*donors*acceptors + index]; a = 180.0 - hbondInfo->seriesAngle[k*donors*acceptors + index]; if (prnlev > 4 ) { printf("series %i, index %i, distance is %.3f, angle is %.3f\n", k, index, d, a); } /* if (hbondInfo->seriesDistance[k*donors*acceptors+index] < hbondInfo->distanceCutoff && (180.0-hbondInfo->seriesAngle[k*donors*acceptors+index])>hbondInfo->angleCutoff) { */ if (d < hbondInfo->distanceCutoff && a > hbondInfo->angleCutoff) { lifetimecounter++; } else { lifetimes[index] += lifetimecounter * hbondInfo->timeinterval; lifetimecounter = 0; lifetime++; } } if (lifetimecounter == hbondInfo->visit) { lifetimes[index] = lifetimecounter * hbondInfo->timeinterval; } else if (lifetime > 0) { lifetimes[index] /= lifetime; } } } sortIndexFloat(lifetimes, sortindex, acceptors*donors); fprintf(stdout, " HBOND LIFETIMES SUMMARY:\n"); fprintf(stdout, " DONOR ACCEPTORH ACCEPTOR\n"); fprintf(stdout, "atom# res# atom res -- atom# res# atom res atom# res# atom res"); fprintf(stdout, " lifetimes\n"); for (k = acceptors*donors-1; k >= 0; k--) { j = sortindex[k]; if (lifetimes[j] > 0.0) { id = j / acceptors; ia = j % acceptors; idonor = hbondInfo->donor[ id ]; iacceptor = hbondInfo->acceptor[ ia ]; iacceptorH = hbondInfo->acceptorH[ ia ]; if ( id <= hbondInfo->numdonor && ia <= hbondInfo->numacceptor && ! (id == hbondInfo->numdonor && ia == hbondInfo->numacceptor) ) { if (id == hbondInfo->numdonor) fprintf(stdout, " solvent donor "); else { fprintf(stdout, "%5i ", idonor+1); printAtom(stdout, idonor, hbondInfo->state); } if (ia == hbondInfo->numacceptor) fprintf(stdout, "-- solvent acceptor "); else { fprintf(stdout, "-- %5i ", iacceptorH+1); printAtom(stdout, iacceptorH, hbondInfo->state); fprintf(stdout, "%5i ", iacceptor+1); printAtom(stdout, iacceptor, hbondInfo->state); } fprintf(stdout," %6.3f\n", lifetimes[j]); } } } safe_free(sortindex); } return 1; } int analyzeSet(analyzeInformation *analyze, stackType *sp, int mode) { stackType **argumentStackPointer; if (mode == PTRAJ_SETUP) { /* * PTRAJ_SETUP: */ argumentStackPointer = (stackType **) analyze->carg1; analyze->carg1 = NULL; if (prnlev > 4) { printStack(argumentStackPointer, printString, NULL); } } else if (mode == PTRAJ_STATUS) { /* * PTRAJ_STATUS: */ } /* * Other possible modes: * * PTRAJ_PRINT -- dump information to output files * PTRAJ_CLEANUP -- free any allocated memory * PTRAJ_NOOP -- NULL mode */ if (mode != PTRAJ_ACTION) return 0; /* * PTRAJ_ACTION: */ return 1; } int analyzeStatistics(analyzeInformation *analyze, stackType *sp, int mode) { /* * usage: * * analyze statistics {name | ALL} [shift value] * * argument usage: * * iarg1: 1 if ALL * darg1: shift value * carg1: the scalarInfo * if iarg1 != 1 */ stackType **argumentStackPointer; char *buffer; scalarInfo *scalar; stackType *stack; double average, stddev, value; int i, periodic; if (mode == PTRAJ_SETUP) { /* * PTRAJ_SETUP: */ argumentStackPointer = (stackType **) analyze->carg1; analyze->carg1 = NULL; if (argumentStackContains(argumentStackPointer, "all")) analyze->iarg1 = 1; else analyze->iarg1 = 0; analyze->darg1 = argumentStackKeyToDouble(argumentStackPointer, "shift", 0.0); if (analyze->iarg1 != 1) { buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { fprintf(stdout, "ptraj(), analyzeStatistics: No name specified\n"); return -1; } scalar = scalarStackGetName(&sp, buffer); if (scalar == NULL) { fprintf(stdout, "ptraj(), analyzeStatistics: Name (%s) not found, ignoring\n", buffer); return -1; } analyze->carg1 = (void *) scalar; } return 0; } scalar = (scalarInfo *) analyze->carg1; if (mode == PTRAJ_STATUS) { /* * PTRAJ_STATUS: */ fprintf(stdout, " ANALYZE STATISTICS: "); if (analyze->iarg1 == 1) fprintf(stdout, "ALL accumulated values "); else fprintf(stdout, "name %s ", scalar->name); if (analyze->darg1 != 0.0) { fprintf(stdout, "shift (about %5.2f) is being applied", analyze->darg1); } fprintf(stdout, "\n"); } /* * Other possible modes: * * PTRAJ_PRINT -- dump information to output files * PTRAJ_CLEANUP -- free any allocated memory * PTRAJ_NOOP -- NULL mode */ if (mode != PTRAJ_ACTION) return 0; /* * PTRAJ_ACTION: */ for (stack = sp; stack != NULL; stack=stack->next) { if (analyze->iarg1 == 1) scalar = (scalarInfo *) stack->entry; average = 0.0; stddev = 0.0; periodic = 0; if (scalar->mode == SCALAR_ANGLE || scalar->mode == SCALAR_TORSION || scalar->mode == SCALAR_PUCKER) periodic = 1; for (i=0; i < scalar->totalFrames; i++) { value = scalar->value[i] - analyze->darg1; if (periodic) { if (value > 180.0) value -= 360.0; else if (value < -180.0) value += 360.0; } average += value; stddev += value*value; } average /= scalar->totalFrames; stddev /= scalar->totalFrames; stddev = stddev - (average*average); if (stddev > 0) stddev = sqrt(stddev); else stddev = 0.0; average += analyze->darg1; fprintf(stdout, " STATISTICS %6s average: %12.4f stddev: %12.4f\n", scalar->name, average, stddev); if (analyze->iarg1 != 1) return 1; } return 1; }