/* _______________________________________________________________________ * * RDPARM/PTRAJ: 2008 * _______________________________________________________________________ * * This file is part of rdparm/ptraj. * * rdparm/ptraj is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * rdparm/ptraj is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You can receive a copy of the GNU General Public License from * http://www.gnu.org or by writing to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA * ________________________________________________________________________ * * CVS tracking: * * $Header: /home/case/cvsroot/amber11/AmberTools/src/ptraj/actions.c,v 10.11 2010/02/25 13:29:54 droe Exp $ * * Revision: $Revision: 10.11 $ * Date: $Date: 2010/02/25 13:29:54 $ * Last checked in by $Author: droe $ * ________________________________________________________________________ * * * CONTACT INFO: To learn who the code developers are, who to contact for * more information, and to know what version of the code this is, refer * to the CVS information and the include files (contributors.h && version.h) * */ #include "contributors.h" #include "version.h" /* ________________________________________________________________________ */ #define ACTION_MODULE #include "ptraj.h" #include #include #include /* * 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 * transformContacts --- calculate number of nearest neighbors of an atom (Holger Gohlke, JWGU) * transformCorr --- perform correlation analysis (Vickie Tsui, Scripps) * transformDihedral --- calculate torsion/dihedral angles * transformDihedralCluster --- calculate torsion/dihedral angles (Dan Roe) * transformDiffusion --- calculate mean squared displacements vs. time * transformDipole --- bin dipoles (Jed Pitera, UCSF) * transformDistance --- compute/store (imaged) distances * transformDNAiontracker --- track ions a' la Hambelberg/Wilson/Williams * transformGrid --- grid atomic densities * transformHBond --- calculate hydrogen bond information * transformImage --- Image molecules outside of a periodic box back in * transformMatrix --- calculate matrices of covariances, cross-correlations, distances (Gohlke) * transformPrincipal --- align coordinates along the principal axis * transformPucker --- compute/store pucker values * transformRadial --- compute radial distribution functions * transformRadiusOfGyration --- compute radius of gyration (Holger Gohlke, JWGU) * transformRMS --- perform RMS fitting * transformScale --- scale the coordinates by a specified amount * transformSecondaryStructure --- determine secondary structure following the method by Kabsch & Sander (Gohlke) * transformStrip --- strip coordinates * transformTranslate --- shift the coordinates by a specified amount * transformTruncOct --- trim/orient a box to make it a truncated octahedron * transformUnwrap --- unwrap atoms (InSuk Joung, Rutgers) * transformVector --- compute/store various vector quantities (IRED, CORR, CORRIRED vectors: Gohlke) * 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) { char *name = "test"; 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 (actionInformation *) 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 * darg3 -- double argument 3 * darg4 -- double argument 4 * carg1 -- (void *) argument 1 (i.e. a pointer to a structure) * carg2 -- (void *) argument 2 * carg3 -- (void *) argument 3 * carg4 -- (void *) argument 4 * * 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. */ #ifdef MPI printParallelError(name); return -1; #endif 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. * * To store information for later processing, some actions use stacks. * For example transformAngle makes use of the scalarStack, which holds * scalarInfo structures. The user passes in a name for the angle command, * and this name is given to the scalarInfo structure, which is then placed * on the scalarStack. The scalarInfo structure can be retrieved from the * stack, which is done when the angle info is added, and when the angles * are printed. */ return 1; } /** ACTION ROUTINE ************************************************************* * * transformAngle --- compute/store angles * ******************************************************************************/ int transformAngle(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { char *name = "angle"; stackType **argumentStackPointer; char *buffer, buffer2[BUFFER_SIZE]; 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; void *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 */ #ifdef MPI #endif 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(stdout, "WARNING: ptraj(), angle: It is necessary to specify a unique name for\n"); fprintf(stdout, "each angle specified. Ignoring command...\n"); safe_free(info); return -1; } else if ( scalarStackGetName(&scalarStack, info->name) != NULL ) { fprintf(stdout, "WARNING: ptraj(), angle: The chosen name (%s) has already been used.\n", info->name); fprintf(stdout, "Ignoring command...\n"); safe_free(info); return -1; } info->state = action->state; /* * grab the output filename, if specified */ info->filename = argumentStackKeyToString(argumentStackPointer, "out", NULL); /* * push this entry onto the scalarStack */ pushBottomStack(&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(stdout, "WARNING in ptraj(), angle: Error in specification of the first mask\n"); fprintf(stdout, "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(stdout, "WARNING in ptraj(), angle: Error in specification of the second mask\n"); fprintf(stdout, "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(stdout, "WARNING in ptraj(), angle: Error in specification of the third mask\n"); fprintf(stdout, "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(stdout, "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(stdout, "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(stdout, "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(stdout, info->mask1, action->state); fprintf(stdout, "\n"); } 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(stdout, info->mask2, action->state); fprintf(stdout, "\n"); } 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(stdout, info->mask3, action->state); fprintf(stdout, "\n"); } 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 ( info->filename != NULL) { outFile = ptrajOpenW(info->filename); if ( outFile == NULL ) { 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)/worldsize; i++) { ptrajfprintf(outFile, "%10.2f %f\n", (i*worldsize+worldrank+1)*action->darg1, info->value[i]); } ptrajCloseFile(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) { char *name = "atomicfluct"; 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 */ #ifdef MPI #endif 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 = worldrank + 1; 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(stdout, action->mask, action->state); fprintf(stdout, "\n"); } else if (mode == PTRAJ_PRINT) { /* * ACTION: PTRAJ_PRINT */ /* * PF - multiptraj * We need to combine the data using 3 MPI_Reduce's, * passing them to rank 0, then let that rank print. * First need to allocate memory for rank 0 recvbuf */ int sets; #ifdef MPI double *xx_local, *yy_local, *zz_local; xx_local = (double *) action->carg2; yy_local = (double *) action->carg3; zz_local = (double *) action->carg4; if (worldrank == 0) { 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); } MPI_Reduce(xx_local, xx, action->state->atoms * 2, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); MPI_Reduce(yy_local, yy, action->state->atoms * 2, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); MPI_Reduce(zz_local, zz, action->state->atoms * 2, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); MPI_Reduce(&action->iarg2, &sets, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); #else xx = (double *) action->carg2; yy = (double *) action->carg3; zz = (double *) action->carg4; sets = action->iarg2; #endif if (worldrank == 0) { 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; for (i=0; i < action->state->atoms; i++) { xx[i] /= sets; yy[i] /= sets; zz[i] /= sets; } for (i=0; i < action->state->atoms; i++) { xi = xx[i+action->state->atoms]/sets - xx[i]*xx[i]; yi = yy[i+action->state->atoms]/sets - yy[i]*yy[i]; zi = zz[i+action->state->atoms]/sets - 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) { if (action->iarg4 > 0) /* * B-factors are (8/3)*PI*PI * **2 hence we do not sqrt the fluctuations!!! */ results[i] = bfactor * fluct; else results[i] = 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 += 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 += 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; } #ifdef MPI safe_free(xx); safe_free(yy); safe_free(zz); #endif } } 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; 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]; } } action->iarg3 += worldsize; return 1; } windowInfoType * setupWindowInfo(windowInfoType *windowInfo, int windowWidth, /* width between the windows */ int windowOffset, /* offset between windows */ int windows, /* number of windows to create or increase */ int size) /* amount of space to allocate per window */ { int i, sz; /* * creates/reallocates a windowInfoType structure to running average * information over multiple windows of width "windowWidth" with an * offset between windows of "windowOffset". * * If windowInfo != NULL reallocate the structures to INCREASE * the size by "windows" * * If windowInfo == NULL do the initial allocation */ if (windowInfo == NULL) { if (windowWidth < 1 || windowOffset < 1 || windows < 1) { warning("setupWindowInfo()", "Window offset (%i), width (%i), or #windows (%i) out of bounds!\n", windowOffset, windowWidth, windows); return NULL; } windowInfo = safe_malloc(sizeof(windowInfoType)); windowInfo->windowWidth = windowWidth; windowInfo->windowOffset = windowOffset; windowInfo->windows = windows; windowInfo->windowStart = (int *) safe_malloc( sizeof(int) * windows ); windowInfo->windowStop = (int *) safe_malloc( sizeof(int) * windows ); windowInfo->windowVisits = (int *) safe_malloc( sizeof(int) * windows ); windowInfo->x = (double *) safe_malloc( sizeof(double) * windows * size ); windowInfo->y = (double *) safe_malloc( sizeof(double) * windows * size ); windowInfo->z = (double *) safe_malloc( sizeof(double) * windows * size ); for (i=0; i < windows; i++) { windowInfo->windowStart[i] = i * windowOffset + 1; windowInfo->windowStop[i] = windowInfo->windowStart[i] + windowWidth; if (prnlev > 3) { fprintf(stdout, "Window %i from %i to %i\n", i+1, windowInfo->windowStart[i], windowInfo->windowStop[i]); } } return windowInfo; } sz = sizeof(int); windowInfo->windowStart = (int *) safe_realloc( (void *) windowInfo->windowStart, sz*windowInfo->windows, sz*windows ); windowInfo->windowStop = (int *) safe_realloc( (void *) windowInfo->windowStop, sz*windowInfo->windows, sz*windows ); windowInfo->windowVisits = (int *) safe_realloc( (void *) windowInfo->windowVisits, sz*windowInfo->windows, sz*windows ); sz = sizeof(double) * size; windowInfo->x = (double *) safe_realloc( (void *) windowInfo->x, sz*windowInfo->windows, sz*windows ); windowInfo->y = (double *) safe_realloc( (void *) windowInfo->y, sz*windowInfo->windows, sz*windows ); windowInfo->z = (double *) safe_realloc( (void *) windowInfo->z, sz*windowInfo->windows, sz*windows ); for (i=0; i < windows; i++) { windowInfo->windowStart[i+windowInfo->windows] = windowInfo->windowStart[windowInfo->windows-1] + (i+1)*windowInfo->windowOffset; windowInfo->windowStop[i+windowInfo->windows] = windowInfo->windowStart[i+windowInfo->windows] + windowInfo->windowWidth; windowInfo->windowVisits[i+windowInfo->windows] = 0; if (prnlev > 3) { fprintf(stdout, "Window %i from %i to %i\n", i+windowInfo->windows+1, windowInfo->windowStart[i+windowInfo->windows], windowInfo->windowStop[i+windowInfo->windows]); } } windowInfo->windows = windowInfo->windows + windows; return(windowInfo); } /** ACTION ROUTINE ************************************************************* * * transformAtomicFluct3D() --- compute atomic positional fluctuations in 3D * ******************************************************************************/ int transformAtomicFluct3D(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { char *name = "atomicfluct3D"; stackType **argumentStackPointer; char *buffer; char *filename; coordinateInfo *info; double *xx, *yy, *zz, *results; double xi, yi, zi, fluct, bfactor; int i, j, w, ind, minwin, maxwin, windowWidth, windowOffset; windowInfoType *windowInfo; /* * USAGE: * * atomicfluct [out filename] [] [start ] [stop ] [offset ] * [byres | byatom | bymask] [bfactor] [window [by ]] * * 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 */ #ifdef MPI printParallelError(name); return -1; #endif 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->carg1 = (void *) info; 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); } /* * goal here is to create a 3D bfactor value that represents a series of * moving windows in time. The width of each window is specified by "window" * and the offset for each window is specified by "by". For example, if you * specify "window 1000 by 100" this will create averages 0-1000, 100-1100, etc. */ windowWidth = argumentStackKeyToInteger(argumentStackPointer, "window", 1000); windowOffset = argumentStackKeyToInteger(argumentStackPointer, "by", 500); windowInfo = setupWindowInfo(NULL, windowWidth, windowOffset, 100, 2*action->state->atoms); action->carg2 = (void *) windowInfo; action->iarg2 = 0; action->iarg3 = 0; return 0; } else if (mode == PTRAJ_STATUS) { /* * ACTION: PTRAJ_STATUS */ info = (coordinateInfo *) action->carg1; windowInfo = (windowInfoType *) action->carg2; fprintf(stdout, " ATOMICFLUCT -3D- : 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"); fprintf(stdout, " WINDOW width = %i, offset = %i\n", windowInfo->windowWidth, windowInfo->windowOffset); 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(stdout, action->mask, action->state); fprintf(stdout, "\n"); } else if (mode == PTRAJ_PRINT) { /* * ACTION: PTRAJ_PRINT */ results = (double *) safe_malloc(sizeof(double) * action->state->atoms); info = (coordinateInfo *) action->carg1; windowInfo = (windowInfoType *) action->carg2; 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 (w=0; w < windowInfo->windows; w++) { if (windowInfo->windowVisits[w] > 0) { printf("WINDOW VISITS = %i for window %i\n", windowInfo->windowVisits[w], w); for (i=0; i < action->state->atoms; i++) { ind = i+w*2*action->state->atoms; windowInfo->x[ind] /= windowInfo->windowVisits[w]; windowInfo->y[ind] /= windowInfo->windowVisits[w]; windowInfo->z[ind] /= windowInfo->windowVisits[w]; } for (i=0; i < action->state->atoms; i++) { ind = i+w*2*action->state->atoms; xi = windowInfo->x[ ind + action->state->atoms] / windowInfo->windowVisits[w] - windowInfo->x[ind] * windowInfo->x[ind]; yi = windowInfo->y[ ind + action->state->atoms] / windowInfo->windowVisits[w] - windowInfo->y[ind] * windowInfo->y[ind]; zi = windowInfo->z[ ind + action->state->atoms] / windowInfo->windowVisits[w] - windowInfo->z[ind] * windowInfo->z[ind]; windowInfo->x[ind] = xi; windowInfo->y[ind] = yi; windowInfo->z[ind] = zi; } for (i=0; i < action->state->atoms; i++) { ind = i+w*2*action->state->atoms; fluct = windowInfo->x[ind] + windowInfo->y[ind] + windowInfo->z[ind]; if (fluct > 0) { if (action->iarg4 > 0) /* * B-factors are (8/3)*PI*PI * **2 hence we do not sqrt the fluctuations!!! */ results[i] = bfactor * fluct; else results[i] = 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 %.2f %f\n", j+1, (windowInfo->windowStart[w] + windowInfo->windowWidth/2.0), 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 += results[j] * action->state->masses[j]; } } if (xi > SMALL) fprintf(info->file, " %6i %.2f %f\n", i+1, (windowInfo->windowStart[w] + windowInfo->windowWidth/2.0), 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 += results[i] * action->state->masses[i]; } } if (xi > SMALL) fprintf(info->file, " %6i %.2f %f\n", 1, (windowInfo->windowStart[w] + windowInfo->windowWidth/2.0), 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; windowInfo = (windowInfoType *) action->carg2; safe_free(info); safe_free(windowInfo->windowStart); safe_free(windowInfo->windowStop); safe_free(windowInfo->windowVisits); safe_free(windowInfo->x); safe_free(windowInfo->y); safe_free(windowInfo->z); safe_free(windowInfo); } if (mode != PTRAJ_ACTION) return 0; /* * ACTION: PTRAJ_ACTION */ info = (coordinateInfo *) action->carg1; windowInfo = (windowInfoType *) action->carg2; action->iarg3++; if (action->iarg3 >= info->start && (info->stop < 0 || action->iarg3 <= info->stop) && (action->iarg3 - info->start)%info->offset == 0) { action->iarg2++; minwin = 0; maxwin = windowInfo->windows; if (maxwin > windowInfo->windows) { setupWindowInfo(windowInfo, -1, -1, 100, 2*action->state->atoms); maxwin = windowInfo->windows; } if (action->iarg2 > windowInfo->windowStop[minwin]) { minwin++; maxwin++; if (maxwin > windowInfo->windows) maxwin = windowInfo->windows; } for (j=minwin; j < maxwin; j++) { if ( action->iarg2 >= windowInfo->windowStart[j] && action->iarg2 < windowInfo->windowStop[j] ) { for (i=0; i < action->state->atoms; i++) { ind=i+j*2*action->state->atoms; windowInfo->x[ind] += x[i]; windowInfo->y[ind] += y[i]; windowInfo->z[ind] += z[i]; windowInfo->x[ind+action->state->atoms] += x[i]*x[i]; windowInfo->y[ind+action->state->atoms] += y[i]*y[i]; windowInfo->z[ind+action->state->atoms] += z[i]*z[i]; } windowInfo->windowVisits[j] += 1; } } } 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) { char *name = "average"; stackType **argumentStackPointer; char *buffer; char *filename; coordinateInfo *info; double *xx, *yy, *zz, *tmp_xx, *tmp_yy, *tmp_zz; double fluctx, flucty, fluctz; int i, atoms, total; ptrajState *state; /* * USAGE: * * average filename [] [append] [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 */ #ifdef MPI #endif 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); info->append= argumentStackContains(argumentStackPointer, "append"); info->isBox = action->state->IFBOX; 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 { /* #ifdef MPI info->type = COORD_PDB; printf("Cannot print average to AMBER trajectory, defaulting to pdb.\n"); #else info->type = COORD_AMBER_TRAJECTORY; #endif */ 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 = worldrank + 1; 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: %s the %s of the coordinates to file %s\n", (info->append ? "appending" : "dumping"), (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(stdout, action->mask, action->state); fprintf(stdout, "\n"); fprintf(stdout, " Output file information:"); printCoordinateInfo( (void *) info ); } else if (mode == PTRAJ_PRINT) { /* * ACTION: PTRAJ_PRINT */ /* * PF - multiptraj * We need to combine the data using 3 MPI_Reduce's, * passing them to rank 0, then let that rank print. * First need to allocate memory for rank 0 recvbuf */ int sets; double allbox[6]; #ifdef MPI double *xx_local, *yy_local, *zz_local; xx_local = (double *) action->carg2; yy_local = (double *) action->carg3; zz_local = (double *) action->carg4; if (action->iarg1) total = action->state->atoms * 2; else total = action->state->atoms; if (worldrank == 0) { xx = (double *) safe_malloc(sizeof(double) * total); yy = (double *) safe_malloc(sizeof(double) * total); zz = (double *) safe_malloc(sizeof(double) * total); } MPI_Reduce(xx_local, xx, total, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); MPI_Reduce(yy_local, yy, total, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); MPI_Reduce(zz_local, zz, total, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD); MPI_Reduce(&action->iarg2, &sets, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); if (action->iarg1 == 0) { MPI_Reduce(action->state->box, allbox, 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); } #else xx = (double *) action->carg2; yy = (double *) action->carg3; zz = (double *) action->carg4; sets = action->iarg2; for (i = 0; i < 6; i++) allbox[i] = action->state->box[i]; #endif if (worldrank == 0) { info = (coordinateInfo *) action->carg1; if (info->start > action->state->maxFrames) return 0; for (i=0; i < action->state->atoms; i++) { xx[i] /= sets; yy[i] /= sets; zz[i] /= sets; } if (action->iarg1 == 0) for (i=0; i<6; i++) allbox[i] /= sets; if (action->iarg1) { for (i=0; i < action->state->atoms; i++) { xx[i+action->state->atoms] /= sets; yy[i+action->state->atoms] /= sets; zz[i+action->state->atoms] /= sets; 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); } } atoms = 0; for (i=0; i < action->state->atoms; i++) if (action->mask[i]) atoms++; if (atoms == action->state->atoms) { ptrajOutputCoordinates(info, action->state, -1, info->append, 1, 0, action->state->atoms, xx, yy, zz, allbox); } 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(info, state, -1, info->append, 1, 0, state->atoms, tmp_xx, tmp_yy, tmp_zz, allbox); ptrajClearState(&state); safe_free(tmp_xx); safe_free(tmp_yy); safe_free(tmp_zz); } #ifdef MPI safe_free(xx); safe_free(yy); safe_free(zz); #endif } } 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; 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]; } } } action->iarg3 += worldsize; 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) { char *name = "center"; 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 */ #ifdef MPI #endif 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(stdout, action->mask, action->state); fprintf(stdout, "\n"); } 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) { char *name = "checkoverlap"; stackType **argumentStackPointer; char *buffer; ptrajState *state; int i, j; double distance, ucell[9], recip[9]; int *around; /* * USAGE: * * checkoverlap [mask] [min ] [max ] [noimage] [around ] * * action argument usage: * * iarg1: 1 implies don't image * darg1: flag distances lower than this * darg2: flag distances greater than this * carg1: around mask */ if (mode == PTRAJ_SETUP) { /* * ACTION: PTRAJ_SETUP */ #ifdef MPI #endif 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 = argumentStackKeyToString(argumentStackPointer, "around", NULL); if (buffer != NULL) { action->carg1 = processAtomMask(buffer, action->state); safe_free(buffer); } else { action->carg1 = NULL; } 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->carg1 != NULL) { around = (int *) action->carg1; if (action->mask == NULL) { fprintf(stdout, " Checking overlap of atoms: * (all atoms)\n"); } else { fprintf(stdout, " Checking overlap of atoms: "); printAtomMask(stdout, action->mask, action->state); fprintf(stdout, "\n"); } fprintf(stdout, " To atoms in: "); printAtomMask(stdout, around, action->state); fprintf(stdout, "\n"); } else { if (action->mask == NULL) { fprintf(stdout, " Atom selection is: * (all atoms)\n"); } else { fprintf(stdout, " Atom selection is "); printAtomMask(stdout, action->mask, action->state); fprintf(stdout, "\n"); } } } 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(stdout, " 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); if (action->carg1 != NULL) { around = (int *) action->carg1; } else { around = NULL; } for (i = 0; i < state->atoms; i++) { for (j = 0; j < state->atoms; j++) { if ( (j > i && action->mask == NULL && around == NULL) || (j > i && around == NULL && action->mask[i] && action->mask[j]) || (i != j && action->mask != NULL && action->mask[i] && around != NULL && around[j]) ) { distance = calculateDistance2(i, j, x, y, z, box, ucell, recip, 0.0, action->iarg1); if (distance < action->darg1) { fprintf(stdout, "OVERLAP: atoms %5i (", i+1); printAtomCompact(stdout, i, action->state); fprintf(stdout, ") and %5i (", 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 %5i (", i+1); printAtomCompact(stdout, i, action->state); fprintf(stdout, ") and %5i (", j+1); printAtomCompact(stdout, j, action->state); fprintf(stdout, ") are too far apart (%8.3f)!\n", sqrt(distance)); } } } } return 1; } scalarInfo * ptrajCopyScalar(scalarInfo **scalarinp) { scalarInfo *scalar, *scalarin; /* * Make a copy of the input scalar pointer. */ scalarin = *scalarinp; scalar = (scalarInfo*)SafeMalloc(__FILE__, __LINE__, sizeof(scalarInfo)); INITIALIZE_scalarInfo(scalar); scalar->mode = scalarin->mode; scalar->totalFrames = scalarin->totalFrames; scalar->frame = scalarin->frame; scalar->mean = scalarin->mean; scalar->stddev = scalarin->stddev; scalar->max = scalarin->max; scalar->min = scalarin->min; scalar->atom1 = scalarin->atom1; scalar->atom2 = scalarin->atom2; scalar->atom3 = scalarin->atom3; scalar->atom4 = scalarin->atom4; scalar->atom5 = scalarin->atom5; scalar->mask1 = scalarin->mask1; scalar->mask2 = scalarin->mask2; scalar->mask3 = scalarin->mask3; scalar->mask4 = scalarin->mask4; scalar->mask5 = scalarin->mask5; scalar->cos = scalarin->cos; scalar->sin = scalarin->sin; scalar->value = scalarin->value; scalar->results = scalarin->results; scalar->action = scalarin->action; scalar->name = scalarin->name; scalar->filename = scalarin->filename; return scalar; } actionInformation * ptrajCopyAction(actionInformation **actioninp) { actionInformation *action, *actionin; /* * Make a copy of the input action pointer. */ actionin = *actioninp; action = (actionInformation*)SafeMalloc(__FILE__, __LINE__, sizeof(actionInformation)); INITIALIZE_actionInformation(action); action->fxn = actionin->fxn; action->type = actionin->type; action->iarg1 = actionin->iarg1; action->iarg2 = actionin->iarg2; action->iarg3 = actionin->iarg3; action->iarg4 = actionin->iarg4; action->iarg5 = actionin->iarg5; action->iarg6 = actionin->iarg6; action->iarg7 = actionin->iarg7; action->darg1 = actionin->darg1; action->darg2 = actionin->darg2; action->darg3 = actionin->darg3; action->darg4 = actionin->darg4; action->suppressProcessing = actionin->suppressProcessing; action->performSecondPass = actionin->performSecondPass; /* The following members are pointers. Just copy from the old action. So you need to allocate new space for the members. */ action->state = actionin->state; action->mask = actionin->mask; action->carg1 = actionin->carg1; action->carg2 = actionin->carg2; action->carg3 = actionin->carg3; action->carg4 = actionin->carg4; action->carg5 = actionin->carg5; action->carg6 = actionin->carg6; action->carg7 = actionin->carg7; return action; } /** ACTION ROUTINE ************************************************************* * * transformCluster() --- organizes trajectory frames into clusters. * * ******************************************************************************/ /* * To use this action, the user should choose a clustering algorithm, and either * a number of clusters or a value of epsilon governing maxiumum cluster size. * * The available clustering algorithms are: * * Hierarchical: Hierarchical divisive clustering first places all points * in one cluster. The algorithm splits, in each iteration, a cluster into * two subclusters; this continues until either the target number of clusters * is reached, or until the maximum eccentricity of any cluster is less than * epsilon. At each iteration, choose the cluster with the largest eccentricity. * Choose two extrema (points with maximal pairwise distance) from this cluster; * these become the centroids of two new clusters, and points are assigned * to one child cluster or the other based on which centroid they are closest to. * The centroid of each child cluster is then calculated, and points reassigned * to child clusters based on these new centroids. Three centroids are chosen in all. * * Single Linkage: Begin with the original collection of points as your cluster * collection. Choose the two closest points, and form a cluster; remove those * two points from the cluster collection, and insert the centroid of this new * cluster. Repeat the process until all points lie in a cluster. To cluster * two centroids, merge their two clusters into a supercluster with new centroid; * to cluster a point and a centroid, add the point to the centroid's cluster. * If an epsilon is provided: Do not cluster two points if their distance is epsilon * or greater. If a cluster count is provided: Continue until every point is clustered; * then, continue merging clusters as long as the number of clusters exceeds the desired * cluster count. * * k-means: If a cluster count n is specified, assign the first n frames to n initial * clusters. Iterate over all the points. Assign each new point to the cluster whose * ncentroid it is nearest to, and recompute the cluster's centroid. If an epsilon value * is specified: Start with 1 cluster. If a point has distance less than epsilon from a * centroid, add it to an existing cluster; otherwise, form a new cluster with the point. * * Centripetal: Based on the CURE clustering algorithm. * * Cobweb: Based on the COBWEB clustering algorithm. See cluster.h for details. * * Bayesian: Based on the AutoClass clustering algorithm. See cluster.h for details. * * SOM: Self-organizing maps. See cluster.h for details * ******************************************************************************/ #define INVALID_ARGUMENTS_RETURN_CODE -1 /* * Save x, y, and z coordinates of atoms into the arrays in trajInfo. * Called once for each frame, normally. Only accumulates atoms according to the action mask. * If FullTrajInfo is not null, it gets *all* the atom information. */ void AccumulateCoordinates(trajectoryInfo* trajInfo, actionInformation* action, double* x, double* y, double* z, trajectoryInfo* FullTrajInfo) { int i; int ValidIndex = 0; int NumberOfFrames; int Sieve = 0; int FrameIndex; PtrajClustering* Clustering; /* * Special case: If we're sieving through the coordinates in a first pass, then * skip over all but every nth frame */ if (action->type == TRANSFORM_CLUSTER) { Clustering = (PtrajClustering*)action->carg5; /* * iarg3 is the sieve size. * iarg4 is the number of frames we'll process in our sieve */ Sieve = action->iarg3; if (!action->iarg4) { if (Sieve> 0) { PtrajClusteringChooseSievePoints(Clustering, action->state->maxFrames); } else { action->iarg4 = action->state->maxFrames; Clustering->PointCount = action->iarg4; } } NumberOfFrames = action->iarg4; } else { NumberOfFrames = action->state->maxFrames; } /* * allocate memory if this is the first visit */ if (trajInfo->x == NULL) { trajInfo->x = (float *) safe_malloc(sizeof(float) * trajInfo->atoms * NumberOfFrames); memset(trajInfo->x, 0, sizeof(float) * trajInfo->atoms * NumberOfFrames); trajInfo->y = (float *) safe_malloc(sizeof(float) * trajInfo->atoms * NumberOfFrames); memset(trajInfo->y, 0, sizeof(float) * trajInfo->atoms * NumberOfFrames); trajInfo->z = (float *) safe_malloc(sizeof(float) * trajInfo->atoms * NumberOfFrames); memset(trajInfo->z, 0, sizeof(float) * trajInfo->atoms * NumberOfFrames); trajInfo->allocated = action->state->maxFrames; trajInfo->current = 0; } if (FullTrajInfo && FullTrajInfo->x==NULL) { FullTrajInfo->x = (float *) safe_malloc(sizeof(float) * FullTrajInfo->atoms * NumberOfFrames); memset(FullTrajInfo->x, 0, sizeof(float) * FullTrajInfo->atoms * NumberOfFrames); FullTrajInfo->y = (float *) safe_malloc(sizeof(float) * FullTrajInfo->atoms * NumberOfFrames); memset(FullTrajInfo->y, 0, sizeof(float) * FullTrajInfo->atoms * NumberOfFrames); FullTrajInfo->z = (float *) safe_malloc(sizeof(float) * FullTrajInfo->atoms * NumberOfFrames); memset(FullTrajInfo->z, 0, sizeof(float) * FullTrajInfo->atoms * NumberOfFrames); FullTrajInfo->allocated = action->state->maxFrames; FullTrajInfo->current = 0; } /* * accumulate coordinates for the current frame - unless we're sieving them out */ if (Sieve==0 || Clustering->SieveFirstPass[trajInfo->current]) { /*fprintf(stdout, "In the first pass, do frame %d\n", trajInfo->current);*/ FrameIndex = trajInfo->current; if (Sieve) { FrameIndex = Clustering->SievedIndex; /*trajInfo->current / Sieve;*/ Clustering->SievedIndex++; } ValidIndex = 0; for (i = 0; i < action->state->atoms; i++) { if (action->mask[i]) { trajInfo->x[FrameIndex * trajInfo->atoms + ValidIndex] = (float)x[i]; trajInfo->y[FrameIndex * trajInfo->atoms + ValidIndex] = (float)y[i]; trajInfo->z[FrameIndex * trajInfo->atoms + ValidIndex] = (float)z[i]; ValidIndex ++; } if (FullTrajInfo) { FullTrajInfo->x[FrameIndex * FullTrajInfo->atoms + i] = (float)x[i]; FullTrajInfo->y[FrameIndex * FullTrajInfo->atoms + i] = (float)y[i]; FullTrajInfo->z[FrameIndex * FullTrajInfo->atoms + i] = (float)z[i]; } } } trajInfo->current++; if (FullTrajInfo) { FullTrajInfo->current++; } } /* * GetOutputFileFormat is a helper routine that parses a trajectory file format from the stack */ int GetOutputFileFormat(char* argument) { if (stringMatch(argument, "pdb")) { return COORD_PDB; } else if (stringMatch(argument, "rest")) { return COORD_AMBER_RESTART; } else if (stringMatch(argument, "binpos")) { return COORD_BINPOS; } else if (stringMatch(argument, "amber")) { return COORD_AMBER_TRAJECTORY; } else if (stringMatch(argument, "none")) { return CLUSTER_OUTPUT_NONE; } else { fprintf(stderr, "WARNING: invalid file format %s!\n", argument); return CLUSTER_OUTPUT_NONE; } } void PrintClusterFileInfo(int FileFormat,char* OutputType) { char* FileFormatString = "No output written"; if (FileFormat == COORD_AMBER_TRAJECTORY) { FileFormatString = "Amber trajectory"; } if (FileFormat == COORD_AMBER_RESTART) { FileFormatString = "Amber restart"; } if (FileFormat == COORD_PDB) { FileFormatString = "PDB"; } if (FileFormat == COORD_BINPOS) { FileFormatString = "Binpos"; } if (FileFormat == COORD_CHARMM_TRAJECTORY) { FileFormatString = "CHARMM trajectory"; } fprintf(stdout," Output of type '%s' will be written to a file of type: '%s'\n",OutputType,FileFormatString); } /* * Wrapper for ClusteringOutputClusterInfo; also FREES the clustering */ /* Wrapper for ClusteringOutputClusterInfo; also FREES the clustering */ void PtrajPrintClustering(actionInformation* action) { char FilePath[512]; int ElapsedTime; PtrajClustering* pClustering; pClustering = (PtrajClustering*)action->carg5; ElapsedTime = time(NULL) - pClustering->StartTime; PtrajClusteringOutputClusterInfo(pClustering,action); ClusteringFree( (Clustering *) pClustering); } /* Check cached pairwise distance file. Check the 9 distance pairs between frames 1, 2, 3 and frames n-2, n-1, n. If the difference between the freshly calcuated distance and the distance read from the is not within == 0, need to calcute or recalculate. != 0, file seems OK, return the Matrix. */ SymmetricMatrix* CheckPairwiseDistance(PtrajClustering* This, trajectoryInfo* trajInfo, actionInformation* action, char *CacheDistanceFile) { int FrameIndex; int OtherFrameIndex; int AtomIndex; int OtherAtomIndex; float* FrameX; float* FrameY; float* FrameZ; float* OtherFrameX; float* OtherFrameY; float* OtherFrameZ; float Distance; float AtomDistance; float OtherAtomDistance; SymmetricMatrix* PairwiseDistance; int FrameCount; FILE* DistanceFile = NULL; int SymMatrixReady = 0; FrameCount = action->iarg4; PairwiseDistance = AllocateSymmetricMatrix(FrameCount); /* If we have a distance cached on disk already, use it! */ DistanceFile = fopen(CacheDistanceFile,"r"); if (DistanceFile) { fclose(DistanceFile); fprintf(stdout, " Read in the existing PairwiseDistances file\n"); SymMatrixReady = ReadSymMatrix(CacheDistanceFile, PairwiseDistance); if (SymMatrixReady) { /* To test if the PairwiseDistance is reusable by comparing 4 values in the table. It does not garrant to be correct!! */ if (FrameCount < 6) fprintf(stderr, "Too few frames.\n"); for (FrameIndex=0; FrameIndex<3; FrameIndex++) { FrameX = trajInfo->x + trajInfo->atoms*FrameIndex; FrameY = trajInfo->y + trajInfo->atoms*FrameIndex; FrameZ = trajInfo->z + trajInfo->atoms*FrameIndex; for (OtherFrameIndex=FrameCount-3; OtherFrameIndexx + trajInfo->atoms*OtherFrameIndex; OtherFrameY = trajInfo->y + trajInfo->atoms*OtherFrameIndex; OtherFrameZ = trajInfo->z + trajInfo->atoms*OtherFrameIndex; Distance=PtrajGetDistance(This, FrameIndex, OtherFrameIndex, FrameX, FrameY, FrameZ, OtherFrameX, OtherFrameY, OtherFrameZ); double error = abs(Distance - GetSymElement(PairwiseDistance,FrameIndex,OtherFrameIndex)) / Distance; if (error > 0.000001){ SymMatrixReady = 0; fprintf(stdout, " Error: distance between %d and %d does not match. difference: %.2f%%. \n", FrameIndex, OtherFrameIndex, error*100); return(NULL); } } } } } if (SymMatrixReady) { This->PairwiseDistances = PairwiseDistance; This->PointCount = PairwiseDistance->Size; return PairwiseDistance; } return(NULL); } #define REPRESENTATIVE_COUNT 5 void PtrajPerformClustering(trajectoryInfo* trajInfo,actionInformation* action) { PtrajClustering* pClustering; int* SeedPoints; int DesiredClusterCount; float Epsilon; char FilePath[1024]; FILE* File; /* */ pClustering = (PtrajClustering*)action->carg5; pClustering->StartTime = time(NULL); pClustering->PointCount = action->iarg4; pClustering->PairwiseDistances = CheckPairwiseDistance(pClustering, trajInfo, action, "PairwiseDistances"); /* * iarg5 is the requested distance metric */ if ( action->iarg1 == CLUSTER_LINKAGE || action->iarg1 == CLUSTER_CENTRIPETAL || action->iarg1 == CLUSTER_EDGELINK || action->iarg1 == CLUSTER_COMPLETELINK || action->iarg1 == CLUSTER_HIERARCHICAL || action->iarg1 == CLUSTER_CENTRIPETAL_COMPLETE || action->iarg1 == CLUSTER_AVERAGELINK || action->iarg1 == CLUSTER_ANOVA || action->iarg1 == CLUSTER_AVERAGELINK_FP || action->iarg1 == CLUSTER_LINKAGE_FP || action->iarg1 == CLUSTER_EDGELINK_FP || action->iarg1 == CLUSTER_COMPLETELINK_FP || action->iarg1 == CLUSTER_SOM2 ) /* action->iarg1 == CLUSTER_READMERGE || action->iarg1 == CLUSTER_DBI || action->iarg1 == CLUSTER_COMPARE || action->iarg1 == CLUSTER_BAYESIAN || action->iarg1 == CLUSTER_COBWEB || action->iarg1 == CLUSTER_COBWEB2 || action->iarg1 == CLUSTER_SOM || action->iarg1 == CLUSTER_SOM2 || action->iarg1 == CLUSTER_BYCHANCE || action->iarg1 == CLUSTER_DECOY || action->iarg1 == CLUSTER_MEANS || action->iarg1 == CLUSTER_MEANS2 || action->iarg1 == CLUSTER_MEANS3 */ { /* * The pairwise Distance matrix is not needed with the SOM, Cobweb and Bayesian algorithms, * but is needed when calculating pSF and dBI. Will optimize later. */ /*pClustering->PairwiseDistances = CheckPairwiseDistance(pClustering, trajInfo, action, "PairwiseDistances");*/ if (! pClustering->PairwiseDistances) { remove("PairwiseDistances"); /* avoid changing the PairwiseDistances file in case that it is a link. */ pClustering->PairwiseDistances = ComputePairwiseDistance(pClustering, trajInfo, action, "PairwiseDistances"); } } if (action->iarg1 == CLUSTER_SOM || action->iarg1 == CLUSTER_SOM2 || action->iarg1 == CLUSTER_COBWEB || action->iarg1 == CLUSTER_COBWEB2 || action->iarg1 == CLUSTER_BAYESIAN) { if (attributeArray != NULL || attributeArrayTorsion != NULL) { pClustering->attributeArray = attributeArray; pClustering->attributeArrayTorsion = attributeArrayTorsion; } if (attributeArray != NULL && attributeArray->length == 0) { fprintf(stderr, "No attribute is detected. Exit ptraj clustering.\n"); fflush(stderr); fflush(stdout); exit(1); } } PtrajClusteringBuildTotalCluster(pClustering); pClustering->MatrixTime = time(NULL); DesiredClusterCount = PtrajClusteringGetDesiredClusterCount(pClustering); Epsilon = PtrajClusteringGetDesiredEpsilon(pClustering); /* * Note: *Hackiness* If you are using epsilon, call PtrajClusteringFoo and not * the base ClusteringFoo; otherwise, the epsilon will get lost along the way. */ /* * Seed random number generator */ srand(time(NULL) + clock()); switch ( action->iarg1 ) { case CLUSTER_HIERARCHICAL: fprintf(stdout, " Clustering using the hierarchical algorithm.\n"); PtrajClusteringClusterHierarchical(pClustering,DesiredClusterCount,Epsilon); break; case CLUSTER_LINKAGE: fprintf(stdout, " Clustering using the single linkage algorithm.\n"); PtrajClusteringClusterLinkage(pClustering,DesiredClusterCount,Epsilon); break; case CLUSTER_LINKAGE_FP: fprintf(stdout, " Clustering using the linkage finger print algorithm.\n"); PtrajClusteringClusterLinkageFingerPrint(pClustering,DesiredClusterCount,Epsilon); break; case CLUSTER_EDGELINK: fprintf(stdout, " Clustering using the edge link algorithm.\n"); PtrajClusteringEdgeLinkage(pClustering,DesiredClusterCount,Epsilon); break; case CLUSTER_EDGELINK_FP: fprintf(stdout, " Clustering using the edge link finger print algorithm.\n"); PtrajClusteringEdgeLinkageFingerPrint(pClustering,DesiredClusterCount,Epsilon); break; case CLUSTER_COMPLETELINK: fprintf(stdout, " Clustering using the complete linkage algorithm.\n"); PtrajClusteringCompleteLinkage(pClustering,DesiredClusterCount,Epsilon); break; case CLUSTER_COMPLETELINK_FP: fprintf(stdout, " Clustering using the complete linkage finger print algorithm.\n"); PtrajClusteringCompleteLinkageFingerPrint(pClustering,DesiredClusterCount,Epsilon); break; case CLUSTER_AVERAGELINK: fprintf(stdout, " Clustering using the average linkage algorithm.\n"); PtrajClusteringAverageLinkage(pClustering,DesiredClusterCount,Epsilon); break; case CLUSTER_AVERAGELINK_FP: fprintf(stdout, " Clustering using the average linkage finger print algorithm.\n"); PtrajClusteringAverageLinkageFingerPrint(pClustering,DesiredClusterCount,Epsilon); break; case CLUSTER_MEANS: fprintf(stdout, " Clustering using the means algorithm.\n"); SeedPoints = ClusteringFindKmeansSeeds( (Clustering *) pClustering,DesiredClusterCount); ClusteringClusterMeans( (Clustering *) pClustering,SeedPoints,DesiredClusterCount,action->iarg7, 1); safe_free(SeedPoints); break; case CLUSTER_MEANS2: fprintf(stdout, " Clustering using the means \"2\" algorithm.\n"); SeedPoints = ClusteringFindKmeansSeeds( (Clustering *) pClustering,DesiredClusterCount); ClusteringClusterMeans( (Clustering *) pClustering,SeedPoints,DesiredClusterCount,action->iarg7, 2); safe_free(SeedPoints); break; case CLUSTER_MEANS3: fprintf(stdout, " Clustering using the means \"3\" algorithm.\n"); SeedPoints = ClusteringFindKmeansSeeds( (Clustering *) pClustering,DesiredClusterCount); PtrajClusteringClusterDecoy( pClustering,SeedPoints,DesiredClusterCount,action->iarg7); safe_free(SeedPoints); break; case CLUSTER_CENTRIPETAL: fprintf(stdout, " Clustering using the centripetal algorithm.\n"); PtrajClusteringClusterCentripetal(pClustering,DesiredClusterCount,Epsilon,REPRESENTATIVE_COUNT); break; case CLUSTER_CENTRIPETAL_COMPLETE: fprintf(stdout, " Clustering using the centripetal complete algorithm.\n"); PtrajClusteringClusterCentripetalComplete(pClustering,DesiredClusterCount,Epsilon,REPRESENTATIVE_COUNT); break; case CLUSTER_COBWEB: fprintf(stdout, " Clustering using the CobWeb algorithm.\n"); PtrajClusteringClusterCobweb( (Clustering *) pClustering,DesiredClusterCount,1); break; case CLUSTER_COBWEB2: fprintf(stdout, " Clustering using the CobWeb \"2\" algorithm.\n"); PtrajClusteringClusterCobweb( (Clustering *) pClustering,DesiredClusterCount,2); break; case CLUSTER_BAYESIAN: fprintf(stdout, " Clustering using the Bayesian algorithm.\n"); PtrajClusteringClusterBayesian(pClustering,DesiredClusterCount); break; case CLUSTER_SOM2: fprintf(stdout, " Clustering using the SOM \"2\" algorithm.\n"); PtrajClusteringAverageLinkage(pClustering,DesiredClusterCount,Epsilon); fprintf(stdout, "Performing SOM2 using the SOM Map.\n"); PtrajClusteringClusterSOM2(pClustering,DesiredClusterCount); break; case CLUSTER_SOM: fprintf(stdout, " Clustering using the SOM algorithm.\n"); PtrajClusteringClusterSOM(pClustering,DesiredClusterCount,(int)action->darg4); break; case CLUSTER_DECOY: { trajectoryInfo *DecoyTrajInfo, *trajInfo; DecoyTrajInfo = PtrajClusteringReadDecoy(pClustering,(char*)action->carg7,DesiredClusterCount); DesiredClusterCount = PtrajClusteringGetDesiredClusterCount(pClustering); PtrajClusteringAssignDecoyToCentroid(pClustering, DesiredClusterCount, DecoyTrajInfo); trajInfo = (trajectoryInfo *) action->carg2; PtrajClusteringClusterDecoy(pClustering,NULL,DesiredClusterCount,action->iarg7); } break; case CLUSTER_BYCHANCE: fprintf(stdout, " Clustering using the by-chance algorithm.\n"); PtrajClusteringClusterByChance(pClustering,DesiredClusterCount); break; case CLUSTER_READMERGE: PtrajClusteringReadFromMergeFile(pClustering, ClusterMergingFile); break; case CLUSTER_READTXT: if (!pClustering->PairwiseDistances) pClustering->PairwiseDistances = CheckPairwiseDistance(pClustering, trajInfo, action, "FullPairwiseDistances"); if (!pClustering->PairwiseDistances) pClustering->PairwiseDistances = CheckPairwiseDistance(pClustering, trajInfo, action, "PairwiseDistances"); if (! pClustering->PairwiseDistances) { remove("PairwiseDistances"); /* avoid changing the PairwiseDistances file in case that it is a link. */ } strcpy(FilePath,(char*)action->carg1); strcat(FilePath,".txt"); pClustering->Head = NULL; pClustering->Tail = NULL; PtrajClusteringReadFromDisk(pClustering,FilePath); ClusteringFindAllCentroids( (Clustering *) pClustering); break; case CLUSTER_DBI: if (!pClustering->PairwiseDistances) pClustering->PairwiseDistances = CheckPairwiseDistance(pClustering, trajInfo, action, "FullPairwiseDistances"); if (!pClustering->PairwiseDistances) pClustering->PairwiseDistances = CheckPairwiseDistance(pClustering, trajInfo, action, "PairwiseDistances"); if (! pClustering->PairwiseDistances) { remove("PairwiseDistances"); /* avoid changing the PairwiseDistances file in case that it is a link. */ } strcpy(FilePath,(char*)action->carg1); strcat(FilePath,".txt"); pClustering->Head = NULL; pClustering->Tail = NULL; PtrajClusteringReadFromDisk(pClustering,FilePath); ClusteringFindAllCentroids( (Clustering *) pClustering); File = fopen(FilePath,"a"); fprintf(File, "# Added by CLUSTER_DBI\n"); OutputClusteringStats(pClustering, File); fclose(File); break; case CLUSTER_ANOVA: if (!pClustering->PairwiseDistances) pClustering->PairwiseDistances = CheckPairwiseDistance(pClustering, trajInfo, action, "FullPairwiseDistances"); if (!pClustering->PairwiseDistances) pClustering->PairwiseDistances = CheckPairwiseDistance(pClustering, trajInfo, action, "PairwiseDistances"); if (! pClustering->PairwiseDistances) { remove("PairwiseDistances"); /* avoid changing the PairwiseDistances file in case that it is a link. */ pClustering->PairwiseDistances = ComputePairwiseDistance(pClustering, trajInfo, action, "PairwiseDistances"); } strcpy(FilePath,(char*)action->carg1); strcat(FilePath,".txt"); pClustering->Head = NULL; pClustering->Tail = NULL; PtrajClusteringReadFromDisk(pClustering,FilePath); ClusteringFindAllCentroids( (Clustering *) pClustering); File = fopen(FilePath,"a"); /* fprintf(File, "# Added by CLUSTER_ANOVA\n"); OutputClusteringStats(pClustering, File); */ ClusteringComputeANOVA(pClustering); fclose(File); /* strcpy(FilePath,(char*)action->carg1); strcat(FilePath,"DBI"); strcpy(action->carg1, FilePath); PtrajClusteringOutputClusterInfo(pClustering,action); */ break; case CLUSTER_COMPARE: if (!pClustering->PairwiseDistances) pClustering->PairwiseDistances = CheckPairwiseDistance(pClustering, trajInfo, action, "FullPairwiseDistances"); if (!pClustering->PairwiseDistances) pClustering->PairwiseDistances = CheckPairwiseDistance(pClustering, trajInfo, action, "PairwiseDistances"); if (! pClustering->PairwiseDistances) { remove("PairwiseDistances"); /* avoid changing the PairwiseDistances file in case that it is a link. */ pClustering->PairwiseDistances = ComputePairwiseDistance(pClustering, trajInfo, action, "PairwiseDistances"); } if ((!((char**)action->carg4)[0]) || (!((char**)action->carg4)[1])) { fprintf(stderr, "Error in ptraj(), cluster comparing: must specify two files for comparing.\n"); } int i; PtrajClustering** Clusterings; #define MAX_COMPARE_FILE 10 Clusterings = (PtrajClustering **) SafeMalloc(__FILE__, __LINE__, sizeof(PtrajClustering *) * (MAX_COMPARE_FILE+1)); for (i=0; icarg4)[i];i++) { Clusterings[i] = PtrajClusteringNew(NULL,trajInfo, action); strcpy(FilePath,((char**)action->carg4)[i]); if (strcmp(FilePath+strlen(FilePath)-4, ".txt")) strcat(FilePath,".txt"); Clusterings[i]->Head = Clusterings[i]->Tail = NULL; Clusterings[i]->PairwiseDistances = pClustering->PairwiseDistances; Clusterings[i]->action = pClustering->action; Clusterings[i]->trajInfo = pClustering->trajInfo; Clusterings[i]->TotalCluster = pClustering->TotalCluster; PtrajClusteringReadFromDisk(Clusterings[i],FilePath); ClusteringFindAllCentroids( (Clustering *) Clusterings[i]); } action->carg6 = (void *) Clusterings; /* PtrajClustering* Clustering2; Clustering2 = PtrajClusteringNew(NULL,pClustering->trajInfo, pClustering->action); strcpy(FilePath,(char*)action->carg6); if (strcmp(FilePath+strlen(FilePath)-4, ".txt")) strcat(FilePath,".txt"); Clustering2->Head = Clustering2->Tail = NULL; Clustering2->PairwiseDistances = pClustering->PairwiseDistances; PtrajClusteringReadFromDisk(Clustering2,FilePath); ClusteringFindAllCentroids(Clustering2); */ FILE* CompareFile; CompareFile = fopen((char *) pClustering->action->carg1, "w"); for (i = 0; ((char**)action->carg4)[i];i++) { fprintf(CompareFile, "Comparefile %2d: %s\n", i+1, ((char**)action->carg4)[i]); } fclose(CompareFile); ClusteringCompare( (Clustering **) Clusterings, (char *) pClustering->action->carg1); break; } PtrajClusteringPrintTransformMap(pClustering); /*printClusters(pClustering);*/ /* FreeMatrix((Matrix*)pClustering->PairwiseDistances);*/ } /** ACTION ROUTINE ************************************************************* * * transformCluster() --- cluster the MD trajectories * * Supplementary routines: **** many!!! *** * * ******************************************************************************/ /* * transformCluster organizes all the frames into clusters. To use this action, * the user should choose a clustering algorithm, and either a number of * clusters or a value of epsilon governing maxiumum cluster size. * * The available clustering algorithms are: * * Hierarchical: Hierarchical divisive clustering first places all points in * one cluster. The algorithm splits, in each iteration, a cluster into two * subclusters; this continues until either the target number of clusters is * reached, or until the maximum eccentricity of any cluster is less than epsilon. * At each iteration, choose the cluster with the largest eccentricity. * Choose two extrema (points with maximal pairwise distance) from this cluster; * these become the centroids of two new clusters, and points are assigned to one child * cluster or the other based on which centroid they are closest to. The * centroid of each child cluster is then calculated, and points * reassigned to child clusters based on these new centroids. Three * centroids are chosen in all. * * Single Linkage: Begin with the original collection of points as your cluster * collection. Choose the two closest points, and form a cluster; remove those two points * from the cluster collection, and insert the centroid of this new cluster. Repeat the * process until all points lie in a cluster. To cluster two centroids, merge their two * clusters into a supercluster with new centroid; to cluster a point and a centroid, add * the point to the centroid's cluster. If an epsilon is provided: Do not cluster two * points if their distance is epsilon or greater. If a cluster count is provided: Continue * until every point is clustered; then, continue merging clusters as long as * the number of clusters exceeds the desired cluster count. * * k-means: If a cluster count n is specified, assign the first n frames to n initial * clusters. Iterate over all the points. Assign each new point to the cluster whose * centroid it is nearest to, and recompute the cluster's centroid. If an epsilon value * is specified: Start with 1 cluster. If a point has distance less than epsilon from a * centroid, add it to an existing cluster; otherwise, form a new cluster with the point. * * Centripetal: Based on the CURE clustering algorithm. * * Cobweb: Based on the COBWEB clustering algorithm. See cluster.h for details. * * Bayesian: Based on the AutoClass clustering algorithm. See cluster.h for details. * * SOM: Self-organizing maps. See cluster.h for details * ******************************************************************************/ #define INVALID_ARGUMENTS_RETURN_CODE -1 extern int transformCluster(actionInformation* action,double* x, double* y, double* z, double* box, int mode) { char *name = "cluster"; stackType** argumentStackPointer; char* buffer; ptrajState* state; trajectoryInfo* trajInfo; trajectoryInfo* FullTrajInfo; SymmetricMatrix* PairwiseDistance; PtrajClustering* TheClustering; int AtomCount; /* USAGE: cluster out [average ] [representative ] [all ] [clusters ] [epsilon ] [sieve [random | start n] ] [verbose] mask Valid clustering algorithms: hierarchical | linkage | means | centripetal | cobweb | bayesian | SOM | edge Valid distance metrics: rms | dme The fileformat argument can be pdb, rest, binpos, none, or amber. Argument usage: iarg1: clustering algorithm iarg2: cluster count iarg3: sieve number. We use this to handle large volumes of data (too many coordinates to cluster all at once because of time and/or memory considerations). If sieve is set to n>0, then we make two passes through the trajectory files. On the first pass, we cluster every nth frame. On the second pass, we add all the remaining frames to whichever of the existing clusters is the best fit. This technique is not guaranteed to work well, especially if the sieve number is large relative to the speed with which the molecule thrashes about. iarg4: frame count - equal to the actual number of frames, unless we're sieving; if sieving, it's equal to the number of sieved frames. iarg5: Distance metric, by number iarg6: mass iarg7: SOM Map darg1: epsilon darg2: Set to non-zero if we're working in verbose mode darg3: sieve mode. less than 0: randomly picking; 0.0|1.0|2.0... the starting frame number darg4: cluster-to-cluster distance. 0.0:best representative; 1.0:centroid (default) darg5: acuity carg1: output filename carg2: trajectory info carg3: trajectory info with unmasked coordinates carg4: file format flags for average, representative, and all carg5: PtrajClustering object carg6: Used in CLUSTER_COMPARE, to store the clusterings. carg7: filename for decoy */ #define C2C_CENTROID 1 #define C2C_BESTREP 0 if (mode == PTRAJ_SETUP) { #ifdef MPI printParallelError(name); return -1; #endif argumentStackPointer = (stackType **) action->carg1; action->carg1 = NULL; /* Parse the output file name */ buffer = argumentStackKeyToString(argumentStackPointer, "out", NULL); if (buffer == NULL) { fprintf(stdout, "WARNING in ptraj(), cluster: you need to specify an output file,\n"); fprintf(stdout, "i.e. \"out cluster\" ...ignoring command.\n"); return INVALID_ARGUMENTS_RETURN_CODE; } action->carg1 = (void *) buffer; /* Parse the clustering algorithm: */ if (argumentStackContains(argumentStackPointer, "hierarchical")) { action->iarg1 = CLUSTER_HIERARCHICAL; } /* CentripetalComplete should be scanned before Centripetal and Complete. */ else if (argumentStackContains(argumentStackPointer, "centripetalcomplete")) { action->iarg1 = CLUSTER_CENTRIPETAL_COMPLETE; } else if (argumentStackContains(argumentStackPointer, "averagelinkage")) { action->iarg1 = CLUSTER_AVERAGELINK; } else if (argumentStackContains(argumentStackPointer, "linkage")) { action->iarg1 = CLUSTER_LINKAGE; } else if (argumentStackContains(argumentStackPointer, "edge")) { action->iarg1 = CLUSTER_EDGELINK; } else if (argumentStackContains(argumentStackPointer, "complete")) { action->iarg1 = CLUSTER_COMPLETELINK; } else if (argumentStackContains(argumentStackPointer, "means1")) { action->iarg1 = CLUSTER_MEANS; } else if (argumentStackContains(argumentStackPointer, "means2")) { action->iarg1 = CLUSTER_MEANS2; } else if (argumentStackContains(argumentStackPointer, "means3")) { action->iarg1 = CLUSTER_MEANS3; } else if (argumentStackContains(argumentStackPointer, "means")) { action->iarg1 = CLUSTER_MEANS; } else if (argumentStackContains(argumentStackPointer, "centripetal")) { action->iarg1 = CLUSTER_CENTRIPETAL; } else if (argumentStackContains(argumentStackPointer, "cobweb1")) { action->iarg1 = CLUSTER_COBWEB; } else if (argumentStackContains(argumentStackPointer, "cobweb2")) { action->iarg1 = CLUSTER_COBWEB2; } else if (argumentStackContains(argumentStackPointer, "cobweb")) { action->iarg1 = CLUSTER_COBWEB; } else if (argumentStackContains(argumentStackPointer, "bayesian")) { action->iarg1 = CLUSTER_BAYESIAN; } else if (argumentStackContains(argumentStackPointer, "som2")) { action->iarg1 = CLUSTER_SOM2; } else if (argumentStackContains(argumentStackPointer, "som")) { action->iarg1 = CLUSTER_SOM; } else if (action->carg7 = (void *) argumentStackKeyToString(argumentStackPointer, "decoy", NULL)) { /* Search for Means, SOM first. */ action->iarg1 = CLUSTER_DECOY; } else if (argumentStackContains(argumentStackPointer, "bychance")) { action->iarg1 = CLUSTER_BYCHANCE; /* Generate random clusters. */ } else if (argumentStackContains(argumentStackPointer, "dbi")) { action->iarg1 = CLUSTER_DBI; /* NO clustering, just read old output and do DBI! */ } else if (argumentStackContains(argumentStackPointer, "anova")) { action->iarg1 = CLUSTER_ANOVA; /* NO clustering, just read old output and do ANOVA! */ } else if (argumentStackContains(argumentStackPointer, "readmerge")) { action->iarg1 = CLUSTER_READMERGE; /* Generate clustering from ClusterMerging.txt */ } else if (argumentStackContains(argumentStackPointer, "readtxt")) { action->iarg1 = CLUSTER_READTXT; /* Generate clustering from previous output .txt */ } else if (argumentStackContains(argumentStackPointer, "clusteringcompare")) { action->iarg1 = CLUSTER_COMPARE; /* Compare clustering from old output files. */ } else { fprintf(stdout,"WARNING in ptraj(), cluster: Unknown clustering algorithm!\n"); return INVALID_ARGUMENTS_RETURN_CODE; } if (!action->carg7) action->carg7 = (void *) argumentStackKeyToString(argumentStackPointer, "decoy", NULL); #define MAX_COMPARE_FILE 10 if (action->iarg1 == CLUSTER_COMPARE) { /* carg4 is borrowed to store filename of clustering to be compared. */ action->carg4 = (char**)SafeMalloc(__FILE__, __LINE__, sizeof(char *) * (MAX_COMPARE_FILE+1)); int i; for (i=0; icarg4)[i] = NULL; char temp[100]; for (i=0; icarg4)[i] = buffer; } sprintf(temp, "comparefile%i", i); buffer = argumentStackKeyToString(argumentStackPointer, temp, NULL); if (buffer != NULL) { fprintf(stdout, "WARNING in ptraj(), cluster comparing: too much files for comparing, ignoring %s and beyond.\n", temp); fprintf(stdout, "i.e. \"comparefile1 clustering1\" ...ignoring command.\n"); return INVALID_ARGUMENTS_RETURN_CODE; } } /* Parse the cluster output formats. The default is to output all clusters as an amber trajectory. */ action->carg4 = (int*)SafeMalloc(__FILE__, __LINE__, sizeof(int) * CLUSTER_FILEFORMAT_LEN); memset(action->carg4,CLUSTER_OUTPUT_NONE,sizeof(int)*CLUSTER_FILEFORMAT_LEN); ((int*)action->carg4)[CLUSTER_FILEFORMAT_ALL] = COORD_AMBER_TRAJECTORY; /* default */ buffer = argumentStackKeyToString(argumentStackPointer, "average", NULL); if (buffer) { ( (int *) action->carg4 )[CLUSTER_FILEFORMAT_AVERAGE] = GetOutputFileFormat(buffer); safe_free(buffer); } buffer = argumentStackKeyToString(argumentStackPointer, "all", NULL); if (buffer) { ((int*)action->carg4)[CLUSTER_FILEFORMAT_ALL] = GetOutputFileFormat(buffer); safe_free(buffer); } buffer = argumentStackKeyToString(argumentStackPointer, "representative", NULL); if (buffer) { ((int*)action->carg4)[CLUSTER_FILEFORMAT_REPRESENTATIVE] = GetOutputFileFormat(buffer); safe_free(buffer); } /* if (argumentStringContains(argumentStackPointer, "average")) { ((int*)action->carg4)[CLUSTER_FILEFORMAT_AVERAGE] = GetOutputFileFormat(argumentStackPointer); } if (argumentStringContains(argumentStackPointer, "representative")) { ((int*)action->carg4)[CLUSTER_FILEFORMAT_REPRESENTATIVE] = GetOutputFileFormat(argumentStackPointer); } if (argumentStringContains(argumentStackPointer, "all")) { ((int*)action->carg4)[CLUSTER_FILEFORMAT_ALL] = GetOutputFileFormat(argumentStackPointer); } */ /* Parse the distance metric */ if (argumentStackContains(argumentStackPointer, "rms")) { action->iarg5 = DISTANCE_METRIC_RMSD; /* default */ } else if (argumentStackContains(argumentStackPointer, "dme")) { action->iarg5 = DISTANCE_METRIC_DME; } else if (argumentStackContains(argumentStackPointer, "mds")) { action->iarg5 = DISTANCE_METRIC_MDS; } /**/ if (argumentStackContains(argumentStackPointer, "fingerprint")) { switch (action->iarg1) { case CLUSTER_AVERAGELINK: action->iarg1 = CLUSTER_AVERAGELINK_FP; break; case CLUSTER_EDGELINK: action->iarg1 = CLUSTER_EDGELINK_FP; break; case CLUSTER_COMPLETELINK: action->iarg1 = CLUSTER_COMPLETELINK_FP; break; case CLUSTER_CENTRIPETAL: action->iarg1 = CLUSTER_CENTRIPETAL_FP; break; case CLUSTER_CENTRIPETAL_COMPLETE: action->iarg1 = CLUSTER_CC_FP; break; case CLUSTER_LINKAGE: action->iarg1 = CLUSTER_LINKAGE_FP; break; } } /* Parse the cluster to cluster measuring method, if it is use the BestRep to BestRep or Centroid to Centroid (default). Actual value is saved in action->darg4 */ char *temp1; temp1 = argumentStackKeyToString(argumentStackPointer, "c2c", "Centroid"); /* default */ if (stringMatch(temp1, "bestrep")) { action->darg4 = C2C_BESTREP; } else if (stringMatch(temp1, "centroid")) { action->darg4 = C2C_CENTROID; } else { fprintf(stdout, "WARNING in ptraj(), cluster: Unknown method for cluster to cluster.\n"); } #define SOM_ISOLATE 0 #define SOM_LOOP 1 #define SOM_STRING 2 #define SOM_BARBLOOP 3 #define SOM_BARBSTRING 4 /* If SOM, The iarg7 will be used for the SOM map. */ if (action->iarg1 == CLUSTER_SOM || action->iarg1 == CLUSTER_SOM2) { if (argumentStackContains(argumentStackPointer, "barbstring")) { action->iarg7 = SOM_BARBSTRING; } else if (argumentStackContains(argumentStackPointer, "barbloop")) { action->iarg7 = SOM_BARBLOOP; } else if (argumentStackContains(argumentStackPointer, "string")) { action->iarg7 = SOM_STRING; } else if (argumentStackContains(argumentStackPointer, "loop")) { action->iarg7 = SOM_LOOP; } else if (argumentStackContains(argumentStackPointer, "isolate")) { action->iarg7 = SOM_ISOLATE; } else { fprintf(stdout,"WARNING in ptraj(), cluster: Unknown SOM map! Use SOM_LOOP map.\n"); action->iarg7 = SOM_LOOP; } } /* If Means or Decoy, The iarg7 will be used for the iteration steps. */ #define ASSIGN_ITERATION_COUNT 100 if (action->iarg1 == CLUSTER_MEANS || action->iarg1 == CLUSTER_MEANS2 || action->iarg1 == CLUSTER_MEANS3 || action->iarg1 == CLUSTER_DECOY) { action->iarg7 = argumentStackKeyToInteger(argumentStackPointer, "iteration", ASSIGN_ITERATION_COUNT); } /* If Cobweb, darg5 is acuity.*/ if (action->iarg1 == CLUSTER_COBWEB || action->iarg1 == CLUSTER_COBWEB2) action->darg5 = argumentStackKeyToDouble(argumentStackPointer, "acuity", 0.1); action->darg2 = argumentStackKeyToDouble(argumentStackPointer, "verbose", 1.0); action->iarg2 = argumentStackKeyToInteger(argumentStackPointer, "clusters", 0); action->iarg3 = argumentStackKeyToInteger(argumentStackPointer, "sieve", 0); /* If we are to perform sieve-clustering, then we demand a second pass through the data! */ if (action->iarg3) { action->performSecondPass = 1; if (argumentStackContains(argumentStackPointer, "random")) { action->darg3 = -1; } else { int temp = argumentStackKeyToInteger(argumentStackPointer, "start", 1); action->darg3 = (double) temp - 1; if (action->darg3 < 0) { fprintf(stdout, "WARNING in ptraj(), cluster: start need to be greater than 0. Now start from frame 1\n"); action->darg3 = 0; } if (action->darg3 > action->iarg3 - 1) { int temp = (int)action->darg3; temp = temp % action->iarg3; action->darg3 = (double) temp; fprintf(stdout, "WARNING in ptraj(), cluster: start need to be less than sieve size. Now start from frame %i\n", temp + 1); } } } if (action->iarg3 == 1) { action->iarg3 = 0; action->performSecondPass = 0; fprintf(stdout, "WARNING in ptraj(), sieve 1 is equivalent to sieve 0 (no sieving).\n"); } if (argumentStackContains(argumentStackPointer, "random")) { fprintf(stdout, "WARNING in ptraj(), cluster. Ignore \"random\" in the cluster command without sieve context.\n"); } if (argumentStackContains(argumentStackPointer, "start")) { fprintf(stdout, "WARNING in ptraj(), cluster. Ignore \"start\" in the cluster command without sieve context, may cause other misintepretations.\n"); } action->darg1 = argumentStackKeyToDouble(argumentStackPointer, "epsilon", 0.0); action->iarg6 = argumentStackContains(argumentStackPointer, "mass"); if (action->iarg2 == 0 && action->darg1 == 0) { fprintf(stdout,"WARNING in ptraj(), cluster: Must specify either 'clusters' or 'epsilon'.\n"); return INVALID_ARGUMENTS_RETURN_CODE; } /* Not all algorithms can handle 'epsilon'. Reject arguments we can't handle */ if (action->darg1 != 0 && \ (action->iarg1 == CLUSTER_COBWEB || action->iarg1 == CLUSTER_BAYESIAN || action->iarg1 == CLUSTER_SOM || action->iarg1 == CLUSTER_SOM2 || action->iarg1 == CLUSTER_MEANS)) { fprintf(stdout,"WARNING in ptraj(), cluster: This clustering algorhtim can't handle \ 'epsilon'. Please specify a number for clusters.\n"); return INVALID_ARGUMENTS_RETURN_CODE; } /* For Edge and Complete linkage, the distance between cluster and cluster is not defined by it centroid of bestrep. So, if Edge or Complete linkage is selected, use bestrep to save alignment time.*/ /* if (action->iarg1 == CLUSTER_EDGELINK || action->iarg1 == CLUSTER_COMPLETELINK) { if (action->darg4 != C2C_BESTREP) { fprintf(stdout,"WARNING in ptraj(), cluster: use Bestrep for edge linkage or complete linkage.\n"); action->darg4 = C2C_BESTREP; } } */ if (action->iarg1 == CLUSTER_CENTRIPETAL) { if (action->darg4 == C2C_BESTREP) { fprintf(stdout,"WARNING in ptraj(), cluster: use Centroid for centripetal linkage.\n"); action->darg4 = C2C_CENTROID; } } /* Most DBI calculation is need by the sieve treatment. So disabling the sieve will accumulate the coordinates of all the frames and get DBI and pSF calculated. */ if (action->iarg1 == CLUSTER_DBI) { action->iarg3 = 0; action->performSecondPass = 0; } AtomCount = action->state->atoms; trajInfo = (trajectoryInfo *) SafeMalloc(__FILE__, __LINE__, sizeof(trajectoryInfo)); INITIALIZE_trajectoryInfo(trajInfo); buffer = getArgumentString(argumentStackPointer, NULL); fprintf(stdout, "MASK = %s\n", buffer); if (buffer != NULL) { action->mask = processAtomMask(buffer, action->state); safe_free(buffer); } else { action->mask = processAtomMask("*", action->state); } if (action->mask==NULL) return INVALID_ARGUMENTS_RETURN_CODE; /* Create a new "substate" trajInfo->state containing only some atoms: */ modifyStateByMask(&trajInfo->state, &action->state, action->mask, 0); trajInfo->atoms = trajInfo->state->atoms; action->carg2 = (void *) trajInfo; /* We want to save *all* the coordinates, so that we can output complete cluster information. So, if the mask excludes some atoms, we use a second trajectoryInfo structure to hold the positions of the full, unmasked trajectory. */ action->carg3 = NULL; /* default */ if (AtomCount != trajInfo->atoms) { FullTrajInfo = (trajectoryInfo *) SafeMalloc(__FILE__, __LINE__, sizeof(trajectoryInfo)); INITIALIZE_trajectoryInfo(FullTrajInfo); FullTrajInfo->atoms = AtomCount; FullTrajInfo->state = action->state; action->carg3 = (void *) FullTrajInfo; } TheClustering = PtrajClusteringNew(NULL,trajInfo,action); TheClustering->Acuity = (float)action->darg5; action->carg5 = TheClustering; return 0; } if (mode == PTRAJ_STATUS) { buffer = (char*)action->carg1; if (action->iarg1 == CLUSTER_COMPARE) { fprintf(stdout, " Comparing clustering files, output file: %s\n", buffer); int i; for (i = 0; ((char**)action->carg4)[i]; i++) { fprintf(stdout, " comparefile%d: %s\n", i+1,((char**)action->carg4)[i]); } return 0; } fprintf(stdout, " CLUSTERING.\n\n"); fprintf(stdout, " Clusters will be assigned and output to files with the prefix \"%s\".\n", buffer); fprintf(stdout, " The clustering algorithm is "); if (action->iarg1 != CLUSTER_SOM && action->iarg1 != CLUSTER_SOM2 && action->iarg1 != CLUSTER_COBWEB && action->iarg1 != CLUSTER_BAYESIAN) { fprintf(stdout,"%s using %s as the distance metric.\n", CLUSTER_ALGORITHM_NAMES[action->iarg1], DISTANCE_METRIC_NAMES[action->iarg5]); if (action->darg4 == C2C_BESTREP) fprintf(stdout," Cluster-to-cluster distance is measured by its best representative.\n"); if (action->darg4 == C2C_CENTROID) fprintf(stdout," Cluster-to-cluster distance is measured by its centroid.\n"); } else { fprintf(stdout,"%s, ", CLUSTER_ALGORITHM_NAMES[action->iarg1]); if (attributeArray != NULL || attributeArrayTorsion != NULL) { fprintf(stdout," using attributes defined by previous ClusterAttribute commands.\n"); } else fprintf(stdout," using XYZ coordinates as attributes.\n"); } /* * If SOM, print the SOM map. */ if (action->iarg1 == CLUSTER_SOM || action->iarg1 == CLUSTER_SOM2) { if (action->iarg7 == SOM_BARBSTRING) { fprintf(stdout," SOM map is Barbed string.\n"); } else if (action->iarg7 == SOM_BARBLOOP) { fprintf(stdout," SOM map is Barbed loop.\n"); } else if (action->iarg7 == SOM_STRING) { fprintf(stdout," SOM map is String (open loop).\n"); } else if (action->iarg7 == SOM_LOOP) { fprintf(stdout," SOM map is Loop.\n"); } else if (action->iarg7 == SOM_ISOLATE) { fprintf(stdout," SOM map is Isolated points.\n"); } } /* * If Means, print the iteration. */ if (action->iarg1 == CLUSTER_MEANS) { fprintf(stdout," Maximum of %d iterations with means.\n", action->iarg7); } if (action->carg7) { fprintf(stdout," Decoy structure will be read from %s.\n", action->carg7); } if (action->iarg2) { fprintf(stdout," The cluster count is set to %d.\n",action->iarg2); } else { fprintf(stdout," Maximum eccentricity (epsilon) is %f.\n",action->darg1); } if (action->iarg3) { fprintf(stdout," Clustering will be performed in two passes (sieve method).\n"); fprintf(stdout," Sieved frames in the first pass will be selected "); if (action->darg3 < 0) fprintf(stdout, "randomly"); if (action->darg3>= 0) fprintf(stdout, "every %i frames starting from frame %i", action->iarg3, (int)action->darg3 +1); fprintf(stdout, ".\n"); } fprintf(stdout, " The atom selection for best fit is "); printAtomMask(stdout, action->mask, action->state); fprintf(stdout, "\n"); PrintClusterFileInfo(((int*)action->carg4)[CLUSTER_FILEFORMAT_ALL],"all"); PrintClusterFileInfo(((int*)action->carg4)[CLUSTER_FILEFORMAT_REPRESENTATIVE],"representative"); PrintClusterFileInfo(((int*)action->carg4)[CLUSTER_FILEFORMAT_AVERAGE],"average"); } if (mode == PTRAJ_PRINT) { if (action->iarg1 == CLUSTER_DBI) return 0; if (action->iarg1 == CLUSTER_ANOVA) return 0; if (action->iarg1 == CLUSTER_COMPARE) return 0; fprintf(stdout, "\nPTRAJ CLUSTER:"); PtrajPrintClustering(action); } /* end of PTRAJ_PRINT */ if (mode == PTRAJ_CLEANUP) { trajInfo = (trajectoryInfo *) action->carg2; safe_free(trajInfo->x); safe_free(trajInfo->y); safe_free(trajInfo->z); trajInfo->x = NULL; trajInfo->y = NULL; trajInfo->z = NULL; safe_free(trajInfo); if (action->carg3) { trajInfo = (trajectoryInfo *) action->carg3; safe_free(trajInfo->x); safe_free(trajInfo->y); safe_free(trajInfo->z); trajInfo->x = NULL; trajInfo->y = NULL; trajInfo->z = NULL; safe_free(trajInfo); } safe_free(action->carg4); safe_free(action->carg7); } if (mode == PTRAJ_FIRSTPASS) { /*This is where the meat of the computation occurs; up until now, we were just accumulating the necessary coordinates. */ fprintf(stdout,"\n\n The first pass through the trajectory is complete.\n"); trajInfo = (trajectoryInfo *) action->carg2; PtrajPerformClustering(trajInfo,action); if (action->iarg1 == CLUSTER_DBI || action->iarg1 == CLUSTER_ANOVA || action->iarg1 == CLUSTER_COMPARE) return 0; AlignBestReps((PtrajClustering*)action->carg5); /* Get ready for the second pass, if we are sieving: */ if (action->iarg3) { /* Our centroids will not necessarily be up-to-date, especially if we're doing an algorithm like Cobweb Clustering that doesn't use them. We need centroids for the second pass, so update them! */ fprintf(stdout," Sieve was active. Preparing for the second pass.\n"); fprintf(stdout," Dumping out current clusters to file EndFirstPass.txt.\n"); FILE* temp; ((PtrajClustering*)action->carg5)->FirstPassTime = time(NULL); temp = fopen("EndFirstPass.txt", "w"); int tempclusters = ((PtrajClustering*)action->carg5)->action->iarg3; /* * to fool OutputClusteringStats to print out the DBI and pSF value for the first pass */ ((PtrajClustering*)action->carg5)->action->iarg3 = 0; PtrajClusteringOutputHeaderToFile((PtrajClustering*)action->carg5,temp); OutputClusteringStats((PtrajClustering*)action->carg5,temp); ClusteringOutputClusterListToFile( (Clustering*) action->carg5,temp,0); fclose(temp); ((PtrajClustering*)action->carg5)->action->iarg3 = tempclusters; /*ClusteringOutputClusterList((PtrajClustering*)action->carg5,"EndFirstPass.txt");*/ fprintf(stdout," Computing centroids"); fflush(stdout); ClusteringFindAllCentroids( (Clustering *) action->carg5); fprintf(stdout," and expanding clusters.\n"); PtrajClusteringExpandClusterListSizes((PtrajClustering*)action->carg5, action->iarg3,action->state->maxFrames); trajInfo->current = 0; PtrajClusteringSetupSecondPassInfo((PtrajClustering*)action->carg5,action); fprintf(stdout," Ready for the second pass through the trajectories...\n"); } } if (mode == PTRAJ_SECONDPASS) { if (action->iarg3 && action->iarg1 != CLUSTER_DBI) { PtrajClusteringSecondPass(action,x,y,z); } } if (mode != PTRAJ_ACTION) { /* We're done: The remainder of the routine handles the case of action */ return 0; } /* ACTION: PTRAJ_ACTION Load up the current coordinates into the trajInfo structure. We keep accumulating until the PRINT event triggers. */ state = (ptrajState *) action->state; trajInfo = (trajectoryInfo *) action->carg2; FullTrajInfo = (trajectoryInfo *) action->carg3; AccumulateCoordinates(trajInfo, action, x, y, z, FullTrajInfo); return 0; } /** ACTION ROUTINE **** transformClusterAttribute Argument usage: iarg1: attribute from stack|file|xyz|backbone iarg2: data manupilation. Nochange/Normalize/Meanshift/Autoscale/Weighted iarg3: iarg4: iarg5: darg1: output time, default 1.0 ps darg2: weight, valid only if defined "weight", default 1.0, darg3: darg4: carg1: temporary for command parameter stack carg2: filename and colume info for attributes input carg3: filename for output carg4: a stackType pointer to the attributes list defined by action carg5: */ #define ATTRIBUTE_STACK 1 #define ATTRIBUTE_FILE 2 #define ATTRIBUTE_BACKBONE 3 #define ATTRIBUTE_XYZ 0 #define ATTRIBUTE_NOCHANGE 0 #define ATTRIBUTE_NORMALIZE 1 #define ATTRIBUTE_MEANSHIFT 2 #define ATTRIBUTE_AUTOSCALE 3 #define ATTRIBUTE_LINEARSCALE 4 #define ATTRIBUTE_WEIGHT 5 #define ATTRIBUTE_EIGHTY_PERCENTILE 8 extern int transformClusterAttribute(actionInformation* action, double* x, double* y, double* z, double* box, int mode) { char *name = "clusterattribute"; stackType** argumentStackPointer; stackType *s, *ss; scalarInfo *info, *match; scalarInfo *attributeInfo; scalarInfo **array; char **attributes, *a, *filename; int *columes; char* buffer; char *line, *col; double value; int invalid; int i, j, k, start, end; double **array1; double **array2; FILE *outFile; FILE *inFile; int torsion; /* define the data is torsion/dihe or not when read data from a file */ int angle; int SkipAttribute = 0; double max, min; double sum, sum2, stdev, mean; double sinX, cosX, sinX2, cosX2; if (mode == PTRAJ_SETUP) { #ifdef MPI printParallelError(name); return -1; #endif argumentStackPointer = (stackType **) action->carg1; action->carg1 = NULL; /* * grab a time interval between frames in ps (for output) */ action->darg1 = argumentStackKeyToDouble(argumentStackPointer, "time", 1.0); /* * grab output filename */ buffer = argumentStackKeyToString(argumentStackPointer, "out", NULL); action->iarg2 = ATTRIBUTE_NOCHANGE; if (argumentStackContains(argumentStackPointer, "normalize")) { action->iarg2 = ATTRIBUTE_NORMALIZE; } else if (argumentStackContains(argumentStackPointer, "meanshift")) { action->iarg2 = ATTRIBUTE_MEANSHIFT; } else if (argumentStackContains(argumentStackPointer, "autoscale")) { action->iarg2 = ATTRIBUTE_AUTOSCALE; } else if (argumentStackContains(argumentStackPointer, "linearscale")) { action->iarg2 = ATTRIBUTE_LINEARSCALE; } else if (argumentStackContains(argumentStackPointer, "percentile")) { action->iarg2 = ATTRIBUTE_EIGHTY_PERCENTILE; } action->darg2 = argumentStackKeyToDouble(argumentStackPointer, "weight", 1.0); if (action->darg2 == 0) fprintf(stdout, "WARNING in ptraj(), clusterAttribute: weight is 0. This attribute will have no effect.\n"); action->carg3 = (void *)buffer; if ( argumentStackContains(argumentStackPointer, "pairwisedistance") ) { action->iarg1 = ATTRIBUTE_STACK; buffer = getArgumentString(argumentStackPointer, NULL); if (buffer != NULL) { action->mask = processAtomMask(buffer, action->state); safe_free(buffer); } else { action->mask = processAtomMask("*", action->state); } char TempCommand[512]; stackType *argumentStack = NULL; for (i = 1; i <= action->state->atoms; i++) { if (action->mask[i-1] == 0) continue; for (j = i+1; j <= action->state->atoms; j++) { if (action->mask[j-1] == 0) continue; sprintf(TempCommand, "distance tempdist_%d_%d @%d @%d\0", i, j, i, j); printf("%s\n", TempCommand); dispatchToken((Token *) &ptrajTokenlist, argumentStack, TempCommand); } } /* Modify the current ClusterAttribute command by pushing "tempdist*" and "stack" (reversed order) into the argumentStack. */ buffer = safe_malloc(sizeof(char)*512); pushStack(argumentStackPointer, strcpy(buffer, "tempdist*")); buffer = safe_malloc(sizeof(char)*512); pushStack(argumentStackPointer, strcpy(buffer, "stack")); /*sprintf(TempCommand, "ClusterAttribute stack tempdist* \0"); dispatchToken((Token *) &ptrajTokenlist, argumentStack, TempCommand);*/ } if (buffer = argumentStackKeyToString(argumentStackPointer, "stack", NULL)) { action->iarg1 = ATTRIBUTE_STACK; action->carg2 = buffer; attributes = stringSplit(buffer, ", "); for (i = 0; attributes[i] != NULL; i++) { for (s = scalarStack; s != NULL; s = s->next) { info = (scalarInfo *) s->entry; attributeInfo = NULL; if ( stringMatchWild(info->name, attributes[i]) == 1 ) { attributeInfo = (scalarInfo *) SafeMalloc(__FILE__, __LINE__, sizeof(scalarInfo)); INITIALIZE_scalarInfo(attributeInfo); attributeInfo->mode = info->mode; attributeInfo->totalFrames = -1; attributeInfo->name = info->name; attributeInfo->results = info; /* If buffer already in the attributeStack, skip*/ match = scalarStackGetName(&attributeStack, attributeInfo->name); if (match != NULL) { safe_free(attributeInfo); attributeInfo = NULL; } if (attributeInfo != NULL) { pushStack(&attributeStack, (void*) attributeInfo); pushStack( (stackType **) &(action->carg4), (void*) attributeInfo); } } } } } else if (filename = argumentStackKeyToString(argumentStackPointer, "file", NULL)) { action->iarg1 = ATTRIBUTE_FILE; buffer = argumentStackKeyToString(argumentStackPointer, "column", NULL); action->carg2 = SafeMalloc(__FILE__, __LINE__, sizeof(char) * (strlen(filename)+ strlen(buffer) + 10)); sprintf(action->carg2, "%s", filename); attributes = stringSplit(buffer,", "); torsion = argumentStackContains(argumentStackPointer, "torsion"); if (!torsion) torsion = argumentStackContains(argumentStackPointer, "dihedral"); angle = argumentStackContains(argumentStackPointer, "angle"); for (i = 0; attributes[i] != NULL; i++) { attributeInfo = NULL; if (strchr(attributes[i], '-')) { if (sscanf(attributes[i], "%i-%i", &start, &end) < 1) { fprintf(stderr, "Scanning integer argument failed in %s\n", attributes[i]); } } else { if (sscanf(attributes[i], "%i", &start) < 1) { fprintf(stderr, "Scanning integer argument failed in %s\n", attributes[i]); } else { end = start; } } for (j = start; j <= end; j++) { attributeInfo = (scalarInfo *) SafeMalloc(__FILE__, __LINE__, sizeof(scalarInfo)); INITIALIZE_scalarInfo(attributeInfo); if (torsion) attributeInfo->mode = SCALAR_TORSION; else if (angle) attributeInfo->mode = SCALAR_ANGLE; else attributeInfo->mode = SCALAR_NULL; attributeInfo->totalFrames = -1; /* attributeInfo->name = SafeMalloc(__FILE__, __LINE__, sizeof(char) * (strlen(filename)+ (log10(j+1)) + 1));*/ attributeInfo->name = SafeMalloc(__FILE__, __LINE__, (size_t)(sizeof(char) * (strlen(filename)+ log10(j+1.0) + 10))); sprintf(attributeInfo->name, "%s:%d", filename, j); attributeInfo->results = NULL; /* If buffer already in the attributeStack, skip*/ match = scalarStackGetName(&attributeStack, attributeInfo->name); if (match != NULL) { safe_free(attributeInfo->name); safe_free(attributeInfo); attributeInfo = NULL; } if (attributeInfo != NULL) { pushBottomStack(&attributeStack, (void*) attributeInfo); pushBottomStack( (stackType **) &(action->carg4), (void*) attributeInfo); } } } } else if (argumentStackContains(argumentStackPointer, "backbone")) { action->iarg1 = ATTRIBUTE_BACKBONE; } else { action->iarg1 = ATTRIBUTE_XYZ; } return 0; } else if (mode == PTRAJ_FIRSTPASS) { /* array1 = (double **)attributeArray->entry; array2 = (double **)attributeArrayTorsion->entry; */ array = (scalarInfo**)attributeArray->entry; if (action->iarg1 == ATTRIBUTE_FILE) { /* Read in the file */ s = (stackType *)action->carg4; /* filename = SafeMalloc(__FILE__, __LINE__, sizeof(char)*strlen((char*)s->entry)); */ if (!action->carg3) { fprintf(stdout, "WARNING in ptraj(), No attribute info from %s will be saved.\n", action->carg2); } if (openFile(&inFile, (char*)action->carg2, "r") == 0) { fprintf(stdout, "WARNING in ptraj(), ClusterAttribute: couldn't open file %s.\n", action->carg2); return 0; } for (s = (stackType *)action->carg4, i = 0, j = 0; s != NULL; s = s->next) { info = (scalarInfo *) s->entry; info->totalFrames = action->state->maxFrames; info->value = (double *) SafeMalloc(__FILE__, __LINE__, sizeof(double) * info->totalFrames); } line = (char *) SafeMalloc(__FILE__, __LINE__, sizeof(char) * BUFFER_SIZE); k = 0; while (fgets(line, BUFFER_SIZE, inFile)) { attributes = stringSplit(line,"\t "); invalid = 0; for (s = (stackType *)action->carg4; s != NULL; s = s->next) { info = (scalarInfo *) s->entry; col = (char *) strrchr( info->name, ':'); *col = '\0'; if (strcmp((char*)action->carg2, info->name)) { *col = ':'; continue; } *col = ':'; col++; j = atoi(col) - 1; value = atof(attributes[j]); if (value == 0 && attributes[j][0] != '0') { invalid = 1; continue; } info->value[k] = value; } if (invalid == 0) { k++; } else { fprintf(stdout, "Line %d in %s is not valid. Skip. \n", k+1, (char*)action->carg2); } for (j = 0; attributes[j]; j++) safe_free(attributes[j]); safe_free(attributes); } for (i = 0; array[i]; i++); for (s = (stackType *)action->carg4, j = 0; s != NULL; s = s->next, i++) { info = (scalarInfo *) s->entry; info->results = info; col = (char *) (strrchr(info->name, ':') - info->name); if (strncmp( (char *) action->carg2, (char *) info->name, (size_t) col)) { continue; } /* if (info->mode != SCALAR_TORSION) array1[i++] = (double *)info->value; else array2[j++] = (double *)info->value; */ if (info->mode != SCALAR_TORSION) array[i] = (scalarInfo *)info; else { info->cos = (double *) SafeMalloc(__FILE__, __LINE__, sizeof(double) * info->totalFrames); info->sin = (double *) SafeMalloc(__FILE__, __LINE__, sizeof(double) * info->totalFrames); for (k = 0; k < info->totalFrames; k++) { if (info->value[k] < 0) info->value[k] += 360; info->cos[k] = cos(info->value[k]/RADDEG); info->sin[k] = sin(info->value[k]/RADDEG); } array[i] = (scalarInfo *)info; } } fclose(inFile); } else if (action->iarg1 == ATTRIBUTE_STACK) { for (i = 0; array[i]; i++); for (s = (stackType *)action->carg4, j = 0; s != NULL; s = s->next, i++) { info = (scalarInfo *) s->entry; info->totalFrames = action->state->maxFrames; info->value = (double *) SafeMalloc(__FILE__, __LINE__, sizeof(double) * info->totalFrames); attributeInfo = (scalarInfo *)info->results; info->results = NULL; if (info->action) { scalarInfo *scalar = ptrajCopyScalar(&info); scalar->totalFrames = 1; scalar->mean = scalar->stddev = scalar->max = scalar->min = 0; scalar->value = SafeMalloc(__FILE__, __LINE__, sizeof(double)*scalar->totalFrames);/* just for one value. */ actionInformation *actionInfo = ptrajCopyAction(&(info->action)); if (actionInfo->type == TRANSFORM_ANGLE || actionInfo->type == TRANSFORM_DIHEDRAL || actionInfo->type == TRANSFORM_DISTANCE || actionInfo->type == TRANSFORM_PUCKER) { /* Although the rms also store one scalarInfo for the rms value, but it is saved on carg2. */ actionInfo->carg1 = (void *)scalar; info->action = actionInfo; } else { fprintf(stderr, "Caution: This attribute can not be updated during clustering!\n"); info->action = NULL; } } /* if (info->mode != SCALAR_TORSION) array1[i++] = (double *)attributeInfo->value; else array2[j++] = (double *)attributeInfo->value; */ if (info->mode != SCALAR_TORSION) { for (k = 0; k < info->totalFrames; k++) { info->value[k] = attributeInfo->value[k]; } array[i] = (scalarInfo *)info; } else { info->cos = (double *) SafeMalloc(__FILE__, __LINE__, sizeof(double) * info->totalFrames); info->sin = (double *) SafeMalloc(__FILE__, __LINE__, sizeof(double) * info->totalFrames); for (k = 0; k < info->totalFrames; k++) { if (attributeInfo->value[k] < 0) attributeInfo->value[k] += 360; info->value[k] = attributeInfo->value[k]; info->cos[k] = cos(attributeInfo->value[k]/RADDEG); info->sin[k] = sin(attributeInfo->value[k]/RADDEG); } array[i] = (scalarInfo *)info; } } } else if (action->iarg1 == ATTRIBUTE_BACKBONE) { } else if(action->iarg1 == ATTRIBUTE_XYZ) { /* # Use the old XYZ functions */ } /* Normalization */ if ( 1 /* Always calculate the mean and stddev. */ ||action->iarg2 == ATTRIBUTE_NORMALIZE || action->iarg2 == ATTRIBUTE_MEANSHIFT || action->iarg2 == ATTRIBUTE_AUTOSCALE || action->iarg2 == ATTRIBUTE_LINEARSCALE || action->iarg2 == ATTRIBUTE_EIGHTY_PERCENTILE) { for (s = (stackType *)action->carg4; s != NULL; s = s->next, i++) { info = (scalarInfo *) s->entry; for (i = 0; strcmp(array[i]->name, info->name);i++); /* col = strrchr(info->name, ':') - info->name; if (strncmp((char*)action->carg2, info->name, (int)col)) { continue; } */ max = min = (array[i]->value)[0]; sum = sum2 = 0; sinX = cosX = sinX2 = cosX2 = 0; for (k = 0; k < info->totalFrames; k++) { value = (array[i]->value)[k]; sum += value; sum2 += value * value; max = max(max, value); min = min(min, value); if (info->mode == SCALAR_TORSION) { sinX += (array[i]->sin)[k]; sinX2 += (array[i]->sin)[k] * (array[i]->sin)[k]; cosX += (array[i]->cos)[k]; cosX2 += (array[i]->cos)[k] * (array[i]->cos)[k]; } } mean = sum / info->totalFrames; stdev = sqrt(sum2 / info->totalFrames - mean * mean ); if (info->mode == SCALAR_TORSION) { mean = atan2(sinX, cosX) * RADDEG; if (mean < 0) mean += 360; stdev = (sinX2 + cosX2)/info->totalFrames - (sinX*sinX+cosX*cosX)/(info->totalFrames*info->totalFrames); stdev = acos(1-stdev/2) *RADDEG; } fprintf(stdout, "attribute %4d: %-20s\tmin is %.2f, max is %.2f, mean is %.3f, stdev is %.3f , ", i, array[i]->name, min, max, mean, stdev); info->max = max; info->min = min; info->mean = mean; info->stddev = stdev; float FitNormal(scalarInfo*); float fitnorm = FitNormal(info); fprintf(stdout, "fitnorm is %f\n", fitnorm); if (fitnorm < -1000) { /* Do not omit any attributes now. */ fprintf(stdout, " !!This attribute shapes like normal distribution, clustering will omit this attribute\n"); attributeArray->length--; for (k = i; k < attributeArray->length; k++) { array[k] = array[k+1]; } } if (info->mode != SCALAR_TORSION && info->mode != SCALAR_ANGLE) { if (action->iarg2 == ATTRIBUTE_NORMALIZE) { for (k = 0; k < info->totalFrames; k++) { info->value[k] = (info->value[k] - mean) / stdev; } fprintf(stdout, "Normalize for scalar value %s.\n", info->name); } if (action->iarg2 == ATTRIBUTE_MEANSHIFT) { for (k = 0; k < info->totalFrames; k++) { info->value[k] = info->value[k] - mean; } fprintf(stdout, "Meanshift for scalar value %s.\n", info->name); } if (action->iarg2 == ATTRIBUTE_AUTOSCALE) { for (k = 0; k < info->totalFrames; k++) { info->value[k] = (info->value[k] - mean) / stdev + mean; } fprintf(stdout, "Autoscale for scalar value %s.\n", info->name); } if (action->iarg2 == ATTRIBUTE_LINEARSCALE) { for (k = 0; k < info->totalFrames; k++) { info->value[k] = (info->value[k] - min)/(max - min); } fprintf(stdout, "Linearscale for scalar value %s.\n", info->name); } if (action->iarg2 == ATTRIBUTE_EIGHTY_PERCENTILE) { float *tempval = SafeMalloc(__FILE__, __LINE__, info->totalFrames * sizeof(float)); memset(tempval, 0, info->totalFrames * sizeof(float)); int l; for (k = 0; k < info->totalFrames; k++) { for (l = k-1; l >= 0; l--) { if (tempval[l] > info->value[k]) tempval[l+1] = tempval[l]; else break; } tempval[l+1] = info->value[k]; } int low, up; low = (int)(info->totalFrames / 10); up = (int)(info->totalFrames * 9 / 10); fprintf(stdout, "80%% percentile is (%.2f, %.2f). its ratio to stdev is %.2f for scalar value %s.\n", tempval[low], tempval[up], (tempval[up]-tempval[low])/stdev, info->name); safe_free(tempval); } } else { fprintf(stdout, "WARNING in ptraj() ClusterAttribute, %s is not a scalar value, so no normalization will be done.\n", info->name); continue; } } } } else if (mode == PTRAJ_STATUS) { for (s = attributeStack, i = 0, j = 0; s != NULL; s = s->next) { info = (scalarInfo *) s->entry; if (info->mode != SCALAR_TORSION) i++; else j++; } /* attributeArray = (arrayType *)SafeMalloc(__FILE__, __LINE__, sizeof(arrayType)); attributeArray->length = i; attributeArray->entry = (double **)SafeMalloc(__FILE__, __LINE__, sizeof(double*) * i); array1 = (double **)attributeArray->entry; attributeArrayTorsion = (arrayType *)SafeMalloc(__FILE__, __LINE__, sizeof(arrayType)); attributeArrayTorsion->length = j; attributeArrayTorsion->entry = (double **)SafeMalloc(__FILE__, __LINE__, sizeof(double*) * j); array2 = (double **)attributeArrayTorsion->entry; */ if (!attributeArray) { attributeArray = (arrayType *)SafeMalloc(__FILE__, __LINE__, sizeof(arrayType)); attributeArray->length = i+j; attributeArray->entry = (scalarInfo **)SafeMalloc(__FILE__, __LINE__, sizeof(scalarInfo *) * (i+j)); memset(attributeArray->entry, 0, sizeof(scalarInfo *) * (i+j)); } fprintf(stdout, " Assign cluster attributes "); if (action->darg2 != 1) fprintf(stdout, "with weight of %.2f ", action->darg2); if (action->iarg2 == ATTRIBUTE_NORMALIZE) fprintf(stdout, "as NORMALIZED attributes\n"); if (action->iarg2 == ATTRIBUTE_MEANSHIFT) fprintf(stdout, "as meanshifting adjusted attributes\n"); if (action->iarg2 == ATTRIBUTE_AUTOSCALE) fprintf(stdout, "as autoscaling adjusted attributes\n"); if (action->iarg2 == ATTRIBUTE_LINEARSCALE) fprintf(stdout, "as linear scaling adjusted attributes\n"); if (action->iarg1 == ATTRIBUTE_STACK) { fprintf(stdout, " Cluster attributes come from %s in stack.\n", action->carg2); fprintf(stdout, " Stack selections are "); if (action->carg4 == NULL) { fprintf(stdout, "No found. "); } else { for (s = (stackType *)action->carg4; s != NULL; s = s->next) { fprintf(stdout, "%s,", ((scalarInfo *)s->entry)->name); } } fprintf(stdout, "\b\n"); } else if (action->iarg1 == ATTRIBUTE_FILE) { fprintf(stdout, " Cluster attributes come from in file %s.\n", action->carg2); } else if (action->iarg1 == ATTRIBUTE_BACKBONE) { } else if(action->iarg1 == ATTRIBUTE_XYZ) { /* # Use the old XYZ functions */ } } else if (mode == PTRAJ_PRINT) { if (!action->carg3) { fprintf(stdout, "WARNING in ptraj(), No attribute info will be saved for %s\n", action->carg2); return 0; } if (openFile(&outFile, action->carg3, "w") == 0) { fprintf(stdout, "WARNING in ptraj(), ClusterAttribute: couldn't open file %s\n", action->carg3); return 0; } fprintf(stdout, "PTRAJ ClusterAttribute dumping attributes into %s\n", action->carg3); /*reverseStack((stackType **)&(action->carg4));*/ fprintf(stdout, " attributes "); for (s = (stackType *)attributeStack; s != NULL; s = s->next) { info = (scalarInfo *) s->entry; fprintf(stdout, "%s,", info->name); } fprintf(stdout, "\b\n"); fprintf(outFile, " Time \t"); for (s = (stackType *)attributeStack; s != NULL; s = s->next) { info = (scalarInfo *) s->entry; fprintf(outFile, "%8s\t", info->name); } fprintf(outFile, "\n"); for (i = 0; i < action->state->maxFrames; i++) { fprintf(outFile, "%10.2f\t", (i + 1) * action->darg1); for (s = (stackType *)attributeStack; s != NULL; s = s->next) { info = (scalarInfo *) s->entry; /*info = (scalarInfo *) info->results;*/ fprintf(outFile, "%f\t", info->value[i]); } fprintf(outFile, "\n"); } safe_fclose(outFile); } else if (mode == PTRAJ_CLEANUP) { safe_free(action->carg2); safe_free(action->carg3); s = (stackType *)action->carg4; while (s != NULL ) { ss = s->next; info = (scalarInfo *) s->entry; if (action->iarg1 != ATTRIBUTE_STACK) safe_free(info->name); /* The name of the attributes in stack will be freed by the operations in stack. Not here. */ safe_free(info->filename); safe_free(info->cos); safe_free(info->sin); safe_free(info); safe_free(s); s = ss; } } else if (mode != PTRAJ_ACTION) { /* We're done: The remainder of the routine handles the case of action */ return 0; } /* ACTION: PTRAJ_ACTION Load up the current coordinates into the trajInfo structure. We keep accumulating until the PRINT event triggers. */ return 0; } /* 1/sqrt(2*pi) */ #define GAUSSIAN_MULTIPLIER 0.39894228040 /* This function is to fit the values of an attribute to the normal distribution. Here, we do not do best fit, just fit the normal distribution using the mean and standard deviation of the values. */ float FitNormal(scalarInfo* info) { double max, min, range; double mean, stdev; float* bins; int bin, binCount; double span; int i,j; float totaldev, dev, var; float x, gpd; max = info->max; min = info->min; mean = info->mean; stdev = info->stddev; range = max - min; max += range/2; min -= range/2; binCount = (int)(info->totalFrames / 5); if (binCount > 100) binCount = 100; if (info->mode == SCALAR_TORSION || info->mode == SCALAR_ANGLE) binCount = 360; bins = (float*)SafeMalloc(__FILE__, __LINE__, sizeof(float*) * binCount); memset(bins,0,sizeof(float*)*binCount); span = (max - min) / binCount; for (i = 0; i < info->totalFrames; i++) { bin = (int) (info->value[i] - min) / span; bins[bin] += 1; } for (i = 0; i < binCount; i++) { bins[i] = bins[i] / info->totalFrames * binCount / (max - min) ; } if (info->mode == SCALAR_TORSION || info->mode == SCALAR_ANGLE) { int lowpos = 0; int low = bins[0]; int lowlength = 1; int newlow = bins[0]; int newlowpos = 0; int newlowlength = 1; for (i = 1; i < binCount * 2; i++) { int currpos = (i>=binCount)?i-binCount:i; if (bins[currpos] < low) { lowpos = currpos; low = bins[currpos]; lowlength = 1; } else if (bins[currpos] == low) { if (lowpos == currpos - lowlength) lowlength++; else if (lowpos == currpos - newlowlength) { newlowlength++; } else { newlow = bins[currpos]; newlowpos = currpos; newlowlength = 1; } if (newlowlength > lowlength) { low = newlow; lowpos = newlowpos; lowlength = newlowlength; } } lowpos = lowpos +lowlength/2; lowpos = lowpos > 360 ? lowpos - 360: lowpos; } totaldev = 0; for (i = 0; i < binCount; i++) { if (i * span > lowpos) x = i*span - lowpos; else x = i*span +360; gpd = GAUSSIAN_MULTIPLIER / stdev * exp (-(x-mean)*(x-mean)/(2*stdev*stdev)); var = (gpd - bins[i]) * (gpd - bins[i]); totaldev += var; } dev = sqrt(totaldev * span / binCount); } else { totaldev = 0; for (i = 0; i < binCount; i++) { x = min + i * span + span / 2; gpd = GAUSSIAN_MULTIPLIER / stdev * exp (-(x-mean)*(x-mean)/(2*stdev*stdev)); var = (gpd - bins[i]) * (gpd - bins[i]); totaldev += var; } dev = sqrt(totaldev * span / binCount); } return(dev); } /** 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) { char *name = "closestwaters"; 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] [noimage] * * 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 */ #ifdef MPI #endif 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(stdout, action->mask, oldstate); fprintf(stdout, "\n"); if (action->iarg3) fprintf(stdout, " Imaging of the coordinates will not be performed\n"); fprintf(stdout, " The current solvent mask is "); printAtomMask(stdout, oldstate->solventMask, oldstate); fprintf(stdout, "\n"); } 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(stdout, " 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 */ if (box[0]!=0.0 && box[1]!=0.0 && box[2]!=0.0) { min = box[0] + box[1] + box[2]; min = min * min; } else { // No box information, set to arbitrarily large max min=DBL_MAX; } 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, box, (double *) ucell, (double *) recip, 0.0, action->iarg3); if ( ! first_only ) { for (ii = oldstate->solventMoleculeStart[i]+1; ii < oldstate->solventMoleculeStop[i]; ii++) { distance = calculateDistance2(ii, j, x, y, z, box, (double *) ucell, (double *) recip, 0.0, action->iarg3); 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 ************************************************************ * * transformContacts() --- perform referenced contact calculation * ******************************************************************************/ int transformContacts(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { char *name = "contacts"; stackType **argumentStackPointer; char *buffer, *buffer2, *buffer3, *buffer2ptr, *buffer3ptr; transformContactsInfo *contactsInfo; double x1,y1,z1,x2,y2,z2,dist; int i, j; int *activeResidues; contactList *list; contactList *current; int *contactNumberList; int *nativeNumberList; int n, nativeNumber, contactNumber,resNum; /* * USAGE * * contacts [first | reference] [byresidue] * [out ] [time ] [distance ] [] * * action argument usage: * * byresidue: calculate number of contacts for every specified atom and save result per residue * iarg1: transformContactsType * CONTACTS_FIRST -- take first structure for reference contacts * CONTACTS_REFERENCE -- take reference structure for reference contacts * iarg2: frame counter * iarg3: stores number of atoms * darg1: time interval * darg2: cutoff distance * carg1: the contactsInfo structure * carg2: list of all reference contacts * carg3: list of active residues */ /* Set up buffers for printing */ buffer2 = safe_malloc( BUFFER_SIZE * sizeof *buffer2); buffer3 = safe_malloc( BUFFER_SIZE * sizeof *buffer3); buffer2ptr = buffer2; buffer3ptr = buffer3; if (mode == PTRAJ_SETUP) { /* * ACTION PTRAJ_SETUP */ #ifdef MPI #endif nativeNumberList = NULL; contactNumberList = NULL; argumentStackPointer = (stackType **) action->carg1; action->carg1 = NULL; contactsInfo = (transformContactsInfo *) safe_malloc(sizeof(transformContactsInfo)); INITIALIZE_transformContactsInfo(contactsInfo); action->iarg1 = (int) CONTACTS_FIRST; if (argumentStackContains(argumentStackPointer, "reference")) action->iarg1 = (int) CONTACTS_REFERENCE; else if (argumentStackContains(argumentStackPointer, "first")) action->iarg1 = (int) CONTACTS_FIRST; contactsInfo->filename = argumentStackKeyToString(argumentStackPointer, "out", NULL); if (contactsInfo->filename != NULL){ contactsInfo->outFile = ptrajOpenW(contactsInfo->filename); } else{ contactsInfo->outFile = stdout; } if (argumentStackContains(argumentStackPointer, "byresidue")) contactsInfo->byResidue = 1; action->iarg2 = 0; action->iarg3 = action->state->atoms; action->carg1 = contactsInfo; action->darg1 = argumentStackKeyToDouble(argumentStackPointer, "time", 1.0); action->darg2 = argumentStackKeyToDouble(argumentStackPointer, "distance", 7.0); /* * Get mask (here, everything else should have been processed from the argumentStack) */ buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { if (contactsInfo->byResidue) action->mask = processAtomMask("@CA", action->state); else action->mask = processAtomMask("*", action->state); } else { action->mask = processAtomMask(buffer, action->state); safe_free(buffer); } /* * Do final setup for which mask information is needed */ nativeNumber = 0; if (action->iarg1 == (int) CONTACTS_REFERENCE) { // Check that a reference structure has been defined if (referenceInfo == NULL) { fprintf(stdout,"Error: No reference coordinates defined.\n"); return -1; } n = action->state->atoms; list = (contactList *) safe_malloc(sizeof(contactList) * n); for(i = 0; i < n; i++){ list[i].index = -1; list[i].name = NULL; list[i].next = NULL; } for (i = 0; i < n; i++) { if (action->mask[i]) { list[i].index = i; list[i].name = action->state->atomName[i]; list[i].next = NULL; for (j = 0; j < n; j++) { if (action->mask[j] && i != j) { x1 = referenceInfo->x[i]; y1 = referenceInfo->y[i]; z1 = referenceInfo->z[i]; x2 = referenceInfo->x[j]; y2 = referenceInfo->y[j]; z2 = referenceInfo->z[j]; dist = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2)); if (dist < action->darg2) { current = (contactList *) safe_malloc(sizeof(contactList)); INITIALIZE_contactList(current); current->name = action->state->atomName[i]; current->index = j; current->next = list[i].next; list[i].next = current; nativeNumber++; } } } } } action->carg2 = (contactList *) list; } if(contactsInfo->byResidue){ contactsInfo->outFile2 = ptrajOpenW(strncat(contactsInfo->filename, ".native", 7)); sprintf(buffer2, "#time"); buffer2 = buffer2ptr + strlen(buffer2ptr); activeResidues = safe_malloc(sizeof(int)*action->state->residues); for (i = 0; i < action->state->residues; i++) activeResidues[i] = 0; for (i = 0; i < action->state->atoms; i++) { if (action->mask[i]) { resNum = atomToResidue(i+1, action->state->residues, action->state->ipres) - 1; activeResidues[resNum] = 1; } } for (i = 0; i < action->state->residues; i++) if (activeResidues[i]) { sprintf(buffer2, "\tresidue %d", i); buffer2 = buffer2ptr + strlen(buffer2ptr); } sprintf(buffer2, "\n"); buffer2 = buffer2ptr; ptrajfprintfone(contactsInfo->outFile, buffer2); ptrajfprintfone(contactsInfo->outFile2, buffer2); action->carg3 = (void *) activeResidues; } else { sprintf(buffer2, "#time\tContacts\tnative Contacts "); buffer2 = buffer2ptr + strlen(buffer2ptr); if (action->iarg1 == (int) CONTACTS_REFERENCE) { sprintf(buffer2, "(number of natives: %d)", nativeNumber); buffer2 = buffer2ptr + strlen(buffer2ptr); } sprintf(buffer2, "\n"); buffer2 = buffer2ptr; ptrajfprintfone(contactsInfo->outFile, buffer2); } } else if (mode == PTRAJ_STATUS) { /* * ACTION PTRAJ_STATUS */ contactsInfo = (transformContactsInfo *) action->carg1; fprintf(stdout, " CONTACTS: Calculating current contacts and comparing results to "); if (action->iarg1 == (int) CONTACTS_FIRST) fprintf(stdout, "first frame\n"); else if (action->iarg1 == (int) CONTACTS_REFERENCE ) fprintf(stdout, "reference structure\n"); fprintf(stdout, " Dumping results to %s,\n", contactsInfo->outFile == NULL ? "STDOUT" : contactsInfo->filename); fprintf(stdout, " using a time interval of %f\n", action->darg1); fprintf(stdout, " and a cutoff of %f A\n", action->darg2); if(contactsInfo->byResidue) fprintf(stdout, " Results are output on a per-residue basis\n"); fprintf(stdout, " Atom selection follows "); printAtomMask(stdout, action->mask, action->state); fprintf(stdout, "\n"); } else if (mode == PTRAJ_ACTION) { /* * ACTION PTRAJ_ACTION */ contactsInfo = (transformContactsInfo *) action->carg1; activeResidues = (int *) action->carg3; n = action->state->atoms; if (action->iarg1 == (int) CONTACTS_FIRST && action->iarg2 == 0) { list = (contactList * )safe_malloc(sizeof(contactList) * n); for(i = 0; i < n; i++){ list[i].index = -1; list[i].name = NULL; list[i].next = NULL; } for (i = 0; i < n; i++) { if (action->mask[i]) { list[i].index = i; list[i].name = action->state->atomName[i]; list[i].next = NULL; for (j = 0; j < n; j++) { if (action->mask[j] && i != j) { dist = sqrt((x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) + (z[i]-z[j])*(z[i]-z[j])); if (dist < action->darg2) { current = (contactList *) safe_malloc(sizeof(contactList)); INITIALIZE_contactList(current); current->name = action->state->atomName[i]; current->index = j; current->next = list[i].next; list[i].next = current; } } } } } action->carg2 = (contactList *) list; } list = (contactList *) action->carg2; contactNumberList = safe_malloc(sizeof(int) * action->state->residues); nativeNumberList = safe_malloc(sizeof(int) * action->state->residues); for (i = 0; i < action->state->residues; i++) { contactNumberList[i] = 0; nativeNumberList[i] = 0; } if (contactsInfo->byResidue) { for (i = 0; i < n; i++) { if (action->mask[i]) { contactNumber = 0; nativeNumber = 0; resNum = atomToResidue(i+1, action->state->residues, action->state->ipres) - 1; for (j = 0; j < n; j++) { if(action->mask[j] && i != j) { dist = sqrt((x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) + (z[i]-z[j])*(z[i]-z[j])); if (dist < action->darg2) { contactNumber++; contactNumberList[resNum]++; current = &list[i]; current = current->next; while (current != NULL) { if (current->index == j) { nativeNumber++; nativeNumberList[resNum]++; break; } current = current->next; } } } } } } sprintf(buffer2, "%10.2f", (double) (action->iarg2*worldsize+worldrank+1) * action->darg1); sprintf(buffer3, "%10.2f", (double) (action->iarg2*worldsize+worldrank+1) * action->darg1); buffer2 = buffer2ptr + strlen(buffer2ptr); buffer3 = buffer3ptr + strlen(buffer3ptr); for (i = 0; i < action->state->residues; i++) if (activeResidues[i]) { sprintf(buffer2, "\t%d", contactNumberList[i]); sprintf(buffer3, "\t%d", nativeNumberList[i]); buffer2 = buffer2ptr + strlen(buffer2ptr); buffer3 = buffer3ptr + strlen(buffer3ptr); } sprintf(buffer2, "%d\n", contactNumber); sprintf(buffer3, "%d\n", nativeNumber); buffer2 = buffer2ptr; buffer3 = buffer3ptr; ptrajfprintf(contactsInfo->outFile, buffer2); ptrajfprintf(contactsInfo->outFile2, buffer3); } else { contactNumber = nativeNumber = 0; for (i = 0; i < n; i++) { if (action->mask[i]) { for (j = 0; j < n; j++) { if(action->mask[j] && i != j) { dist = sqrt((x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) + (z[i]-z[j])*(z[i]-z[j])); if (dist < action->darg2) { contactNumber++; current = &list[i]; current = current->next; while (current != NULL) { if (current->index == j) { nativeNumber++; break; } current = current->next; } } } } } } sprintf(buffer2, "%10.2f\t%d\t%d\n", (action->iarg2*worldsize+worldrank+1) * action->darg1, contactNumber, nativeNumber); ptrajfprintf(contactsInfo->outFile, buffer2); } action->iarg2++; /* Clean up contactNumberList and nativeNumberList */ safe_free(contactNumberList); safe_free(nativeNumberList); } else if (mode == PTRAJ_PRINT) { /* * ACTION PTRAJ_PRINT */ /* Nothing to do here */ } else if (mode == PTRAJ_CLEANUP) { /* * ACTION PTRAJ_CLEANUP */ contactsInfo = (transformContactsInfo *) action->carg1; if(contactsInfo != NULL){ if (contactsInfo->byResidue){ if(contactsInfo->filename != NULL) { ptrajCloseFile(contactsInfo->outFile2); } safe_free((int *) action->carg3); } if(contactsInfo->filename != NULL){ ptrajCloseFile(contactsInfo->outFile); safe_free(contactsInfo->filename); } INITIALIZE_transformContactsInfo(contactsInfo); safe_free(contactsInfo); } /* * This should be done each time in PTRAJ_ACTION * * safe_free(contactNumberList); * safe_free(nativeNumberList); */ list = (contactList *) action->carg2; if(list != NULL){ for(i = 0; i < action->iarg3; i++){ current = list + i; current = current->next; while (current != NULL) { contactList *next = current->next; INITIALIZE_contactList(current); safe_free(current); current = next; } } safe_free(list); } } 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) { char *name = "correlation"; 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 */ #ifdef MPI printParallelError(name); return -1; #endif 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(stdout, corrInfo->mask, action->state); fprintf(stdout, "\n"); fprintf(stdout, " Atom selection 2 is "); printAtomMask(stdout, corrInfo->mask2, action->state); fprintf(stdout, "\n"); } 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) { char *name = "dihedral"; stackType **argumentStackPointer; char *buffer, buffer2[BUFFER_SIZE]; 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; void *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 */ #ifdef MPI #endif 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(stdout, "WARNING: ptraj(), dihedral: It is necessary to specify a unique name\n"); fprintf(stdout, "for each angle specified. Ignoring command...\n"); safe_free(info); return -1; } else if ( scalarStackGetName(&scalarStack, info->name) != NULL ) { fprintf(stdout, "WARNING: ptraj(), dihedral: The chosen name (%s) has already been used.\n", info->name); fprintf(stdout, "Ignoring command...\n"); safe_free(info); return -1; } info->state = action->state; /* * grab the type if present */ buffer = argumentStackKeyToString(argumentStackPointer, "type", NULL); if (buffer != NULL) { if (strcmp(buffer, "alpha") == 0) info->type = SCALAR_TYPE_ALPHA; else if (strcmp(buffer, "beta") == 0) info->type = SCALAR_TYPE_BETA; else if (strcmp(buffer, "gamma") == 0) info->type = SCALAR_TYPE_GAMMA; else if (strcmp(buffer, "delta") == 0) info->type = SCALAR_TYPE_DELTA; else if (strcmp(buffer, "epsilon") == 0) info->type = SCALAR_TYPE_EPSILON; else if (strcmp(buffer, "zeta") == 0) info->type = SCALAR_TYPE_ZETA; else if (strcmp(buffer, "chi") == 0) info->type = SCALAR_TYPE_CHI; else if (strcmp(buffer, "c2p") == 0) info->type = SCALAR_TYPE_C2P; else if (strcmp(buffer, "h1p") == 0) info->type = SCALAR_TYPE_H1P; else if (strcmp(buffer, "phi") == 0) info->type = SCALAR_TYPE_PHI; else if (strcmp(buffer, "psi") == 0) info->type = SCALAR_TYPE_PSI; else if (strcmp(buffer, "pchi") == 0) info->type = SCALAR_TYPE_PCHI; safe_free(buffer); } /* * grab the output filename, if specified */ info->filename = argumentStackKeyToString(argumentStackPointer, "out", NULL); /* * push the distance info on to the scalar stack */ pushBottomStack(&scalarStack, (void *) info); action->darg1 = argumentStackKeyToDouble(argumentStackPointer, "time", 1.0); /* * process mask1 --> mask4 */ buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { fprintf(stdout, "WARNING in ptraj(), dihedral: Error in specification of the first mask: name %s\n", info->name); fprintf(stdout, "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(stdout, "WARNING in ptraj(), dihedral: Error in specification of the second mask: name %s\n", info->name); fprintf(stdout, "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(stdout, "WARNING in ptraj(), dihedral: Error in specification of the third mask: name %s\n", info->name); fprintf(stdout, "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(stdout, "WARNING in ptraj(), dihedral: Error in specification of the fourth mask: name %s\n", info->name); fprintf(stdout, "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(stdout, "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(stdout, "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(stdout, "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(stdout, "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) { if (prnlev < 1) return 0; /* * 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(stdout, info->mask1, action->state); fprintf(stdout, "\n"); } 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(stdout, info->mask2, action->state); fprintf(stdout, "\n"); } 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(stdout, info->mask3, action->state); fprintf(stdout, "\n"); } 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(stdout, info->mask4, action->state); fprintf(stdout, "\n"); } 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 ( info->filename != NULL ) { outFile = ptrajOpenW(info->filename); if ( outFile == NULL ) { fprintf(stdout, "WARNING in ptraj(), dihedral: couldn't open file %s\n", info->filename); return 0; } if (prnlev > 2) fprintf(stdout, "PTRAJ DIHEDRAL dumping named values %s\n", info->name); for (i=0; i < action->state->maxFrames/worldsize; i++) { ptrajfprintf(outFile, "%10.2f %f\n", (i*worldsize+worldrank+1) * action->darg1, info->value[i]); } ptrajCloseFile(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 **************************************************** * * transformDihedralCluster() --- Cluster trajectory by dihedral angles * * Written by Dan Roe, Stonybrook U. * ******************************************************************************/ /*NOTE: Place checks after the realloc statements*/ /* DCbin: recursive function to place diedral bin combinations in a Tree * structure - this conserves memory and has better search performance than * a linear array. */ int DCbin(int* Bins,int level, DCnodetype* T, int max, double FRAME,long int* numCluster) { long int i,j; FILE *outfile; /*DEBUG Output*/ if (prnlev>1) { if (level==0) printf("DIHEDRAL CLUSTER, FRAME %lf\n",FRAME); for (i=0; i<(level*2); i++) printf(" "); printf("Level %i Bin %i ",level,Bins[level]); } /*Is this the first time bin has been called?*/ if (T->numBranch==0) { if (prnlev>1) printf("Creating. "); T->numBranch=1; T->bin=(int*) safe_malloc(sizeof(int)); T->bin[0]=Bins[level]; /*If this is the lowest level allocate count and frames*/ if (level==max) { T->branch=NULL; T->count=(long int*) safe_malloc(sizeof(long int)); T->count[0]=1; (*numCluster)++; T->frames=(double**) safe_malloc(sizeof(double*)); T->frames[0]=(double*) safe_malloc(sizeof(double)); T->frames[0][0]=FRAME; if (prnlev>1) printf("Count= %li, Frame= %.0lf\n",T->count[0],T->frames[0][0]); return 0; /*Otherwise create a branch node for the next bin*/ } else { T->count=NULL; T->frames=NULL; T->branch=(DCnodetype**) safe_malloc(sizeof(DCnodetype*)); T->branch[0]=(DCnodetype*) safe_malloc(sizeof(DCnodetype)); INITIALIZE_transformDihedralCluster(T->branch[0]); if (prnlev>1) printf("Next.\n"); DCbin(Bins,level+1,T->branch[0],max,FRAME,numCluster); return 0; } } else { /*If not first call, does this Bin already exist?*/ if (prnlev>1) printf("Searching. "); for (i=0; inumBranch; i++) { if (T->bin[i]==Bins[level]) { /*If it does and we're at lowest level, increment count, record frame*/ if (level==max) { T->count[i]++; j=T->count[i]; T->frames[i]=(double*) realloc(T->frames[i],j*sizeof(double)); T->frames[i][j-1]=FRAME; if (prnlev>1) printf("Count= %li, Frame= %.0lf\n",j,T->frames[i][j-1]); return 0; } else { /*If it does and we're not at lowest level, continue search*/ if (prnlev>1) printf("Next.\n"); DCbin(Bins,level+1,T->branch[i],max,FRAME,numCluster); return 0; } } } /*Bin doesnt exist, create a new branch*/ if (prnlev>1) printf("Newbranch. "); T->numBranch++; i=T->numBranch; T->bin=(int*) realloc(T->bin,i*sizeof(int)); T->bin[i-1]=Bins[level]; if (level==max) { /*If lowest level, set count and frame number*/ T->branch=NULL; T->count=(long int*) realloc(T->count,i*sizeof(long int)); T->count[i-1]=1; (*numCluster)++; T->frames=(double**) realloc(T->frames,i*sizeof(double*)); T->frames[i-1]=(double*) safe_malloc(sizeof(double)); T->frames[i-1][0]=FRAME; if (prnlev>1) printf("Count= %li, Frame= %.0lf\n",T->count[i-1],T->frames[i-1][0]); return 0; } else { /*Otherwise, continue down the branch*/ T->count=NULL; T->frames=NULL; T->branch=(DCnodetype**) realloc(T->branch,i*sizeof(DCnodetype*)); T->branch[i-1]=(DCnodetype*) safe_malloc(sizeof(DCnodetype)); INITIALIZE_transformDihedralCluster(T->branch[i-1]); if (prnlev>1) printf("Next.\n"); DCbin(Bins,level+1,T->branch[i-1],max,FRAME,numCluster); return 0; } } return 1; } /* freeDCbin: recursive function to free memory allocated to Tree */ void freeDCbin(DCnodetype* T) { long int i; if (T->branch!=NULL) for (i=0; inumBranch; i++) freeDCbin(T->branch[i]); if (T->count!=NULL) safe_free(T->count); if (T->frames!=NULL) { for (i=0; inumBranch; i++) safe_free(T->frames[i]); safe_free(T->frames); } safe_free(T->bin); safe_free(T->branch); return; } /* DCtree2array: recursive function to put contents of tree into array for sorting. */ void DCtree2array(int* Bins, int level, DCnodetype* T,long int* numCluster,DCarray** A) { long int i,j,k; for (i=0; inumBranch; i++) { Bins[level]=T->bin[i]; if (T->branch!=NULL) DCtree2array(Bins,level+1,T->branch[i],numCluster,A); else { (*numCluster)++; k=(*numCluster); /*Store cluster info in array for sorting later*/ if (prnlev>1) printf("DEBUG: Writing cluster %li\n",k); /*Assign Bins*/ if (prnlev>1) printf("DEBUG: Bins= "); for (j=0; j<=level; j++) { A[k-1]->Bins[j]=Bins[j]; if (prnlev>1) printf("%i ",A[k-1]->Bins[j]); } if (prnlev>1) printf("\n"); /*Assign count*/ A[k-1]->count=T->count[i]; if (prnlev>1) printf("DEBUG: Count= %li\n",A[k-1]->count); /*Assign frames*/ A[k-1]->frames=(double*) safe_malloc(T->count[i]*sizeof(double)); if (prnlev>1) printf("DEBUG: frames= "); for (j=0; jcount[i]; j++){ A[k-1]->frames[j]=T->frames[i][j]; if (prnlev>1) printf("%.0lf ",A[k-1]->frames[j]); } if (prnlev>1) printf("\n"); /*printf("DEBUG: Current Counts\n"); for (j=0; jcount);*/ } } return; } /*compare: used by the qsort function in sorting cluster array */ int compare(const void *a, const void *b) { DCarray *A; DCarray *B; A=*(DCarray**) a; B=*(DCarray**) b; /*printf(" QSORT DEBUG: Comparing %i to %i\n",A->count,B->count);*/ return ( B->count - A->count ); } int transformDihedralCluster(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { char *name = "clusterdihedral"; stackType **argumentStackPointer; ptrajState *state; char *buffer = NULL; /* Usage: clusterdihedral [phibins N] number of bins in phi direction [psibins M] number of bins in psi direction Note: phibins and psibins only used if dihedralfile not specified [out ] file to print cluster information to [cut CUT] only print clusters with pop > CUT [framefile ] file to print frame-cluster info [clusterinfo ] print cluster info in format sander can read for NB weighted RREMD [dihedralfile] read dihedral definitions from FILE. Format is ATOM#1 ATOM#2 ATOM#3 ATOM#4 BINS Note: This functionality treats atom numbers as starting from 1 to be consistent with sander. [] if not reading dihedrals from file Backbone dihedrals will be searched for within MASK. Action argument usage: iarg1 = number of dihedral angles to keep track of iarg2 = if 1, dihedral angles were read from a file, otherwise backbone dihedrals were searched for in MASK. iarg4 = store CUT darg3 = keep track of number of frames, used for memory allocation at end darg4 = number of clusters carg1 = int array; atom masks of each dihedral and the number of bins carg2 = hold the DCnodetype tree carg3 = int array; during run, store which bin each dihedral falls into Note: Even though no information from carg3 needs to be carried over from frame to frame, it is stored here to avoid reallocating the memory each time. carg4 = char* array to hold output filenames */ DCnodetype* T; DCarray** A; int** DCmasks; int* Bins; long int* framecluster; long int i,j,k; int C1,N2,CA,C2; int phibins,psibins; int numDihedral,CUT; long int numCluster; double FRAME; double phistep, PHI, temp; double cx1, cy1, cz1; double cx2, cy2, cz2; double cx3, cy3, cz3; double cx4, cy4, cz4; char** filenames; FILE* outfile; state = (ptrajState *) action->state; /* * ---=== PTRAJ_SETUP ===--- */ if (mode == PTRAJ_SETUP) { #ifdef MPI printParallelError(name); return -1; #endif argumentStackPointer = (stackType **) action->carg1; action->carg1 = NULL; /*Parse Command Line*/ phibins=argumentStackKeyToInteger(argumentStackPointer, "phibins", 10); if ((phibins>360)||(phibins<=1)) { fprintf(stdout,"clusterdihedral Error: phibins out of range 360 < x <= 1 (%i)\n",phibins); return -1; } psibins=argumentStackKeyToInteger(argumentStackPointer, "psibins", 10); if ((psibins>360)||(psibins<=1)) { fprintf(stdout,"clusterdihedral Error: psibins out of range 360 < x <= 1 (%i)\n",psibins); return -1; } /*Cluster cutoff*/ CUT=argumentStackKeyToInteger(argumentStackPointer, "cut",0); /*Output Files*/ filenames=(char**) safe_malloc(4*sizeof(char*)); buffer=argumentStackKeyToString(argumentStackPointer, "out",NULL); if (buffer!=NULL) { filenames[0]=(char*) safe_malloc((strlen(buffer)+1)*sizeof(char)); strcpy(filenames[0],buffer); } else filenames[0]=NULL; buffer=argumentStackKeyToString(argumentStackPointer, "framefile",NULL); if (buffer!=NULL) { filenames[1]=(char*) safe_malloc((strlen(buffer)+1)*sizeof(char)); strcpy(filenames[1],buffer); } else filenames[1]=NULL; buffer=argumentStackKeyToString(argumentStackPointer, "clusterinfo",NULL); if (buffer!=NULL) { filenames[2]=(char*) safe_malloc((strlen(buffer)+1)*sizeof(char)); strcpy(filenames[2],buffer); } else filenames[2]=NULL; buffer=argumentStackKeyToString(argumentStackPointer, "clustervtime",NULL); if (buffer!=NULL) { filenames[3]=(char*) safe_malloc((strlen(buffer)+1)*sizeof(char)); strcpy(filenames[3],buffer); } else filenames[3]=NULL; action->carg4=(void*) filenames; /*Input Dihedral file*/ action->iarg2=0; numDihedral=0; buffer=argumentStackKeyToString(argumentStackPointer, "dihedralfile",NULL); if (buffer!=NULL) { if ( (outfile=safe_fopen(buffer,"r"))==NULL ) { fprintf(stdout,"WARNING: Could not open dihedralfile %s",buffer); } else { /*Read Dihedrals and bins from a file*/ printf("Reading Dihedrals from %s\n",buffer); DCmasks=(int**) safe_malloc(sizeof(int*)); while ( fscanf(outfile,"%i %i %i %i %i",&C1,&N2,&CA,&C2,&phibins) != EOF) { numDihedral++; DCmasks=(int**) realloc(DCmasks,numDihedral*sizeof(int*)); if (DCmasks==NULL) { fprintf(stdout,"clusterdihedral Error: Memory reallocation for masks failed.\n"); return -1; } DCmasks[numDihedral-1]=(int*) safe_malloc(5*sizeof(int)); /*amber atom numbers start from 1*/ DCmasks[numDihedral-1][0]=C1-1; DCmasks[numDihedral-1][1]=N2-1; DCmasks[numDihedral-1][2]=CA-1; DCmasks[numDihedral-1][3]=C2-1; DCmasks[numDihedral-1][4]=phibins; printf("(%i)-(%i)-(%i)-(%i) Bins=%i\n",C1,N2,CA,C2,phibins); } printf("Read %i dihedrals\n",numDihedral); action->iarg2=1; safe_fclose(outfile); } } /*Allocate Memory to hold cluster info*/ action->darg4=0; T=(DCnodetype*) safe_malloc(sizeof(DCnodetype)); INITIALIZE_transformDihedralCluster(T); FRAME=0; /* Process Mask */ buffer = getArgumentString(argumentStackPointer, NULL); if (buffer == NULL) { action->mask = processAtomMask("*", action->state); } else { action->mask = processAtomMask(buffer, action->state); safe_free(buffer); } /*Set up backbone dihedral angles if none were read*/ if (numDihedral==0) { printf("Scanning for backbone dihedrals.\n"); C1=-1; N2=-1; CA=-1; C2=-1; DCmasks=(int**) safe_malloc(sizeof(int*)); for (i=0; i < state->atoms; i++) { if (action->mask[i]==1) { /*printf(" DEBUG: Atom %i: %s\n",i,state->atomName[i]);*/ if (C2>-1) { /* If we have already found the last C in phi dihedral, this N is * the last atom in psi dihedral - store both. */ if ( strcmp(state->atomName[i],"N ")==0 ) { numDihedral+=2; /* dynamically grow array */ DCmasks=(int**) realloc(DCmasks,numDihedral*sizeof(int*)); if (DCmasks==NULL) { fprintf(stdout,"clusterdihedral Error: Memory reallocation for masks failed.\n"); return -1; } DCmasks[numDihedral-2]=(int*) safe_malloc(5*sizeof(int)); DCmasks[numDihedral-2][0]=C1; DCmasks[numDihedral-2][1]=N2; DCmasks[numDihedral-2][2]=CA; DCmasks[numDihedral-2][3]=C2; DCmasks[numDihedral-2][4]=phibins; DCmasks[numDihedral-1]=(int*) safe_malloc(5*sizeof(int)); DCmasks[numDihedral-1][0]=N2; DCmasks[numDihedral-1][1]=CA; DCmasks[numDihedral-1][2]=C2; DCmasks[numDihedral-1][3]=i; DCmasks[numDihedral-1][4]=psibins; if (prnlev>0) printf("DIHEDRAL PAIR FOUND: C1= %i, N2= %i, CA= %i, C2= %i, N3= %li\n", C1,N2,CA,C2,i); /* Since the carbonyl C/amide N probably starts a new dihedral, * reset to those */ C1=C2; N2=i; C2=-1; CA=-1; } } else if (C1>-1) { /* If we've already found the first carbonyl, look for other atoms * in the dihedral pair. */ if ( strcmp(state->atomName[i],"N ")==0 ) N2=i; if ( strcmp(state->atomName[i],"CA ")==0 ) CA=i; if ( strcmp(state->atomName[i],"C ")==0 ) C2=i; } else if ( strcmp(state->atomName[i],"C ")==0 ) C1=i; /*1st carbon*/ } } } if ( numDihedral == 0 ) { fprintf(stdout,"clusterdihedral Error: No Backbone Dihedral Angles Found!\n"); return -1; } if (prnlev>0) printf("FOUND %i DIHEDRAL ANGLES.\n",numDihedral); /* Allocate memory to store dihedral bins so we don't continuously * reallocate during the ptraj run. */ Bins=(int*) safe_malloc(numDihedral*sizeof(int)); /* Assign action pointers */ action->iarg1 = numDihedral; action->iarg4 = CUT; action->carg1 = (void*) DCmasks; action->carg2 = (void*) T; action->carg3 = (void*) Bins; action->darg3 = FRAME; return 0; } /* * ---=== PTRAJ_CLEANUP ===--- */ if (mode == PTRAJ_CLEANUP) { numDihedral=action->iarg1; DCmasks=(int**) action->carg1; T=(DCnodetype*) action->carg2; Bins=(int*) action->carg3; filenames=(char**) action->carg4; if (filenames[0]!=NULL) safe_free(filenames[0]); if (filenames[1]!=NULL) safe_free(filenames[1]); if (filenames[2]!=NULL) safe_free(filenames[2]); if (filenames[3]!=NULL) safe_free(filenames[3]); if (prnlev>1) printf("DIHEDRAL: Freeing %i dihedral masks.\n",numDihedral); for (i=0; iiarg1; DCmasks=(int**) action->carg1; CUT=action->iarg4; filenames=(char**) action->carg4; fprintf(stdout,"\n DIHEDRAL CLUSTERING on %i dihedral angles\n",numDihedral); if (filenames[0]!=NULL) fprintf(stdout," Cluster data will be output to %s\n",filenames[0]); if (CUT>0) fprintf(stdout," Only clusters with pop > %i will be printed.\n",CUT); if (filenames[1]!=NULL) fprintf(stdout," Frame-Cluster data will be output to %s\n",filenames[1]); if (filenames[2]!=NULL) fprintf(stdout," Cluster information (pop. & ID) will be output to %s\n",filenames[2]); if (filenames[3]!=NULL) fprintf(stdout," Number of clusters v time will be output to %s\n",filenames[3]); if (action->iarg2==0) { fprintf(stdout, " Looked for dihedrals within atom selection= "); printAtomMask(stdout, action->mask, action->state); fprintf(stdout, "\n"); } else fprintf(stdout," Read dihedrals from an input file.\n"); for (i=0; iatomName[C1]); } fprintf(stdout," [Bins=%i]",DCmasks[i][4]); phistep=360/DCmasks[i][4]; fprintf(stdout," [Step=%lf]\n",phistep); } fprintf(stdout,"\n"); return 0; } /* * ---=== PTRAJ_PRINT ===--- */ if (mode == PTRAJ_PRINT) { DCmasks=(int**) action->carg1; filenames=(char**) action->carg4; numCluster=(long int) action->darg4; numDihedral=action->iarg1; CUT=action->iarg4; T=(DCnodetype*) action->carg2; Bins=(int*) action->carg3; /*Set up output file*/ if (filenames[0]!=NULL) { if ( (outfile=safe_fopen(filenames[0],"w"))==NULL ) { fprintf(stdout,"Could not open %s, reverting to stdout.",filenames[0]); outfile=stdout; } } else outfile=stdout; /*Print Bin information*/ printf("Printing Dihedral Clustering Results.\n"); fprintf(outfile,"DIHEDRAL CLUSTER RESULTS"); if (action->iarg2==0) { fprintf(outfile," for "); printAtomMask(outfile, action->mask, action->state); fprintf(outfile,"\n"); } else fprintf(outfile,"\n"); if (outfile!=stdout) { for (i=0; iatomName[C1],C1+1); } fprintf(outfile," [Bins=%i]\n",DCmasks[i][4]); } } fprintf(outfile,"%li clusters.\n",numCluster); if (CUT>0) fprintf(outfile,"Only printing clusters with pop > %i\n",CUT); /* fprintf(outfile,"Phi Bin Ranges:\n"); for (i=0; iBins=(int*) safe_malloc(numDihedral*sizeof(int)); } /*Place tree elements into an array for sorting, then sort*/ i=0; DCtree2array(Bins,0,T,&i,A); if (prnlev>1) printf("%li clusters written.\n",i); qsort(A,numCluster,sizeof(DCarray*),compare); /*Allocate memory for printing frame/cluster */ i=(long int)action->darg3; framecluster=(long int*) safe_malloc((i+1)*sizeof(long int)); /*Print sorted cluster array*/ for (j=0; jcount>CUT) { fprintf(outfile,"Cluster %10li %10li [ ",j,A[j]->count); for (i=0; iBins[i]); fprintf(outfile," ]\n"); for (i=0; icount; i++) { fprintf(outfile,"%.0lf ",A[j]->frames[i]); /* store which cluster each frame belongs to. Not neccesary if user * didn't specify this option, but avoids a second loop if they did. */ k=(long int) A[j]->frames[i]; framecluster[k]=j; } fprintf(outfile,"\n"); } } if (buffer!=NULL) safe_fclose(outfile); /*cluster for each frame*/ if (filenames[1]!=NULL) { printf("Printing cluster number for each frame.\n"); if ( (outfile=safe_fopen(filenames[1],"w"))==NULL ) { fprintf(stdout,"WARNING: Could not open framefile %s",filenames[1]); } else { i=(long int) action->darg3; for (j=1; j<=i; j++) { C1=framecluster[j]; fprintf(outfile,"%10li %10i %10li ",j,C1,A[C1]->count); for (k=0; kBins[k]); fprintf(outfile,"\n"); } safe_fclose(outfile); } } /*cluster information file*/ if (filenames[2]!=NULL) { printf("Printing cluster information.\n"); if ( (outfile=safe_fopen(filenames[2],"w"))==NULL ) { fprintf(stdout,"WARNING: Could not open clusterinfo file %s",filenames[2]); } else { fprintf(outfile,"%i\n",numDihedral); for (i=0; icount); for (k=0; kBins[k]); fprintf(outfile,"\n"); } safe_fclose(outfile); } } /*Cleanup arrays*/ for (j=0; jBins); safe_free(A[j]->frames); } safe_free(A); safe_free(framecluster); return 0; } /* * ---=== PTRAJ_ACTION ===--- */ if (mode == PTRAJ_ACTION) { numDihedral=action->iarg1; DCmasks=(int**) action->carg1; T=(DCnodetype*) action->carg2; Bins=(int*) action->carg3; FRAME=action->darg3; FRAME++; numCluster=(long int) action->darg4; /*For each dihedral, calculate which bin it should go into and store bin#*/ j=0; for (i=0; i0) printf("%9s%8.2lf ","Dihedral=",PHI); PHI+=180; phistep=360/DCmasks[i][4]; PHI/=phistep; modf(PHI,&temp); phibins=temp; if (prnlev>0) printf("%4s%3i\n","Bin=",phibins); Bins[j]=phibins; j++; } /* At this point j=numDihedral */ if (prnlev>0) { printf("["); for (i=0; idarg4=(double) numCluster; action->carg2=(void*) T; action->darg3=FRAME; return 1; } return 0; } /*DAN ROE*/ /** ACTION ROUTINE ************************************************************* * * transformDiffusion() --- calculate mean squared displacements vs. time * ******************************************************************************/ int transformDiffusion(actionInformation *action, double *x, double *y, double *z, double *box, int mode) { char *name = "diffusion"; 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