/* _______________________________________________________________________ * * RDPARM/PTRAJ: APR 2000 * _______________________________________________________________________ * * $Header: /thr/gamow/cvsroot/amber7/src/ptraj/actions.c,v 1.17 2001/07/13 21:59:18 cheatham Exp $ * * Revision: $Revision: 1.17 $ * Date: $Date: 2001/07/13 21:59:18 $ * 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) * _______________________________________________________________________ */ #include #include #include #define ACTION_MODULE #include "ptraj.h" /* * The code in this file implements the various "actions" of ptraj to * perform various coordinate manipulations or analyses. * * The following routines (along with supplemental routines as necessary) * are defined: * * actionTest --- the prototypical action routine (see comments below) * transformAngle --- compute angles and place on the scalarStack * transformAtomicFluct --- compute atomic positional fluctuations or B-factors * transformAverage --- average the coordinates and output to a file * transformCenter --- center the coordinates to the origin or box center * transformCheckOverlap --- check for close contacts * transformClosestWaters --- find the closest set of "n" waters * transformCorr --- perform correlation analysis (Vickie Tsui, Scripps) * transformDihedral --- calculate torsion/dihedral angles * transformDiffusion --- calculate mean squared displacements vs. time * transformDipole --- bin dipoles (Jed Pitera, UCSF) * transformDistance --- compute/store (imaged) distances * transformGrid --- grid atomic densities * transformHBond --- calculate hydrogen bond information * transformImage --- Image molecules outside of a periodic box back in * transformPrincipal --- align coordinates along the principal axis * transformPucker --- compute/store pucker values * transformRadial --- compute radial distribution functions * transformRMS --- perform RMS fitting * transformScale --- scale the coordinates by a specified amount * transformStrip --- strip coordinates * transformTranslate --- shift the coordinates by a specified amount * transformTruncOct --- trim/orient a box to make it a truncated octahedron * transformVector --- compute/store various vector quantities * transformWatershell --- calculate the number of waters in a given shell * * Each of these routines performs its function by being called in a series of * modes defined as follows: * * PTRAJ_SETUP --- perform initialization, * parse arguments off the argumentStack, * complete setup of the actionInformation structure. * NOTE: this works in conjunction with ptrajSetup() * in ptraj.c * * PTRAJ_STATUS --- print useful information about this action, such as * a summary of the arguments. This is called at * startup in ptraj() to print a summary. * * PTRAJ_CLEANUP --- clean up this action, freeing any associated memory. * This mode is applied upon error detection or at the * and of the current round of trajectory processing. * * PTRAJ_ACTION --- perform the actual action on this set of coordinates. * This is called repeatedly for each set of coordinates. * * PTRAJ_PRINT --- Print out any data as necessary; this is currently called * after all of the coordinates have been processed. * * The most complicated mode, and the only mode not directly called by * ptraj() is the PTRAJ_SETUP mode, which involves a detailed obfuscation * from user input file processing to the actual setup of the "transformActionStack" * which contains the list of actions to be performed on each coordinate set. * * This obfuscation is as follows: * * dispatchToken() searches the ptrajTokenlist (both defined in dispatch.c * for a match to a trigger string typed by the user. If this is found, the * remaining text typed on the line is placed into a stack of white space * separated strings called the argumentStack. The subroutine ptrajSetup() * defined in ptraj.c is then called. This routine allocates and initializes * an actionInformation structure, sets this "action" structure to point to the * appropriate function (defined in this file), places the argumentStack into * the complex argument 1 slot of the action structure (action->carg1), and calls * the function in the PTRAJ_SETUP mode. In this mode, it is the functions * responsibility to parse the arguments off of the argument stack and to finish * setup of the actionInformation structure. Upon successful return (i.e. -1 is not * returned), ptrajSetup() places this action onto the stack of actions that will * be processed for every coordinate set, the globally accessible * "transformActionStack". * * The other modes are called during coordinate processing by traversing * though this stack and executing each function. This is performed in * ptraj.c as follows: * * for (actionStackTemp = transformActionStack; * actionStackTemp != NULL; * actionStackTemp = actionStackTemp->next) { * * action = (actionInformation *) actionStackTemp->entry; * if (action->type != TRANSFORM_NOOP && * action->type != TRANSFORM_TRANSFORM) { * * action->fxn(action, X, Y, Z, box, PTRAJ_STATUS); * } * } * * By reading the prototypical "action" routine actionTest below, or any of the * other transformXXX routines defined in this file, this should all become * clearer... * * So, in summary if you would like to add a new "action" that involves * processing of coordinate sets: * * (1) add an entry to the ptrajTokenlist structure in dispatch.c using * the same format as the current entries. Set the expected number * of arguments negative (to create the argumentStack) and set the * function to be called on dispatch, i.e. ptrajSetup() in most cases * (not the function you are defining). You will also define an integer * enumerated type to correspond to this function (i.e. TRANSFORM_ACTION); * this will require modifying the actionType enumerated type defined in * actions.h * * (2) Modify ptrajSetup() in ptraj.c to recognize the new actionType and * to set the action->fxn to point to the new routine you are defining. * * (3) Write the actual routine with the same argument structure as the * other action routines (i.e. action, X, Y, Z, box, mode) and * minimally setup code to handle the PTRAJ_SETUP mode * (for argumentStack processing) and the PTRAJ_ACTION mode for * working on the coordinates. See actionTest() below for detailed * comments. * * NOTES * * o transformDiffusion is currently only setup to image distances in * orthorhombic lattices!!! */ /** ACTION ROUTINE ************************************************************* * * actionTest() --- the prototypical action routine * * This routine demonstrates how an action routine is used and information * can be obtained from the actionInformation structure. It is currently * activated by the keyword "test" in the ptraj command list and therefore * is functional. Information in the actionInformation structure is setup * and initialized in the routine ptrajSetup() in ptraj.c and filled out via * a call to this function in the PTRAJ_SETUP mode. * ******************************************************************************/ int actionTest(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { stackType **argumentStackPointer; ptrajState *state; char *buffer; /* This is a prototype routine for processing the series of * coordinates. Passed into this routine are: * * action -- a structure that contains global control information * (such as the "state") and local information for this * routine. This gets used and setup in this routine. * * x,y,z -- the coordinates * * box -- the box information * * mode -- the current mode (see below). * * The mode controls what is done within the routine and currently this * routine will be called multiple times, each time with a different * mode. The current order is... * * initially, call once at setup, each of: * * PTRAJ_SETUP -- process arguments and setup the action structure * PTRAJ_STATUS -- print a summary of what is going to happen * * call repeatedly for each coordinate frame: * * PTRAJ_ACTION -- do the actual work on each coordinate * * finally, at the end call once, each of: * * PTRAJ_PRINT -- print out any results * PTRAJ_CLEANUP -- clean up memory allocated, etc. * * * The action (transformAction *) structure is a complex beast. * * It contains or will contain all the information necessary for * processing the trajectory files. Initially on the first call * (at PTRAJ_SETUP) it comes in partially setup via the following * code in ptrajSetup(): * * statep = ptrajCurrentState(); * state = *statep; * * action = (actionInformation *) * safe_malloc(sizeof(actionInformation)); * INITIALIZE_actionInformation(action); * * action->state = ptrajCopyState(ptrajCurrentState()); * * action->type = TRANSFORM_TEST; * action->fxn = (actionFunction) actionTest; * * In summary, this: * * -- sets the current state * -- sets the typedef name for this transform function * -- sets a function pointer to this routine. * * Based on this and the processing of command line information, we can * set up the rest of the action structure by setting the values of * various placeholders: * * iarg1 -- integer argument 1 * iarg2 -- integer argument 2 * iarg3 -- integer argument 3 * iarg4 -- integer argument 4 * darg1 -- double argument 1 * darg2 -- double argument 2 * carg1 -- (void *) argument 1 (i.e. a pointer to a structure) * carg2 -- (void *) argument 2 * * Initially, command line arguments come in on action->carg1; these * should be processed at the PTRAJ_SETUP stage as shown below. */ state = (ptrajState *) action->state; if (prnlev > 2) { fprintf(stdout, "In actionTest: currently set to process %i frames\n", state->maxFrames); } if (mode == PTRAJ_SETUP) { /* * ACTION: PTRAJ_SETUP * * Parse arguments off the stack and fill in any necessary * information in the actionInformation "action" structure. * * This mode is invoked by ptrajSetup(). The current * argumentStack is passed in as action->carg1. */ argumentStackPointer = (stackType **) action->carg1; action->carg1 = NULL; /* * To process user input to this "action" command (which is in the form * of text typed in by the user that was placed onto the "argumentStack" * as a series of white space separated strings) it is necessary to * parse the information placed on the "argumentStack" that came into * this routine. There are a variety of routines that can aid in the * processing of these "strings" from the argumentStack (these are * defined in dispatch.c): * * (int ) argumentStackContains(&argumentStack, "foo"); * * Check to see if the string "foo" is a substring of any of the strings * currently on the stack. If it is, remove it from the stack and return * 1 (true), otherwise return 0 (false). * * (char *) argumentStackKeyToString( &argumentStack, "foo", default); * (int) argumentStackKeyToInteger(&argumentStack, "foo", default); * (float) argumentStackKeyToFloat( &argumentStack, "foo", default); * (double) argumentStackKeyToDouble( &argumentStack, "foo", default); * * If the string "foo" is a substring of any of the strings currently * on the stack, remove this entry and grab the next entry converted * to a string, integer, floating point or double precision value, * as appropriate. If "foo" is not found, return the default value. * * (char *) getArgumentString( &argumentStack, default); * (int) getArgumentInteger(&argumentStack, default); * (float) getArgumentFloat( &argumentStack, default); * (double) getArgumentDouble( &argumentStack, default); * * Grab the next string off of the argumentStack and return its * value converted to a string, integer, float or double * depending on the function called. If the string on the stack * is null, return the default value. * * As an example, lets say the possible commands to test are as * follows: * * test mask [alpha | beta | gamma] [sets 100]" * * This can be parsed with the following, setting up the action "mask", * iarg1 and iarg2. */ buffer = getArgumentString(argumentStackPointer, NULL); action->mask = processAtomMask(buffer, action->state); safe_free(buffer); if (argumentStringContains(argumentStackPointer, "alpha")) action->iarg1 = 1; else if (argumentStringContains(argumentStackPointer, "beta")) action->iarg1 = 2; else if (argumentStringContains(argumentStackPointer, "gamma")) action->iarg1 = 3; action->iarg2 = argumentStackKeyToInteger(argumentStackPointer, "set", 0); /* * For more examples, see the various transformSetup routines below. * * * After argument processing, assuming successful return * (i.e. return 0 is performed by this routine), the action is * placed on the transformActionStack by ptrajSetup() with * the following code: * * pushBottomStack( &transformActionStack, (void *) action ); * * If an error was encountered during argument processing, return -1 * and this action will be freed and not placed on the * transformActionStack */ return 0; } if (mode == PTRAJ_STATUS) { /* * Print out a summary of information about this action. * Basically you would output a summary of what specific * options were turned on, etc. */ } /* * 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. */ return 1; } /** ACTION ROUTINE ************************************************************* * * transformAngle --- compute/store angles * ******************************************************************************/ int transformAngle(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { stackType **argumentStackPointer; char *buffer; scalarInfo *info; ptrajState *state; int i; int mask1tot; int mask2tot; int mask3tot; double cx1, cy1, cz1, total_mass1; double cx2, cy2, cz2, total_mass2; double cx3, cy3, cz3, total_mass3; FILE *outFile; /* * USAGE: * * angle name mask1 mask2 mask3 [out ] [time ] * * action argument usage: * * darg1: time interval in ps (for output) * carg1: * a scalarInfo structure */ if (mode == PTRAJ_SETUP) { /* * ACTION: PTRAJ_SETUP */ argumentStackPointer = (stackType **) action->carg1; action->carg1 = NULL; /* * set up the information necessary to place this on the scalarStack */ info = (scalarInfo *) safe_malloc(sizeof(scalarInfo)); INITIALIZE_scalarInfo(info); info->mode = SCALAR_ANGLE; info->totalFrames = -1; info->name = getArgumentString(argumentStackPointer, NULL); if (info->name == NULL) { fprintf(stderr, "WARNING: ptraj(), angle: It is necessary to specify a unique name for\n"); fprintf(stderr, "each angle specified. Ignoring command...\n"); safe_free(info); return -1; } else if ( scalarStackGetName(&scalarStack, info->name) != NULL ) { fprintf(stderr, "WARNING: ptraj(), angle: The chosen name (%s) has already been used.\n", info->name); fprintf(stderr, "Ignoring command...\n"); safe_free(info); return -1; } /* * grab the output filename, if specified */ info->filename = argumentStackKeyToString(argumentStackPointer, "out", NULL); /* * push this entry onto the scalarStack */ pushStack(&scalarStack, (void *) info); /* * grab a time interval between frames in ps (for output) */ action->darg1 = argumentStackKeyToDouble(argumentStackPointer, "time", 1.0); /* * process mask1 --> mask3 */ buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { fprintf(stderr, "WARNING in ptraj(), angle: Error in specification of the first mask\n"); fprintf(stderr, "Ignoring command\n"); safe_free(info); return -1; } else { info->mask1 = processAtomMask(buffer, action->state); safe_free(buffer); } buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { fprintf(stderr, "WARNING in ptraj(), angle: Error in specification of the second mask\n"); fprintf(stderr, "Ignoring command\n"); safe_free(info); return -1; } else { info->mask2 = processAtomMask(buffer, action->state); safe_free(buffer); } buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { fprintf(stderr, "WARNING in ptraj(), angle: Error in specification of the third mask\n"); fprintf(stderr, "Ignoring command\n"); safe_free(info); return -1; } else { info->mask3 = processAtomMask(buffer, action->state); safe_free(buffer); } /* * check to see if each mask only represents a single atom or not * (to save on memory and speed up calculation) */ mask1tot = 0; info->atom1 = -1; mask2tot = 0; info->atom2 = -1; mask3tot = 0; info->atom3 = -1; for (i=0; i < action->state->atoms; i++) { if (info->mask1[i] == 1) { mask1tot++; info->atom1 = i; } if (info->mask2[i] == 1) { mask2tot++; info->atom2 = i; } if (info->mask3[i] == 1) { mask3tot++; info->atom3 = i; } } if (mask1tot == 0) { fprintf(stderr, "WARNING in ptraj(), angle: No atoms selected in mask1, ignoring command\n"); safe_free(info->mask1); safe_free(info); return -1; } else if (mask1tot == 1) { safe_free(info->mask1); info->mask1 = NULL; } else info->atom1 = -1; if (mask2tot == 0) { fprintf(stderr, "WARNING in ptraj(), angle: No atoms selected in mask2, ignoring command\n"); safe_free(info->mask2); safe_free(info); return -1; } else if (mask2tot == 1) { safe_free(info->mask2); info->mask2 = NULL; } else info->atom2 = -1; if (mask3tot == 0) { fprintf(stderr, "WARNING in ptraj(), angle: No atoms selected in mask3, ignoring command\n"); safe_free(info->mask3); safe_free(info); return -1; } else if (mask3tot == 1) { safe_free(info->mask3); info->mask3 = NULL; } else info->atom3 = -1; action->carg1 = (void *) info; return 0; } info = (scalarInfo *) action->carg1; if (mode == PTRAJ_STATUS) { /* * ACTION: PTRAJ_STATUS */ fprintf(stdout, " ANGLE: saved to array named %s\n", info->name); if (info->atom1 == -1) { fprintf(stdout, " Atom selection 1 is "); printAtomMask(info->mask1, action->state); } else { fprintf(stdout, " Atom selection 1 is :%i@%s\n", atomToResidue(info->atom1+1, action->state->residues, action->state->ipres), action->state->atomName[info->atom1]); } if (info->atom2 == -1) { fprintf(stdout, " Atom selection 2 is "); printAtomMask(info->mask2, action->state); } else { fprintf(stdout, " Atom selection 2 is :%i@%s\n", atomToResidue(info->atom2+1, action->state->residues, action->state->ipres), action->state->atomName[info->atom2]); } if (info->atom3 == -1) { fprintf(stdout, " Atom selection 3 is "); printAtomMask(info->mask3, action->state); } else { fprintf(stdout, " Atom selection 3 is :%i@%s\n", atomToResidue(info->atom3+1, action->state->residues, action->state->ipres), action->state->atomName[info->atom3]); } if (info->filename != NULL) { fprintf(stdout, " Data will be dumped to a file named %s\n", info->filename); } } else if (mode == PTRAJ_PRINT) { /* * ACTION: PTRAJ_PRINT */ if ( openFile(&outFile, info->filename, "w") == 0 ) { fprintf(stdout, "WARNING in ptraj(), angle: couldn't open file %s\n", info->filename); return 0; } fprintf(stdout, "PTRAJ ANGLE dumping named values %s\n", info->name); for (i=0; i < action->state->maxFrames; i++) { fprintf(outFile, "%10.2f %f\n", (i+1)*action->darg1, info->value[i]); } safe_fclose(outFile); } else if (mode == PTRAJ_CLEANUP) { /* * ACTION: PTRAJ_CLEANUP */ safe_free(info->name); safe_free(info->filename); safe_free(info->mask1); safe_free(info->mask2); safe_free(info->mask3); safe_free(info->value); INITIALIZE_scalarInfo(info); safe_free(info); } if (mode != PTRAJ_ACTION) return 0; /* * ACTION: PTRAJ_ACTION */ state = (ptrajState *) action->state; /* * update local state information */ for (i=0; i<6; i++) state->box[i] = box[i]; if (info->totalFrames < 0) { info->totalFrames = state->maxFrames; info->value = (double *) safe_malloc(sizeof(double) * info->totalFrames); } if (info->frame > info->totalFrames) { warning("transformAngle()", "Blowing array; too many frames!!\n"); return 0; } cx1 = 0.0; cy1 = 0.0; cz1 = 0.0; total_mass1 = 0.0; cx2 = 0.0; cy2 = 0.0; cz2 = 0.0; total_mass2 = 0.0; cx3 = 0.0; cy3 = 0.0; cz3 = 0.0; total_mass3 = 0.0; if (info->atom1 == -1) { for (i=0; i < state->atoms; i++) { if (info->mask1[i]) { cx1 += state->masses[i] * x[i]; cy1 += state->masses[i] * y[i]; cz1 += state->masses[i] * z[i]; total_mass1 += state->masses[i]; } } cx1 = cx1 / total_mass1; cy1 = cy1 / total_mass1; cz1 = cz1 / total_mass1; } else { cx1 = x[info->atom1]; cy1 = y[info->atom1]; cz1 = z[info->atom1]; } if (info->atom2 == -1) { for (i=0; i < state->atoms; i++) { if (info->mask2[i]) { cx2 += state->masses[i] * x[i]; cy2 += state->masses[i] * y[i]; cz2 += state->masses[i] * z[i]; total_mass2 += state->masses[i]; } } cx2 = cx2 / total_mass2; cy2 = cy2 / total_mass2; cz2 = cz2 / total_mass2; } else { cx2 = x[info->atom2]; cy2 = y[info->atom2]; cz2 = z[info->atom2]; } if (info->atom3 == -1) { for (i=0; i < state->atoms; i++) { if (info->mask3[i]) { cx3 += state->masses[i] * x[i]; cy3 += state->masses[i] * y[i]; cz3 += state->masses[i] * z[i]; total_mass3 += state->masses[i]; } } cx3 = cx3 / total_mass3; cy3 = cy3 / total_mass3; cz3 = cz3 / total_mass3; } else { cx3 = x[info->atom3]; cy3 = y[info->atom3]; cz3 = z[info->atom3]; } info->value[info->frame] = angle(cx1,cy1,cz1,cx2,cy2,cz2,cx3,cy3,cz3); info->frame++; return 1; } /** ACTION ROUTINE ************************************************************* * * transformAtomicFluct() --- compute atomic positional fluctuations * ******************************************************************************/ int transformAtomicFluct(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { stackType **argumentStackPointer; char *buffer; char *filename; coordinateInfo *info; double *xx, *yy, *zz, *results; double xi, yi, zi, fluct, bfactor; int i, j; /* * USAGE: * * atomicfluct [out filename] [] [start ] [stop ] [offset ] * [byres | byatom | bymask] [bfactor] * * action argument usage: * * iarg1: * 0 -- by atom * 1 -- by residue * 2 -- by mask * iarg2: * the number of sets in the average * iarg3: * the number of visits * iarg4: * 0 -- atomic positional fluctuations * 1 -- B factors * carg1: * a coordinate info structure * carg2, carg3, carg4: * the x, y and z coordinates (accumulated) * */ if (mode == PTRAJ_SETUP) { /* * ACTION: PTRAJ_SETUP */ argumentStackPointer = (stackType **) action->carg1; action->carg1 = NULL; filename = argumentStackKeyToString( argumentStackPointer, "out", NULL ); info = (coordinateInfo *) safe_malloc(sizeof(coordinateInfo)); INITIALIZE_coordinateInfo(info); info->file = NULL; info->filename = filename; info->option1 = 0; info->option2 = 0; info->isVelocity = 0; info->info = NULL; info->mask = NULL; info->start = argumentStackKeyToInteger(argumentStackPointer, "start", 1); info->stop = argumentStackKeyToInteger(argumentStackPointer, "stop", -1); if (info->stop == -1) { info->stop = argumentStackKeyToInteger(argumentStackPointer, "end", -1); } info->offset= argumentStackKeyToInteger(argumentStackPointer, "offset", 1); action->iarg1 = 0; if ( argumentStackContains( argumentStackPointer, "byres" ) ) action->iarg1 = 1; else if ( argumentStackContains( argumentStackPointer, "bymask" ) ) action->iarg1 = 2; else if ( argumentStackContains( argumentStackPointer, "byatom" ) ) action->iarg1 = 0; else if ( argumentStackContains( argumentStackPointer, "byatm" ) ) action->iarg1 = 0; action->iarg4 = argumentStackContains( argumentStackPointer, "bfactor" ); buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { action->mask = processAtomMask( (char *) "*", action->state); } else { action->mask = processAtomMask(buffer, action->state); safe_free(buffer); } action->carg1 = (void *) info; xx = (double *) safe_malloc(sizeof(double) * action->state->atoms * 2); yy = (double *) safe_malloc(sizeof(double) * action->state->atoms * 2); zz = (double *) safe_malloc(sizeof(double) * action->state->atoms * 2); for (i=0; i < action->state->atoms*2; i++) { xx[i] = 0.0; yy[i] = 0.0; zz[i] = 0.0; } action->carg2 = (void *) xx; action->carg3 = (void *) yy; action->carg4 = (void *) zz; action->iarg2 = 0; action->iarg3 = 0; return 0; } else if (mode == PTRAJ_STATUS) { /* * ACTION: PTRAJ_STATUS */ info = (coordinateInfo *) action->carg1; fprintf(stdout, " ATOMICFLUCT: dumping %s %s %s", (action->iarg4 ? "B factors" : "atomic positional fluctuations"), (action->iarg1 == 2 ? "by mask" : (action->iarg1 == 1 ? "by residue" : "by atom")), (info->filename == NULL ? "to standard output" : "to file ")); if (info->filename != NULL) fprintf(stdout, "%s\n", info->filename); else fprintf(stdout, "\n"); if (info->start != 1 && info->stop != -1 && info->offset != 1) { fprintf(stdout, " start: %i", info->start); if (info->stop > 0) fprintf(stdout, " stop: %i", info->stop); else fprintf(stdout, " stop [at final frame]"); fprintf(stdout, " offset: %i\n", info->offset); } fprintf(stdout, " Atom selection "); printAtomMask(action->mask, action->state); } else if (mode == PTRAJ_PRINT) { /* * ACTION: PTRAJ_PRINT */ results = (double *) safe_malloc(sizeof(double) * action->state->atoms); for (i=0; i < action->state->atoms; i++) { results[i] = 0.0; } info = (coordinateInfo *) action->carg1; xx = (double *) action->carg2; yy = (double *) action->carg3; zz = (double *) action->carg4; for (i=0; i < action->state->atoms; i++) { xx[i] /= action->iarg2; yy[i] /= action->iarg2; zz[i] /= action->iarg2; } for (i=0; i < action->state->atoms; i++) { xi = xx[i+action->state->atoms]/action->iarg2 - xx[i]*xx[i]; yi = yy[i+action->state->atoms]/action->iarg2 - yy[i]*yy[i]; zi = zz[i+action->state->atoms]/action->iarg2 - zz[i]*zz[i]; xx[i] = xi; yy[i] = yi; zz[i] = zi; } if (info->filename) info->file = safe_fopen(info->filename, "w"); else info->file = stdout; if (info->file == NULL) { fprintf(stdout, "WARNING in ptraj(), atomicfluct: error on opening %s for output\n", (info->filename == NULL ? "(stdout)" : info->filename)); return 0; } fprintf(stdout, "PTRAJ ATOMICFLUCT: Dumping atomic positional fluctuations\n"); if (info->file == stdout) { if (action->iarg1 == 0) fprintf(stdout, " ATOM FLUCTUATION\n"); else if (action->iarg1 == 1) fprintf(stdout, " RES FLUCTUATION\n"); else if (action->iarg1 == 2) fprintf(stdout, " MASK FLUCTUATION\n"); } if (action->iarg4 > 0) bfactor = (8.0/3.0)*PI*PI; else bfactor = 1.0; for (i=0; i < action->state->atoms; i++) { fluct = xx[i] + yy[i] + zz[i]; if (fluct > 0) results[i] = bfactor * sqrt(fluct); else results[i] = 0.0; } if (action->iarg1 == 0) { /* * byatom print out */ j = 0; for (i=0; i < action->state->atoms; i++) { if (action->mask[i]) { if (action->iarg1 == 0) { fprintf(info->file, " %6i %f\n", j+1, results[i]); j++; } } } } else if (action->iarg1 == 1) { /* * byres print out */ for (i=0; i < action->state->residues; i++) { xi = 0.0; fluct = 0.0; for (j=action->state->ipres[i]-1; jstate->ipres[i+1]-1; j++) { if (action->mask[j]) { xi += action->state->masses[j]; fluct += bfactor * results[j] * action->state->masses[j]; } } if (xi > SMALL) fprintf(info->file, " %6i %f\n", i+1, fluct/xi); } } else if (action->iarg1 == 2) { /* * bymask print out */ xi = 0.0; fluct = 0.0; for (i=0; i < action->state->atoms; i++) { if (action->mask[i]) { xi += action->state->masses[i]; fluct += bfactor * results[i] * action->state->masses[i]; } } if (xi > SMALL) fprintf(info->file, " %6i %f\n", 1, fluct/xi); } if (info->file != stdout) { safe_fclose(info->file); info->file = NULL; } } else if (mode == PTRAJ_CLEANUP) { /* * ACTION: PTRAJ_CLEANUP */ info = (coordinateInfo *) action->carg1; xx = (double *) action->carg2; yy = (double *) action->carg3; zz = (double *) action->carg4; safe_free(info); safe_free(xx); safe_free(yy); safe_free(zz); } if (mode != PTRAJ_ACTION) return 0; /* * ACTION: PTRAJ_ACTION */ info = (coordinateInfo *) action->carg1; xx = (double *) action->carg2; yy = (double *) action->carg3; zz = (double *) action->carg4; action->iarg3++; if (action->iarg3 >= info->start && (info->stop < 0 || action->iarg3 <= info->stop) && (action->iarg3 - info->start)%info->offset == 0) { action->iarg2++; for (i=0; i < action->state->atoms; i++) { xx[i] += x[i]; yy[i] += y[i]; zz[i] += z[i]; xx[i+action->state->atoms] += x[i]*x[i]; yy[i+action->state->atoms] += y[i]*y[i]; zz[i+action->state->atoms] += z[i]*z[i]; } } return 1; } /** ACTION ROUTINE ************************************************************* * * transformAverage() --- average (or standard deviate) the coordinates * ******************************************************************************/ int transformAverage(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { stackType **argumentStackPointer; char *buffer; char *filename; coordinateInfo *info, *infotmp; double *xx, *yy, *zz, *tmp_xx, *tmp_yy, *tmp_zz; double fluctx, flucty, fluctz; int i, atoms, total; ptrajState *state; /* * USAGE: * * average filename [] [start ] [stop ] [offset ] * [pdb [parse | dumpq] | binpos | rest] [nobox] [stddev] * * action argument usage: * * iarg1: * 0 -- dump out averages * 1 -- dump out standard deviations * iarg2: * the number of sets in the average * iarg3: * the number of visits * carg1: * the coordinate info structure * carg2, carg3, carg4: * the x, y and z coordinates (accumulated) * */ if (mode == PTRAJ_SETUP) { /* * ACTION: PTRAJ_SETUP */ argumentStackPointer = (stackType **) action->carg1; action->carg1 = NULL; filename = getArgumentString( argumentStackPointer, NULL ); if (filename == NULL) { fprintf(stdout, "WARNING in ptraj(): average command lacks a filename!\n"); safe_free(filename); return -1; } info = (coordinateInfo *) safe_malloc(sizeof(coordinateInfo)); INITIALIZE_coordinateInfo(info); info->filename = filename; info->start = argumentStackKeyToInteger(argumentStackPointer, "start", 1); info->stop = argumentStackKeyToInteger(argumentStackPointer, "stop", -1); info->offset= argumentStackKeyToInteger(argumentStackPointer, "offset", 1); /* * check to see if we do not want box coordinates dumped */ info->isBox = 1; if (argumentStackContains(argumentStackPointer, "nobox" )) info->isBox = 0; /* * check to see if a format other than amber trajectory is wanted, * and note that CHARMM binary is currently not supported... */ if (argumentStackContains( argumentStackPointer, "pdb" )) info->type = COORD_PDB; else if (argumentStackContains( argumentStackPointer, "rest" )) info->type = COORD_AMBER_RESTART; else if (argumentStackContains( argumentStackPointer, "binpos" )) info->type = COORD_BINPOS; else info->type = COORD_AMBER_TRAJECTORY; /* * check to see if we want charges/radii dumped to pdb */ if (argumentStackContains( argumentStackPointer, "dumpq" )) info->option1 = 1; else if (argumentStackContains( argumentStackPointer, "parse" )) info->option1 = 2; if (argumentStackContains( argumentStackPointer, "nowrap" )) info->option2 = 1; if (info->type != COORD_PDB) { info->option1 = 0; info->option2 = 0; } action->iarg1 = argumentStackContains( argumentStackPointer, "stddev" ); buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { action->mask = processAtomMask( (char *) "*", action->state); } else { action->mask = processAtomMask(buffer, action->state); safe_free(buffer); } action->carg1 = (void *) info; if (action->iarg1) total = action->state->atoms * 2; else total = action->state->atoms; xx = (double *) safe_malloc(sizeof(double) * total); yy = (double *) safe_malloc(sizeof(double) * total); zz = (double *) safe_malloc(sizeof(double) * total); for (i=0; i < total; i++) { xx[i] = 0.0; yy[i] = 0.0; zz[i] = 0.0; } action->carg2 = (void *) xx; action->carg3 = (void *) yy; action->carg4 = (void *) zz; action->iarg2 = 0; action->iarg3 = 0; for (i=0; i<6; i++) action->state->box[i] = 0.0; return 0; } else if (mode == PTRAJ_STATUS) { /* * ACTION: PTRAJ_STATUS */ info = (coordinateInfo *) action->carg1; fprintf(stdout, " AVERAGE: dumping the %s of the coordinates to file %s\n", (action->iarg1 ? "standard deviation" : "average"), info->filename); fprintf(stdout, " start: %i", info->start); if (info->stop > 0) fprintf(stdout, " Stop: %i", info->stop); else fprintf(stdout, " Stop [at final frame]"); fprintf(stdout, " Offset: %i\n", info->offset); fprintf(stdout, " Atom selection "); printAtomMask(action->mask, action->state); fprintf(stdout, " Output file information:"); printCoordinateInfo( (void *) info ); } else if (mode == PTRAJ_PRINT) { /* * ACTION: PTRAJ_PRINT */ info = (coordinateInfo *) action->carg1; xx = (double *) action->carg2; yy = (double *) action->carg3; zz = (double *) action->carg4; for (i=0; i < action->state->atoms; i++) { xx[i] /= action->iarg2; yy[i] /= action->iarg2; zz[i] /= action->iarg2; } if (action->iarg1 == 0) for (i=0; i<6; i++) action->state->box[i] /= action->iarg2; if (action->iarg1) { for (i=0; i < action->state->atoms; i++) { xx[i+action->state->atoms] /= action->iarg2; yy[i+action->state->atoms] /= action->iarg2; zz[i+action->state->atoms] /= action->iarg2; fluctx = xx[i+action->state->atoms] - xx[i]*xx[i]; flucty = yy[i+action->state->atoms] - yy[i]*yy[i]; fluctz = zz[i+action->state->atoms] - zz[i]*zz[i]; xx[i] = (fluctx > 0 ? sqrt(fluctx) : 0.0); yy[i] = (flucty > 0 ? sqrt(flucty) : 0.0); zz[i] = (fluctz > 0 ? sqrt(fluctz) : 0.0); } } infotmp = outInfo; outInfo = info; atoms = 0; for (i=0; i < action->state->atoms; i++) if (action->mask[i]) atoms++; if (atoms == action->state->atoms) { ptrajOutputCoordinates(action->state, 1, 0, 1, action->state->atoms, xx, yy, zz, action->state->box); } else { /* * only a subset of the atoms should be dropped; make a copy of the current state, * modify it and also create the corresponding subset of active coordinates. */ tmp_xx = (double *) safe_malloc(sizeof(double) * atoms); tmp_yy = (double *) safe_malloc(sizeof(double) * atoms); tmp_zz = (double *) safe_malloc(sizeof(double) * atoms); atoms = 0; for (i=0; i < action->state->atoms; i++) { if (action->mask[i]) { tmp_xx[atoms] = xx[i]; tmp_yy[atoms] = yy[i]; tmp_zz[atoms] = zz[i]; atoms++; } } modifyStateByMask(&state, &action->state, action->mask, 0); ptrajOutputCoordinates(state, 1, 0, 1, state->atoms, tmp_xx, tmp_yy, tmp_zz, action->state->box); ptrajClearState(&state); safe_free(tmp_xx); safe_free(tmp_yy); safe_free(tmp_zz); } outInfo = infotmp; } else if (mode == PTRAJ_CLEANUP) { /* * ACTION: PTRAJ_CLEANUP */ info = (coordinateInfo *) action->carg1; xx = (double *) action->carg2; yy = (double *) action->carg3; zz = (double *) action->carg4; safe_free(info->filename); INITIALIZE_coordinateInfo(info); safe_free(info); safe_free(x); safe_free(y); safe_free(z); } if (mode != PTRAJ_ACTION) return 0; /* * ACTION: PTRAJ_ACTION */ info = (coordinateInfo *) action->carg1; xx = (double *) action->carg2; yy = (double *) action->carg3; zz = (double *) action->carg4; action->iarg3++; if (action->iarg3 >= info->start && (info->stop < 0 || action->iarg3 <= info->stop) && (action->iarg3 - info->start)%info->offset == 0) { action->iarg2++; if (action->iarg1 == 0) { for (i=0; i<6; i++) action->state->box[i] += box[i]; } for (i=0; i < action->state->atoms; i++) { xx[i] += x[i]; yy[i] += y[i]; zz[i] += z[i]; } if (action->iarg1) { for (i=0; i < action->state->atoms; i++) { xx[i+action->state->atoms] += x[i]*x[i]; yy[i+action->state->atoms] += y[i]*y[i]; zz[i+action->state->atoms] += z[i]*z[i]; } } } return 1; } /** ACTION ROUTINE ************************************************************* * * transformCenter() --- center the coordinates to the origin or box center * ******************************************************************************/ int transformCenter(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { stackType **argumentStackPointer; char *buffer; ptrajState *state; int *mask; int i, origin, com; double cx, cy, cz; double bx, by, bz; double total_mass; /* * USAGE: * * center [origin] [mass] mask * * action argument usage: * * mask: * move the center of mass/geometry of these atoms * iarg1: * 1 -- center to origin (0.0, 0.0, 0.0) * 0 -- center to box center (if there is box information) * iarg2: * 1 -- use center of mass * 0 -- use center of geometry */ if (mode == PTRAJ_SETUP) { /* * ACTION: PTRAJ_SETUP */ argumentStackPointer = (stackType **) action->carg1; action->carg1 = NULL; action->iarg1 = 0; /* center of box */ action->iarg2 = 0; /* center of geometry */ if (argumentStackContains(argumentStackPointer, "origin")) action->iarg1 = 1; if (argumentStackContains(argumentStackPointer, "mass")) action->iarg2 = 1; buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { action->mask = processAtomMask( (char *) "*", action->state); } else { action->mask = processAtomMask(buffer, action->state); safe_free(buffer); } } else if (mode == PTRAJ_STATUS) { /* * ACTION: PTRAJ_STATUS */ fprintf(stdout, " CENTER to %s via center of %s, atom selection follows ", (action->iarg1 ? "origin" : "box center"), (action->iarg2 ? "mass" : "geometry")); printAtomMask(action->mask, action->state); } if (mode != PTRAJ_ACTION) return 0; /* * ACTION: PTRAJ_ACTION */ state = (ptrajState *) action->state; /* * update local state information with current box information */ for (i=0; i<6; i++) state->box[i] = box[i]; /* * process arguments */ mask = action->mask; origin = action->iarg1; com = action->iarg2; if (mask == NULL) return 0; /* * accumulate center of mass or geometry... */ cx = 0.0; cy = 0.0; cz = 0.0; total_mass = 0.0; if ( com ) { for (i=0; i < state->atoms; i++) { if (mask[i]) { cx += state->masses[i] * x[i]; cy += state->masses[i] * y[i]; cz += state->masses[i] * z[i]; total_mass += state->masses[i]; } } } else { for (i=0; i < state->atoms; i++) { if (mask[i]) { cx += x[i]; cy += y[i]; cz += z[i]; total_mass += 1.0; } } } cx /= total_mass; cy /= total_mass; cz /= total_mass; if ( state->IFBOX && origin == 0 ) { bx = state->box[0] / 2.0; by = state->box[1] / 2.0; bz = state->box[2] / 2.0; for (i=0; i < state->atoms; i++) { x[i] = x[i] - cx + bx; y[i] = y[i] - cy + by; z[i] = z[i] - cz + bz; } } else { for (i=0; i < state->atoms; i++) { x[i] -= cx; y[i] -= cy; z[i] -= cz; } } return 1; } /** ACTION ROUTINE ************************************************************* * * transformCheckOverlap --- check for bad atom overlaps * ******************************************************************************/ int transformCheckOverlap(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { stackType **argumentStackPointer; char *buffer; ptrajState *state; int i, j; double distance, ucell[9], recip[9]; /* * USAGE: * * checkoverlap [mask] [min ] [max ] [noimage] * * action argument usage: * * iarg1: 1 implies don't image * darg1: flag distances lower than this * darg2: flag distances greater than this */ if (mode == PTRAJ_SETUP) { /* * ACTION: PTRAJ_SETUP */ argumentStackPointer = (stackType **) action->carg1; action->carg1 = NULL; /* * set up the information necessary to place this on the scalarStack */ /* * grab the min and max */ action->darg1 = argumentStackKeyToDouble(argumentStackPointer, "min", 0.95); action->darg1 = action->darg1 * action->darg1; action->darg2 = argumentStackKeyToDouble(argumentStackPointer, "max", -1.0); if (action->darg2 > 0) { action->darg2 = action->darg2 * action->darg2; } /* * check to see if we want imaging disabled */ action->iarg1 = argumentStackContains(argumentStackPointer, "noimage"); /* * process the atom masks */ buffer = getArgumentString(argumentStackPointer, NULL); if (buffer != NULL) { action->mask = processAtomMask(buffer, action->state); safe_free(buffer); } return 0; } if (mode == PTRAJ_STATUS) { /* * ACTION: PTRAJ_STATUS */ fprintf(stdout, " CHECK OVERLAP: Will flag close contact between atoms if they are\n"); fprintf(stdout, " less than %6.2f angstroms apart\n", sqrt(action->darg1)); if (action->darg2 > 0) { fprintf(stdout, " or if they are greater than %6.2f angstroms apart\n", sqrt(action->darg2)); } if (action->iarg1) fprintf(stdout, " Imaging has been disabled\n"); if (action->mask == NULL) { fprintf(stdout, " Atom selection is: * (all atoms)\n"); } else { fprintf(stdout, " Atom selection is "); printAtomMask(action->mask, action->state); } } if (mode != PTRAJ_ACTION) return 0; /* * ACTION: PTRAJ_ACTION */ state = (ptrajState *) action->state; /* * update local state information */ for (i=0; i<6; i++) state->box[i] = box[i]; if (box[3] <= 0.0 && action->iarg1 == 0) { action->iarg1 = 1; fprintf(stderr, " DISTANCE: box angles are zero, disabling imaging!\n"); } if (action->iarg1 == 0 && (box[3] != 90.0 || box[4] != 90.0 || box[5] != 90.0)) boxToRecip(box, ucell, recip); for (i = 0; i < state->atoms; i++) { for (j = i+1; j < state->atoms; j++) { if (action->mask == NULL || (action->mask[i] && action->mask[j])) { distance = calculateDistance2(i, j, x, y, z, (action->iarg1 ? NULL : box), ucell, recip, 0.0); if (distance < action->darg1) { fprintf(stdout, "OVERLAP: atoms %i (", i+1); printAtomCompact(stdout, i, action->state); fprintf(stdout, ") and %i (", j+1); printAtomCompact(stdout, j, action->state); fprintf(stdout, ") are too close (%8.3f)!\n", sqrt(distance)); } if (action->darg2 > 0 && distance > action->darg2) { fprintf(stdout, "OVERLAP: atoms %i (", i+1); printAtomCompact(stdout, i, action->state); fprintf(stdout, ") and %i (", j+1); printAtomCompact(stdout, j, action->state); fprintf(stdout, ") are too far apart (%8.3f)!\n", sqrt(distance)); } } } } return 1; } /** ACTION ROUTINE ************************************************************* * * transformClosestWaters() --- find the closest set of "n" solvent molecules * * Supplementary routines: * * calculateDistance2 --- calculate a distance**2 and do imaging (ptraj.c) * calculateMinImagedDistance2 --- find shortest distance**2 between images * ******************************************************************************/ int transformClosestWaters(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { stackType **argumentStackPointer, *sp; char *buffer; ptrajState *oldstate, *newstate, **statep; double *minDistance; double *coords; double distance, min; int *sortindex; int ii, i, j; int *mask, closestWaters; int first_only; int *closestWatersMask; double ucell[9], recip[9]; coordinateInfo *refInfo; /* * USAGE: * * closestwaters total [mask] [oxygen | first] * * action argument usage: * * mask: the atom selection around which the closest solvent molecules are saved * iarg1: the number of solvent molecules to save (total) * iarg2: only use the first atom if > 0 (set when oxygen or first specified) */ if (mode == PTRAJ_SETUP) { /* * ACTION: PTRAJ_SETUP */ argumentStackPointer = (stackType **) action->carg1; action->carg1 = NULL; if (action->state->solventMolecules == 0) { fprintf(stdout, "WARNING in ptraj(), closestwaters: This command only works if solvent\n"); fprintf(stdout, "information has been specified. See the \"solvent\" command.\n"); fprintf(stdout, "Ignoring the closestwaters command.\n"); return -1; } action->iarg1 = getArgumentInteger(argumentStackPointer, -1); if (action->iarg1 < 0) { fprintf(stdout, "WARNING in ptraj(), closestwaters: the number of solvent molecules to save is missing\n"); return -1; } else if (action->iarg1 > action->state->solventMolecules) { fprintf(stdout, "WARNING in ptraj(), closestwaters: more solvent molecules are to be saved (%i)\n", action->iarg1); fprintf(stdout, "than are presently listed in the solvent information (%i), %s\n", action->state->solventMolecules, "ignoring command."); return -1; } action->iarg2 = 0; if (argumentStackContains(argumentStackPointer, "oxygen") || argumentStackContains(argumentStackPointer, "first")) { action->iarg2 = 1; } action->iarg3 = 0; if (argumentStackContains(argumentStackPointer, "noimage")) action->iarg3 = 1; /* * check the solvent information to make sure that each solvent listed has the * same number of atoms in each molecule; otherwise a uniform trajectory is not * possible and therefore this command will be ignored... */ j = action->state->solventMoleculeStop[0] - action->state->solventMoleculeStart[0]; for (i=1; i < action->state->solventMolecules; i++) { if (j != (action->state->solventMoleculeStop[i] - action->state->solventMoleculeStart[i])) { fprintf(stdout, "WARNING in ptraj(), closestwaters: the solvent molecules are not of uniform\n"); fprintf(stdout, "size hence this command will be ignored. [Try resetting the solvent\n"); fprintf(stdout, "information with the \"solvent\" command...\n"); return -1; } } buffer = getArgumentString(argumentStackPointer, NULL); action->mask = processAtomMask(buffer, action->state); safe_free(buffer); if (action->mask == NULL) { fprintf(stdout, "WARNING in ptraj(), closestwaters: NULL mask for the solute specified\n"); return -1; } /* * Like with the strip command, we modify the global state pointer at this * point to effectively reduce the size of the trajectory. The oldstate is * kept in complex argument 1. */ oldstate = action->state; closestWatersMask = (int *) safe_malloc(sizeof(int) * oldstate->atoms); /* * set the mask to select all atoms except for the solvent */ for (i=0; i < oldstate->atoms; i++) if (oldstate->solventMask[i]) closestWatersMask[i] = 0; else closestWatersMask[i] = 1; /* * now turn on the total number of solvents that will be saved */ for (i=0; i < action->iarg1; i++) for (j = oldstate->solventMoleculeStart[i]; j < oldstate->solventMoleculeStop[i]; j++) closestWatersMask[j] = 1; /* * modify the state to delete those atoms not selected */ modifyStateByMask(&newstate, &oldstate, closestWatersMask, 0); action->state = newstate; action->carg1 = (void *) oldstate; safe_free(closestWatersMask); /* * update global state variable!!! */ statep = ptrajCurrentState(); *statep = newstate; /* * modify reference structures if set via a recursive call to this routine */ for (sp=transformReferenceStack; sp != NULL; sp=sp->next) { refInfo = (coordinateInfo *) sp->entry; transformClosestWaters(action, refInfo->x, refInfo->y, refInfo->z, action->state->box, PTRAJ_ACTION); } } else if (mode == PTRAJ_STATUS) { /* * ACTION: PTRAJ_STATUS */ oldstate = (ptrajState *) action->carg1; fprintf(stdout, " CLOSESTWATERS: saving the %i closest solvent molecules around atoms ", action->iarg1); printAtomMask(action->mask, oldstate); if (action->iarg3) fprintf(stdout, " Imaging of the coordinates will not be performed\n"); fprintf(stdout, " The current solvent mask is "); printAtomMask(oldstate->solventMask, oldstate); } else if (mode == PTRAJ_CLEANUP) { oldstate = (ptrajState *) action->carg1; ptrajClearState(&oldstate); action->carg1 = NULL; } if (mode != PTRAJ_ACTION) return 0; /* * ACTION: PTRAJ_ACTION */ oldstate = (ptrajState *) action->carg1; newstate = action->state; /* * update local state information */ for (i=0; i<6; i++) newstate->box[i] = box[i]; mask = action->mask; if (mask == NULL) return 0; closestWaters = action->iarg1; first_only = action->iarg2; minDistance = (double *) safe_malloc(sizeof(double) * oldstate->solventMolecules); sortindex = (int *) safe_malloc(sizeof(int) * oldstate->solventMolecules); /* * If non-orthorhomic, find the minimum possible distance between images in * this unit cell; this will save calculation by avoiding the need to * necessarilly calculate distances over all possible images */ if (action->iarg3 == 0 && box[3] == 0.0) { action->iarg3 = 1; fprintf(stderr, " CLOSESTWATERS: box angles are zero, disabling imaging!\n"); } if (action->iarg3 == 0 && (box[3] != 90.0 || box[4] != 90.0 || box[5] != 90.0)) boxToRecip(box, ucell, recip); /* * set the minimum distance to the solute to be larger than * any possible (imaged!!!) distance */ min = box[0] + box[1] + box[2]; min = min * min; for (i=0; i < oldstate->solventMolecules; i++) { minDistance[i] = min; sortindex[i] = i; } /* * loop over all solvent molecules */ for (i=0; i < oldstate->solventMolecules; i++) { /* * loop over all atoms */ for (j=0; j < oldstate->atoms; j++) if (mask[j]) { min = calculateDistance2(oldstate->solventMoleculeStart[i], j, x, y, z, (action->iarg3 ? NULL : box), (double *) ucell, (double *) recip, 0.0); if ( ! first_only ) { for (ii = oldstate->solventMoleculeStart[i]+1; ii < oldstate->solventMoleculeStop[i]; ii++) { distance = calculateDistance2(ii, j, x, y, z, (action->iarg3 ? NULL : box), (double *) ucell, (double *) recip, 0.0); if (distance < min) min = distance; } } if ( minDistance[i] > min ) minDistance[i] = min; } } if (prnlev > 2) { printf("Minimum distances from waters to mask atoms...\n"); for (i = 0; i < oldstate->solventMolecules; i++) { printf("%8.3f ", minDistance[i]); if (i > 0 && (i+1) % 10 == 0) printf("\n"); } printf("\n"); } /* * Sort the minimum distances calculated */ sortIndex(minDistance, sortindex, oldstate->solventMolecules); if (prnlev > 2) { printf("Order of minimum above...\n"); for (i=0; i < oldstate->solventMolecules; i++) { printf("%8i ", sortindex[i]); if (i > 0 && (i+1) % 10 == 0) printf("\n"); } printf("\n"); } coords = (double *) safe_malloc(sizeof(double) * oldstate->atoms); ii = newstate->solventMoleculeStop[0]-newstate->solventMoleculeStart[0]; /* * X */ for (i=0; i < newstate->atoms; i++) coords[i] = x[i]; for (i=0; i < closestWaters; i++) for (j = 0; j < ii; j++) coords[ newstate->solventMoleculeStart[i]+j ] = x[ oldstate->solventMoleculeStart[sortindex[i]]+j ]; for (i=0; i < oldstate->atoms; i++) x[i] = coords[i]; /* * Y */ for (i=0; i < newstate->atoms; i++) coords[i] = y[i]; for (i=0; i < closestWaters; i++) for (j = 0; j < ii; j++) coords[ newstate->solventMoleculeStart[i]+j ] = y[ oldstate->solventMoleculeStart[sortindex[i]]+j ]; for (i=0; i < oldstate->atoms; i++) y[i] = coords[i]; /* * Z */ for (i=0; i < newstate->atoms; i++) coords[i] = z[i]; for (i=0; i < closestWaters; i++) for (j = 0; j < ii; j++) coords[ newstate->solventMoleculeStart[i]+j ] = z[ oldstate->solventMoleculeStart[sortindex[i]]+j ]; for (i=0; i < oldstate->atoms; i++) z[i] = coords[i]; safe_free(coords); safe_free(sortindex); safe_free(minDistance); return 1; } /** ACTION ROUTINE ************************************************************* * * transformCorr() --- perform correlation analysis (Vickie Tsui, Scripps) * * Supplementary routines: * compute_corr() * ******************************************************************************/ void compute_corr(char *outfile, double *x, double *y, double *z, int ts, int tc, int tm, int tf) #define MAXFRAME 5000 { typedef struct _complex { double real; double imaginary; } complex; int i,j, current_frame; double *dipcrd1, *dipcrd2, *dipcrd3; double rmag0, rmagt; double x0, y0, z0, xt, yt, zt; double *p2, *corr, *rcorr; double r6ave, r3ave, rave, avecrd[4], rrig; double dot, y2asy, y20; double th0, phi0; complex y21, y21c, y22, y22c; int *cfind, npts, ncorr, ind0, indt, ntot; int doit, jmax; FILE *ifp; /* allocate space */ jmax = tm/ts + 1; if (tf+1 > jmax) jmax = tf+1; dipcrd1 = (double *) safe_malloc( sizeof(double) * jmax ); dipcrd2 = (double *) safe_malloc( sizeof(double) * jmax ); dipcrd3 = (double *) safe_malloc( sizeof(double) * jmax ); p2 = (double *) safe_malloc( sizeof(double) * jmax ); corr = (double *) safe_malloc( sizeof(double) * jmax ); rcorr = (double *) safe_malloc( sizeof(double) * jmax ); cfind = (int *) safe_malloc( sizeof(int) * jmax ); /* initialize */ for (i=tf; i>=1; --i) { x[i]=x[i-1]; y[i]=y[i-1]; z[i]=z[i-1]; } for (i=1; i<=tf; ++i) { corr[i]=0.0; p2[i]=0.0; rcorr[i]=0.0; cfind[i]=i; } ntot=tm/ts; npts=ntot; ncorr=0; current_frame=ntot; for (i=1; i<=ntot; ++i) { dipcrd1[i]=x[i]; dipcrd2[i]=y[i]; dipcrd3[i]=z[i]; } jmax=ntot; r6ave=0.0; r3ave=0.0; avecrd[1]=0.0; avecrd[2]=0.0; avecrd[3]=0.0; rave=0.0; y2asy=0.0; y20=0.0; y21.real=0.0; y21.imaginary=0.0; y21c.real=0.0; y21c.imaginary=0.0; y22.real=0.0; y22.imaginary=0.0; y22c.real=0.0; y22c.imaginary=0.0; /* main loop for calculating correlation functions */ doit=1; while (doit > 0) { ncorr=ncorr+1; ind0=cfind[1]; rmag0=pow(dipcrd1[ind0],2)+pow(dipcrd2[ind0],2)+pow(dipcrd3[ind0],2); rmag0=sqrt(rmag0); x0=dipcrd1[ind0]/rmag0; y0=dipcrd2[ind0]/rmag0; z0=dipcrd3[ind0]/rmag0; r6ave=r6ave+1/pow(rmag0,6); r3ave=r3ave+1/pow(rmag0,3); rave=rave+rmag0; avecrd[1]=avecrd[1]+dipcrd1[ind0]; avecrd[2]=avecrd[2]+dipcrd2[ind0]; avecrd[3]=avecrd[3]+dipcrd3[ind0]; th0=acos(z0); phi0=atan2(y0,x0); y22.real=y22.real+sqrt(3.0/4.0)*pow((sin(th0)),2)*(cos(2*phi0))/pow(rmag0,3); y22.imaginary=y22.imaginary+sqrt(3.0/4.0)*pow((sin(th0)),2)*(sin(2*phi0))/pow(rmag0,3); y22c.real=y22c.real+sqrt(3.0/4.0)*pow((sin(th0)),2)*(cos(2*phi0))/pow(rmag0,3); y22c.imaginary=y22c.imaginary+sqrt(3.0/4.0)*pow((sin(th0)),2)*(-sin(2*phi0))/pow(rmag0,3); y21.real=y21.real+sqrt(3.0)*cos(th0)*sin(th0)*cos(phi0)/pow(rmag0,3); y21.imaginary=y21.imaginary+sqrt(3.0)*cos(th0)*sin(th0)*sin(phi0)/pow(rmag0,3); y21c.real=y21c.real+sqrt(3.0)*cos(th0)*sin(th0)*cos(phi0)/pow(rmag0,3); y21c.imaginary=y21c.imaginary+sqrt(3.0)*cos(th0)*sin(th0)*(-sin(phi0))/pow(rmag0,3); y20=y20+(pow((3*(cos(th0))),2)-1)/(2.0*pow(rmag0,3)); for (j=1; j<=jmax; ++j) { indt=cfind[j]; rmagt=pow(dipcrd1[indt],2)+pow(dipcrd2[indt],2)+pow(dipcrd3[indt],2); rmagt=sqrt(rmagt); xt=dipcrd1[indt]/rmagt; yt=dipcrd2[indt]/rmagt; zt=dipcrd3[indt]/rmagt; dot=(3*pow((x0*xt+y0*yt+z0*zt),2)-1)/2.0; corr[j]=corr[j]+dot/pow((rmag0*rmagt),3); p2[j]=p2[j]+dot; rcorr[j]=rcorr[j]+1/pow((rmag0*rmagt),3); } if (ncorr != npts) { for (j=1; j<=jmax-1; ++j) { cfind[j]=cfind[j+1]; } cfind[jmax]=ind0; current_frame=current_frame+1; if (current_frame < tf) { dipcrd1[current_frame]=x[current_frame]; dipcrd2[current_frame]=y[current_frame]; dipcrd3[current_frame]=z[current_frame]; npts=npts+1; } else { jmax=jmax-1; } } else { doit=0; } } /* normalize correlation functions */ r6ave=r6ave/npts; r3ave=r3ave/npts; rave=rave/npts; avecrd[1]=avecrd[1]/npts; avecrd[2]=avecrd[2]/npts; avecrd[3]=avecrd[3]/npts; rrig=pow(avecrd[1],2)+pow(avecrd[2],2)+pow(avecrd[3],2); rrig=sqrt(rrig); y2asy=(y22.real*y22c.real+y21.real*y21c.real)+pow(y20,2); y2asy=y2asy/(npts*npts*r6ave); for (i=1; i<=ntot; ++i) { corr[i]=corr[i]/((npts-i+1)*r6ave); rcorr[i]=rcorr[i]/((npts-i+1)*r6ave); p2[i]=p2[i]/(npts-i+1); } /* output correlation functions */ ifp=safe_fopen(outfile, "w"); if (ifp == NULL) { warning("ptraj(), correlation: cannot open output file %s\n", outfile); } else { fprintf(ifp, "# Rrigid= %lf Rave= %lf \n", rrig, rave); fprintf(ifp, "# <1/r^3>= %lf <1/r^6>= %lf\n", r3ave, r6ave); /* rfac = r6ave*pow(rave,6); qfac = y2asy*rfac; */ fprintf(ifp, "# time C(t) P2(t) R(t)\n"); i=tc/ts; for (j=1; j<=i; ++j) { fprintf(ifp, "%d %lf %lf %lf\n", (j-1)*ts, corr[j], p2[j], rcorr[j]); } safe_fclose(ifp); } /* deallocate space */ safe_free(dipcrd1); safe_free(dipcrd2); safe_free(dipcrd3); safe_free(p2); safe_free(corr); safe_free(rcorr); safe_free(cfind); } int transformCorr(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { stackType **argumentStackPointer; char *buffer; transformCorrInfo *corrInfo; ptrajState *state; int i; double cx, cy, cz, total_mass; double vx, vy, vz; /* * USAGE: * * correlation name mask1 mask2 tmin tcorr tmax [out filename] * * action argument usage: * * carg1: * a transformCorrInfo structure */ if (mode == PTRAJ_SETUP) { /* * ACTION: PTRAJ_SETUP */ argumentStackPointer = (stackType **) action->carg1; action->carg1 = NULL; /* * set up complex argument */ corrInfo = (transformCorrInfo *) safe_malloc(sizeof(transformCorrInfo)); INITIALIZE_transformCorrInfo(corrInfo); corrInfo->totalFrames = -1; corrInfo->name = getArgumentString(argumentStackPointer, NULL); buffer = getArgumentString(argumentStackPointer, NULL); corrInfo->mask = processAtomMask(buffer, action->state); safe_free(buffer); buffer = getArgumentString(argumentStackPointer, NULL); corrInfo->mask2 = processAtomMask(buffer, action->state); safe_free(buffer); corrInfo->mode = VECTOR_MASK; corrInfo->tmin = getArgumentInteger(argumentStackPointer, 1.0); corrInfo->tcorr = getArgumentInteger(argumentStackPointer, 1.0); corrInfo->tmax = getArgumentInteger(argumentStackPointer, 1.0); /* * assume "out" may be missing */ corrInfo->filename = argumentStackKeyToString(argumentStackPointer, "out", NULL); if (corrInfo->filename == NULL) { corrInfo->filename = getArgumentString(argumentStackPointer, NULL); if (corrInfo->filename == NULL) { error("ptraj()", "correlation, no out file specified\n"); } } if (corrInfo->name == NULL || corrInfo->mask == NULL || corrInfo->mask2 == NULL) { error("ptraj()", "correlation arguments\n"); } action->carg1 = (void *) corrInfo; return 0; } corrInfo = (transformCorrInfo *) action->carg1; if (mode == PTRAJ_STATUS) { /* * ACTION: PTRAJ_STATUS */ fprintf(stdout, " CORRELATION: storage to array named %s", corrInfo->name); fprintf(stdout, " -- tmin: %i tcorr: %i tmax: %i\n", corrInfo->tmin, corrInfo->tcorr, corrInfo->tmax); fprintf(stdout, " Atom selection 1 is "); printAtomMask(corrInfo->mask, action->state); fprintf(stdout, " Atom selection 2 is "); printAtomMask(corrInfo->mask2, action->state); } else if (mode == PTRAJ_PRINT) { /* * ACTION: PTRAJ_PRINT */ fprintf(stdout, "PTRAJ CORRELATION: calculating correlation %s\n", corrInfo->name); if (corrInfo != NULL) { compute_corr(corrInfo->filename, corrInfo->vx, corrInfo->vy, corrInfo->vz, corrInfo->tmin, corrInfo->tcorr, corrInfo->tmax, corrInfo->totalFrames); } return 0; } else if (mode == PTRAJ_CLEANUP) { /* * ACTION: PTRAJ_CLEANUP */ safe_free(corrInfo->cx); safe_free(corrInfo->cy); safe_free(corrInfo->cz); safe_free(corrInfo->vx); safe_free(corrInfo->vy); safe_free(corrInfo->vz); safe_free(corrInfo->mask); safe_free(corrInfo->mask2); safe_free(corrInfo->name); INITIALIZE_transformCorrInfo(corrInfo); safe_free(corrInfo); } if (mode != PTRAJ_ACTION) return 0; /* * ACTION: PTRAJ_ACTION */ state = (ptrajState *) action->state; if (corrInfo->totalFrames < 0) { corrInfo->totalFrames = state->maxFrames; corrInfo->cx = (double *) safe_malloc(sizeof(double) * corrInfo->totalFrames); corrInfo->cy = (double *) safe_malloc(sizeof(double) * corrInfo->totalFrames); corrInfo->cz = (double *) safe_malloc(sizeof(double) * corrInfo->totalFrames); corrInfo->vx = (double *) safe_malloc(sizeof(double) * (corrInfo->totalFrames+1)); corrInfo->vy = (double *) safe_malloc(sizeof(double) * (corrInfo->totalFrames+1)); corrInfo->vz = (double *) safe_malloc(sizeof(double) * (corrInfo->totalFrames+1)); } if (corrInfo->frame > corrInfo->totalFrames) { warning("transformCorrr()", "Blowing array; too many frames!!\n"); return 0; } corrInfo->mode = CORR_MASK; total_mass = 0.0; cx = 0.0; cy = 0.0; cz = 0.0; for (i=0; i < state->atoms; i++) { if (corrInfo->mask[i]) { cx += state->masses[i] * x[i]; cy += state->masses[i] * y[i]; cz += state->masses[i] * z[i]; total_mass += state->masses[i]; } } cx = cx / total_mass; cy = cy / total_mass; cz = cz / total_mass; total_mass = 0.0; vx = 0.0; vy = 0.0; vz = 0.0; for (i=0; i < state->atoms; i++) { if (corrInfo->mask2[i]) { vx += state->masses[i] * x[i]; vy += state->masses[i] * y[i]; vz += state->masses[i] * z[i]; total_mass += state->masses[i]; } } vx = vx / total_mass; vy = vy / total_mass; vz = vz / total_mass; corrInfo->vx[corrInfo->frame] = vx - cx; corrInfo->vy[corrInfo->frame] = vy - cy; corrInfo->vz[corrInfo->frame] = vz - cz; corrInfo->cx[corrInfo->frame] = cx; corrInfo->cy[corrInfo->frame] = cy; corrInfo->cz[corrInfo->frame] = cz; corrInfo->frame++; return 1; } /** ACTION ROUTINE ************************************************************* * * transformDihedral --- compute/store dihedral angles * ******************************************************************************/ int transformDihedral(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { stackType **argumentStackPointer; char *buffer; scalarInfo *info; ptrajState *state; int i; int mask1tot; int mask2tot; int mask3tot; int mask4tot; double cx1, cy1, cz1, total_mass1; double cx2, cy2, cz2, total_mass2; double cx3, cy3, cz3, total_mass3; double cx4, cy4, cz4, total_mass4; FILE *outFile; /* * USAGE: * * dihedral name mask1 mask2 mask3 mask4 [out ] [time ] * * action argument usage: * * darg1: time interval in ps (for output) * carg1: * a scalarInfo structure */ if (mode == PTRAJ_SETUP) { /* * ACTION: PTRAJ_SETUP */ argumentStackPointer = (stackType **) action->carg1; action->carg1 = NULL; /* * set up the information necessary to place this on the scalarStack */ info = (scalarInfo *) safe_malloc(sizeof(scalarInfo)); INITIALIZE_scalarInfo(info); info->mode = SCALAR_TORSION; info->totalFrames = -1; info->name = getArgumentString(argumentStackPointer, NULL); if (info->name == NULL) { fprintf(stderr, "WARNING: ptraj(), dihedral: It is necessary to specify a unique name\n"); fprintf(stderr, "for each angle specified. Ignoring command...\n"); safe_free(info); return -1; } else if ( scalarStackGetName(&scalarStack, info->name) != NULL ) { fprintf(stderr, "WARNING: ptraj(), dihedral: The chosen name (%s) has already been used.\n", info->name); fprintf(stderr, "Ignoring command...\n"); safe_free(info); return -1; } /* * grab the output filename, if specified */ info->filename = argumentStackKeyToString(argumentStackPointer, "out", NULL); /* * push the distance info on to the scalar stack */ pushStack(&scalarStack, (void *) info); action->darg1 = argumentStackKeyToDouble(argumentStackPointer, "time", 1.0); /* * process mask1 --> mask4 */ buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { fprintf(stderr, "WARNING in ptraj(), dihedral: Error in specification of the first mask\n"); fprintf(stderr, "Ignoring command\n"); safe_free(info); return -1; } else { info->mask1 = processAtomMask(buffer, action->state); safe_free(buffer); } buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { fprintf(stderr, "WARNING in ptraj(), dihedral: Error in specification of the second mask\n"); fprintf(stderr, "Ignoring command\n"); safe_free(info); return -1; } else { info->mask2 = processAtomMask(buffer, action->state); safe_free(buffer); } buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { fprintf(stderr, "WARNING in ptraj(), dihedral: Error in specification of the third mask\n"); fprintf(stderr, "Ignoring command\n"); safe_free(info); return -1; } else { info->mask3 = processAtomMask(buffer, action->state); safe_free(buffer); } buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { fprintf(stderr, "WARNING in ptraj(), dihedral: Error in specification of the fourth mask\n"); fprintf(stderr, "Ignoring command\n"); safe_free(info); return -1; } else { info->mask4 = processAtomMask(buffer, action->state); safe_free(buffer); } /* * check to see if each mask only represents a single atom or not * (to save on memory and speed up calculation) */ mask1tot = 0; info->atom1 = -1; mask2tot = 0; info->atom2 = -1; mask3tot = 0; info->atom3 = -1; mask4tot = 0; info->atom4 = -1; for (i=0; i < action->state->atoms; i++) { if (info->mask1[i] == 1) { mask1tot++; info->atom1 = i; } if (info->mask2[i] == 1) { mask2tot++; info->atom2 = i; } if (info->mask3[i] == 1) { mask3tot++; info->atom3 = i; } if (info->mask4[i] == 1) { mask4tot++; info->atom4 = i; } } if (mask1tot == 0) { fprintf(stderr, "WARNING in ptraj(), dihedral: No atoms selected in mask1, ignoring command\n"); safe_free(info->mask1); safe_free(info); return -1; } else if (mask1tot == 1) { safe_free(info->mask1); info->mask1 = NULL; } else info->atom1 = -1; if (mask2tot == 0) { fprintf(stderr, "WARNING in ptraj(), dihedral: No atoms selected in mask2, ignoring command\n"); safe_free(info->mask2); safe_free(info); return -1; } else if (mask2tot == 1) { safe_free(info->mask2); info->mask2 = NULL; } else info->atom2 = -1; if (mask3tot == 0) { fprintf(stderr, "WARNING in ptraj(), dihedral: No atoms selected in mask3, ignoring command\n"); safe_free(info->mask3); safe_free(info); return -1; } else if (mask3tot == 1) { safe_free(info->mask3); info->mask3 = NULL; } else info->atom3 = -1; if (mask4tot == 0) { fprintf(stderr, "WARNING in ptraj(), dihedral: No atoms selected in mask4, ignoring command\n"); safe_free(info->mask4); safe_free(info); return -1; } else if (mask4tot == 1) { safe_free(info->mask4); info->mask4 = NULL; } else info->atom4 = -1; action->carg1 = (void *) info; return 0; } info = (scalarInfo *) action->carg1; if (mode == PTRAJ_STATUS) { /* * ACTION: PTRAJ_STATUS */ fprintf(stdout, " DIHEDRAL: saved to array named %s\n", info->name); if (info->atom1 == -1) { fprintf(stdout, " Atom selection 1 is "); printAtomMask(info->mask1, action->state); } else { fprintf(stdout, " Atom selection 1 is :%i@%s\n", atomToResidue(info->atom1+1, action->state->residues, action->state->ipres), action->state->atomName[info->atom1]); } if (info->atom2 == -1) { fprintf(stdout, " Atom selection 2 is "); printAtomMask(info->mask2, action->state); } else { fprintf(stdout, " Atom selection 2 is :%i@%s\n", atomToResidue(info->atom2+1, action->state->residues, action->state->ipres), action->state->atomName[info->atom2]); } if (info->atom3 == -1) { fprintf(stdout, " Atom selection 3 is "); printAtomMask(info->mask3, action->state); } else { fprintf(stdout, " Atom selection 3 is :%i@%s\n", atomToResidue(info->atom3+1, action->state->residues, action->state->ipres), action->state->atomName[info->atom3]); } if (info->atom4 == -1) { fprintf(stdout, " Atom selection 4 is "); printAtomMask(info->mask4, action->state); } else { fprintf(stdout, " Atom selection 4 is :%i@%s\n", atomToResidue(info->atom4+1, action->state->residues, action->state->ipres), action->state->atomName[info->atom4]); } if (info->filename != NULL) { fprintf(stdout, " Data will be dumped to a file named %s\n", info->filename); } } else if (mode == PTRAJ_PRINT) { /* * ACTION: PTRAJ_PRINT */ if ( openFile(&outFile, info->filename, "w") == 0 ) { fprintf(stdout, "WARNING in ptraj(), dihedral: couldn't open file %s\n", info->filename); return 0; } fprintf(stdout, "PTRAJ DIHEDRAL dumping named values %s\n", info->name); for (i=0; i < action->state->maxFrames; i++) { fprintf(outFile, "%10.2f %f\n", (i+1) * action->darg1, info->value[i]); } safe_fclose(outFile); } else if (mode == PTRAJ_CLEANUP) { /* * ACTION: PTRAJ_CLEANUP */ safe_free(info->name); safe_free(info->filename); safe_free(info->mask1); safe_free(info->mask2); safe_free(info->mask3); safe_free(info->mask4); safe_free(info->value); INITIALIZE_scalarInfo(info); safe_free(info); } if (mode != PTRAJ_ACTION) return 0; /* * ACTION: PTRAJ_ACTION */ state = (ptrajState *) action->state; /* * update local state information */ for (i=0; i<6; i++) state->box[i] = box[i]; if (info->totalFrames < 0) { info->totalFrames = state->maxFrames; info->value = (double *) safe_malloc(sizeof(double) * info->totalFrames); } if (info->frame > info->totalFrames) { warning("transformDihedral()", "Blowing array; too many frames!!\n"); return 0; } cx1 = 0.0; cy1 = 0.0; cz1 = 0.0; total_mass1 = 0.0; cx2 = 0.0; cy2 = 0.0; cz2 = 0.0; total_mass2 = 0.0; cx3 = 0.0; cy3 = 0.0; cz3 = 0.0; total_mass3 = 0.0; cx4 = 0.0; cy4 = 0.0; cz4 = 0.0; total_mass4 = 0.0; if (info->atom1 == -1) { for (i=0; i < state->atoms; i++) { if (info->mask1[i]) { cx1 += state->masses[i] * x[i]; cy1 += state->masses[i] * y[i]; cz1 += state->masses[i] * z[i]; total_mass1 += state->masses[i]; } } cx1 = cx1 / total_mass1; cy1 = cy1 / total_mass1; cz1 = cz1 / total_mass1; } else { cx1 = x[info->atom1]; cy1 = y[info->atom1]; cz1 = z[info->atom1]; } if (info->atom2 == -1) { for (i=0; i < state->atoms; i++) { if (info->mask2[i]) { cx2 += state->masses[i] * x[i]; cy2 += state->masses[i] * y[i]; cz2 += state->masses[i] * z[i]; total_mass2 += state->masses[i]; } } cx2 = cx2 / total_mass2; cy2 = cy2 / total_mass2; cz2 = cz2 / total_mass2; } else { cx2 = x[info->atom2]; cy2 = y[info->atom2]; cz2 = z[info->atom2]; } if (info->atom3 == -1) { for (i=0; i < state->atoms; i++) { if (info->mask3[i]) { cx3 += state->masses[i] * x[i]; cy3 += state->masses[i] * y[i]; cz3 += state->masses[i] * z[i]; total_mass3 += state->masses[i]; } } cx3 = cx3 / total_mass3; cy3 = cy3 / total_mass3; cz3 = cz3 / total_mass3; } else { cx3 = x[info->atom3]; cy3 = y[info->atom3]; cz3 = z[info->atom3]; } if (info->atom4 == -1) { for (i=0; i < state->atoms; i++) { if (info->mask4[i]) { cx4 += state->masses[i] * x[i]; cy4 += state->masses[i] * y[i]; cz4 += state->masses[i] * z[i]; total_mass4 += state->masses[i]; } } cx4 = cx4 / total_mass4; cy4 = cy4 / total_mass4; cz4 = cz4 / total_mass4; } else { cx4 = x[info->atom4]; cy4 = y[info->atom4]; cz4 = z[info->atom4]; } info->value[info->frame] = torsion(cx1,cy1,cz1,cx2,cy2,cz2,cx3,cy3,cz3,cx4,cy4,cz4); info->frame++; return 1; } /** ACTION ROUTINE ************************************************************* * * transformDiffusion() --- calculate mean squared displacements vs. time * ******************************************************************************/ int transformDiffusion(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { stackType **argumentStackPointer; char *buffer; ptrajState *state; transformDiffusionInfo *diffusionInfo; int i; int currentAtom; double delx, dely, delz; double xx, yy, zz; double avgx, avgy, avgz; double average; double time; /* * USAGE: * * diffusion mask [average] [time