/* _______________________________________________________________________ * * RDPARM/PTRAJ: APR 2000 * _______________________________________________________________________ * * $Header: /thr/gamow/cvsroot/amber7/src/ptraj/ptraj.c,v 1.28 2001/07/13 21:59:19 cheatham Exp $ * * Revision: $Revision: 1.28 $ * Date: $Date: 2001/07/13 21:59:19 $ * Last checked in by $Author: cheatham $ * * This source is now archived under CVS at Scripps in the amber7 tree. * * NOTE: this is a beta, pre-release version, is under constant development, * probably contains some bugs and is under constant revision; therefore * take care. Please report bugs (and fixes!) to cheatham@chpc.utah.edu or * cheatham@cgl.ucsf.edu * * Do not distribute this code without explicit permission. * Do not incorporate into other codes without explicit permission. * * Current contact information: * * Thomas Cheatham, III * 2000 East, 30 South, Skaggs Hall 201 * Department of Medicinal Chemistry * University of Utah * Salt Lake City, UT 84112-5820 * cheatham@chpc.utah.edu * (801) 587-9653 * FAX: (801) 585-5366 * * Other contributors: * * David Case (Scripps) * Michael Crowley (Scripps) * Jed Pitera (UCSF) * Vickie Tsui (Scripps) * _______________________________________________________________________ */ #include #include #include #define PTRAJ_MODULE #include "ptraj.h" /* * This is the main routine for the ptraj functionality and in general this * code is likely not to be modified by general users/developers, except: * * (*) Adding new "actions" or routines to process/analyze coordinate sets. * In this case, ptrajSetup() will likely need modification for any * "action" on the coordinates or ptrajSetupAnalyze() for any routine that * proprocesses data accumulated by the actions. * (*) Adding new coordinate formats. Modify checkCoordinates(), * ptrajSetupIO(), printCoordinateInfo(), ptrajPreProcessInputCoordinates(), * and ptrajProcessInputCoordinates() * (*) Modifying the code to parse "mask" strings: all the *mask* routines * * For more information on adding new actions, see the details comments in * actions.c. */ int atomToResidue(int atom, int residues, int *ipres) { int i; if (ipres == NULL) return -1; for (i = 0; i < residues; i++) if (atom >= ipres[i] && atom < ipres[i+1]) return (i+1); return -1; } int atomToSolventMolecule(int atom, int molecules, int *start, int *stop) { int i; if (start == NULL || stop == NULL) return -1; for (i = 0; i < molecules; i++) if (atom <= start[i]) return -1; else if (atom > start[i] && atom <= stop[i]) return (i+1); return -1; } int atomToMolecule(int atom, int molecules, int *mol) { int i, a; if (mol == NULL) return -1; a = 0; for (i = 0; i < molecules; i++) { a += mol[i]; if (atom <= a) return (i+1); } return -1; } int isActiveDetailed(int atom, int residue, int *mask, int atoms, int residues, Name *atomName, Name *residueName, int *ipres) { int i; if (residue >= residues || residue < 0) { printf("WARNING: residue out of range in isActiveDetailed, res %i (total %i)\n", residue, residues); return 0; } for (i = ipres[residue]-1; i < ipres[residue+1]-1; i++) if ( mask[i] && strcmp(atomName[i], atomName[atom]) == 0 ) return 1; return 0; } int isActive(int atom, int residue, int *mask, ptrajState *state) { return( isActiveDetailed(atom, residue, mask, state->atoms, state->residues, state->atomName, state->residueName, state->ipres) ); } int isActiveResidueDetailed(int residue, int *mask, int atoms, int residues, Name *atomName, Name *residueName, int *ipres) { int i, total; total = 0; if (residue >= residues || residue < 0) { printf("WARNING: residue out of range in isActiveResidueDetailed(), res %i (total %i)\n", residue, residues); return 0; } for (i = ipres[residue]-1; i < ipres[residue+1]-1; i++) if ( ! mask[i] ) return 0; else total++; return total; } int isActiveResidue(int residue, int *mask, ptrajState *state) { return( isActiveResidueDetailed(residue, mask, state->atoms, state->residues, state->atomName, state->residueName, state->ipres) ); } void printAtomMaskDetailed(int *mask, int atoms, int residues, Name *atomName, Name *residueName, int *ipres) { int i, j, curres; char tmpatom[20]; int *resactive, *ressimilar; int printed, numactive, numresactive, numressimilar; int incurres, innextres; printed = 0; numactive = 0; numresactive = 0; numressimilar = 0; /* * This routine is kind of junky since in general I want to avoid printing * as much detail as possible. Therefore, we check to see is certain ranges * residues have all atoms active, etc. to avoid printing each atom in a residue. * This makes it ugly and obtuse. */ if (mask == NULL) { printf("[No atoms are selected]\n"); return; } j=0; for (i=0; i < atoms; i++) if (mask[i]) j++; if (j == 0) { printf("[No atoms are selected]\n"); return; } /* * check if all atoms are active and if so print an asterisk */ j = 0; for (i=0; i < atoms; i++) { j += mask[i]; } if ( j == atoms ) { fprintf(stdout, " * (All atoms are selected)\n"); return; } numactive = j; /* * determine which residues have all the atoms in that residue active */ resactive = (int *) safe_malloc(sizeof(int) * residues); for (i=0; i < residues; i++) { resactive[i] = 0.0; } curres = 0; j = 0; for (i=0; i < atoms; i++) { if (i == ipres[curres+1]-1) { if (j == ipres[curres+1] - ipres[curres]) { resactive[curres] = 1.0; numresactive++; } j = 0; curres++; } if (mask[i]) j++; } if (j == ipres[curres+1] - ipres[curres]) { resactive[curres] = 1.0; } /* * determine the range over which the residues are fully active */ for (curres = residues-2; curres >= 0; curres--) { if (resactive[curres]) { resactive[curres] += resactive[curres+1]; numresactive--; } } /* * determine ranges over which residues have the same atoms active * as the next residue */ ressimilar = (int *) safe_malloc(sizeof(int) * residues); for (i=0; i < residues; i++) { ressimilar[i] = 0.0; } for (curres = residues-2; curres >=0; curres--) { incurres = 0; innextres = 0; for (i = ipres[curres]-1; i < ipres[curres+1]-1; i++) { if ( mask[i] ) { incurres++; } if (isActiveDetailed(i, curres+1, mask, atoms, residues, atomName, residueName, ipres)) innextres++; } if (incurres && innextres == incurres) { ressimilar[curres] = ressimilar[curres+1] + 1; } else { numressimilar++; } } /* * do the actual printing */ j = 0; for (curres = 0; curres < residues; curres++) { if (resactive[curres] ) { /* * If all of the atoms are active in this residue print either the * residue number or range as appropriate */ if (resactive[curres] > 2) { if (j!=0 && j%10 != 0) fprintf(stdout, ","); fprintf(stdout, ":%i-%i", curres+1, curres+resactive[curres]); curres += resactive[curres]-1; } else { if (j!=0 && j%10 != 0) fprintf(stdout, ","); fprintf(stdout, ":%i", curres+1); } j++; if (j != 0 && j % 10 == 0) { fprintf(stdout, "\n "); j = 0; } } else if (ressimilar[curres]) { /* * If there is a set of residues with a similar atom selection... */ if (ressimilar[curres] > 2) { if (j!=0 && j%10 != 0) fprintf(stdout, ","); fprintf(stdout, ":%i-%i", curres+1, curres+ressimilar[curres]+1); curres += ressimilar[curres]; } else { if (j!=0 && j%10 != 0) fprintf(stdout, ","); fprintf(stdout, ":%i", curres+1); } for (i = ipres[curres]-1; i < ipres[curres+1]-1; i++) { if ( mask[i] ) { if (printed) fprintf(stdout, ","); else { fprintf(stdout, "@"); printed = 1; } strcpy(tmpatom, atomName[i]); fprintf(stdout, "%s", strtok(tmpatom, " ")); } } j++; if (j != 0 && j % 10 == 0) { fprintf(stdout, "\n "); j = 0; } } else { /* * Print individual atoms */ if (numactive > 10 && numressimilar > 10 && numresactive > 10) { fprintf(stdout, "\n "); numactive = 0; } for (i = ipres[curres]-1; i < ipres[curres+1]-1; i++) { if ( mask[i] ) { if (j!=0 && j%10 != 0) fprintf(stdout, ","); strcpy(tmpatom, atomName[i]); fprintf(stdout, ":%i@%s", curres+1, strtok(tmpatom, " ")); j++; } if (j != 0 && j % 10 == 0) { fprintf(stdout, "\n "); j = 0; } } } } if (j!=0 && j%10 != 0) fprintf(stdout, "\n"); safe_free(resactive); safe_free(ressimilar); } void printAtom(FILE *fpout, int atom, ptrajState *state) { int curres; char buffer[50]; curres = atomToResidue(atom+1, state->residues, state->ipres)-1; sprintf(buffer, ":%i ", curres+1); sprintf(buffer+6, "%s %s", state->atomName[atom], state->residueName[curres]); fprintf(fpout, "%s", buffer); } void printAtomCompact(FILE *fpout, int atom, ptrajState *state) { int curres; char buffer[50], *bufferp; curres = atomToResidue(atom+1, state->residues, state->ipres)-1; sprintf(buffer, ":%i", curres+1); bufferp = buffer+strlen(buffer); sprintf(bufferp, "@%s", state->atomName[atom]); fprintf(fpout, "%10s", buffer); } parseEntry * parseToken(char **textp, int operator) { char *text; int i, j; parseEntry *p; text = *textp; for (i=0;; i++) { if (text[i] == (char) 0 || text[i] == ',' || text[i] == ':' || text[i] == '@' || (i>0 && text[i] == '-' && !isalpha(text[i-1])) || isspace(text[i])) break; } p = (parseEntry *) safe_malloc(sizeof(parseEntry)); p->isnumber = 0; if ( i > 0 ) p->token = (char *) safe_malloc(sizeof(char) * (i+1)); for (j=0; j < i; j++) { p->token[j] = text[j]; } p->token[i] = (char) 0; p->operator = PARSE_NOOP; switch ( text[i] ) { case '-': p->operator = PARSE_CONTINUATION; case ',': text++; } if ( isdigit( p->token[0] ) ) { p->isnumber = 1; for (j=1; j < i; j++) { if ( isdigit(p->token[j]) == 0 ) { p->isnumber = 0; } } } text = text+i; skipWhitespace(text); /* * extra special check to handle extra spacing between continuation * operators */ if (text[0] == '-') { p->operator = PARSE_CONTINUATION; text += 1; skipWhitespace(text); } else if (text[0] == ',') { text += 1; skipWhitespace(text); } *textp = text; return(p); } int isMatch(char *s1, char *s2) { int i; /* * straight match */ if ( strcmp(s1, s2) == 0 ) return 1; /* * fuzzy match: NOTE this will break if s1 > s2 */ if ( strchr(s1, '?') != NULL || strchr(s1, '*') != NULL ) { /* * look over the minimal map between the two strings */ for (i=0; i < strlen(s1) && i < strlen(s2); i++) { if ( s1[i] != s2[i] ) { switch( s1[i] ) { case '*': /* wild card, multiple characters */ return 1; case '?': /* wild card, single character */ if ( isspace(s2[i]) ) return 0; break; default: /* mismatch */ return 0; } } } return 1; } return 0; } void sortIndex(double *values, int *index, int size) { #ifdef SHELL_SORT double *sortValue; double approx_ln_to_log2, ftmp; int logsplit, itmp, psort; int i, j, k; approx_ln_to_log2 = 1.0 / log( (double) 2.0 ) + 0.000001; sortValue = (double *) safe_malloc(sizeof(double) * size); for (i=0; i < size; i++) sortValue[i] = values[i]; logsplit = ( log( (double) size ) * approx_ln_to_log2 ); i = size; for (psort = 1; psort <= logsplit; psort++) { i >>= 1; for (j = i+1; j <= size; j++) { k = j - i; ftmp = sortValue[j-1]; itmp = index[j-1]; while (k >= 1 && sortValue[k-1] > ftmp ) { sortValue[k+i-1] = sortValue[k-1]; index[k+i-1] = index[k-1]; k -= i; } sortValue[k+i-1] = ftmp; index[k+i-1] = itmp; } } safe_free(sortValue); #else double *sortValue; double ftmp; int itmp; int i,j; sortValue = (double *) safe_malloc(sizeof(double) * size); for (i=0; i < size; i++) sortValue[i] = values[i]; for (i=1; i < size; i++) { ftmp = sortValue[i]; itmp = index[i]; j = i-1; while (j >= 0 && sortValue[j] > ftmp) { sortValue[j+1] = sortValue[j]; index[j+1] = index[j]; j--; } sortValue[j+1] = ftmp; index[j+1] = itmp; } safe_free(sortValue); #endif } void sortIndexFloat(float *values, int *index, int size) { #ifdef SHELL_SORT float *sortValue; float approx_ln_to_log2, ftmp; int logsplit, itmp, psort; int i, j, k; approx_ln_to_log2 = 1.0 / log( (double) 2.0 ) + 0.000001; sortValue = (float *) safe_malloc(sizeof(float) * size); for (i=0; i < size; i++) sortValue[i] = values[i]; logsplit = ( log( (double) size ) * approx_ln_to_log2 ); i = size; for (psort = 1; psort <= logsplit; psort++) { i >>= 1; for (j = i+1; j <= size; j++) { k = j - i; ftmp = sortValue[j-1]; itmp = index[j-1]; while (k >= 1 && sortValue[k-1] > ftmp ) { sortValue[k+i-1] = sortValue[k-1]; index[k+i-1] = index[k-1]; k -= i; } sortValue[k+i-1] = ftmp; index[k+i-1] = itmp; } } safe_free(sortValue); #else float *sortValue; float ftmp; int itmp; int i,j; sortValue = (float *) safe_malloc(sizeof(float) * size); for (i=0; i < size; i++) sortValue[i] = values[i]; for (i=1; i < size; i++) { ftmp = sortValue[i]; itmp = index[i]; j = i-1; while (j >= 0 && sortValue[j] > ftmp) { sortValue[j+1] = sortValue[j]; index[j+1] = index[j]; j--; } sortValue[j+1] = ftmp; index[j+1] = itmp; } safe_free(sortValue); #endif } /* * The goal of this routine is to create a new ptrajState (newstate) * based on the old ptrajState (oldstate) deleting atoms that are * set in the mask array. This will modify the value of newstate. */ void modifyStateByMask(ptrajState **newstatep, ptrajState **oldstatep, int *mask, int strip) { int i, ires, isol, imol; int j, jres, jsol, jmol; int curres, cursol, curmol; int k; Name *atomName, *residueName; double *charges, *masses; int *startsol, *stopsol, *ipres, *moleculeInfo; ptrajState *newstate, *oldstate; oldstate = *oldstatep; /* * allocate space for the new state */ newstate = (ptrajState *) safe_malloc(sizeof(ptrajState)); INITIALIZE_ptrajState(newstate); /* * allocate space for temporary arrays and perform initialization */ atomName = (Name *) safe_malloc(sizeof(Name) * oldstate->atoms); charges = (double *) safe_malloc(sizeof(double) * oldstate->atoms); masses = (double *) safe_malloc(sizeof(double) * oldstate->atoms); residueName = (Name *) safe_malloc(sizeof(Name) * oldstate->residues); ipres = (int *) safe_malloc(sizeof(int) * (oldstate->residues+1)); if (oldstate->molecules) moleculeInfo = (int *) safe_malloc(sizeof(int) * oldstate->molecules); if (oldstate->solventMolecules) { startsol = (int *) safe_malloc(sizeof(int) * oldstate->solventMolecules); stopsol = (int *) safe_malloc(sizeof(int) * oldstate->solventMolecules); for (i=0; i < oldstate->solventMolecules; i++) startsol[i] = -1; } j = 0; jres = -1; jsol = -1; jmol = -1; ires = -1; isol = -1; imol = -1; /* * loop over all atoms and set up information for the newstate if the atom is * not to be deleted... */ for (i=0; i < oldstate->atoms; i++) { curres = atomToResidue(i+1, oldstate->residues, oldstate->ipres)-1; if ( mask[i] == (strip ? 0 : 1) ) { /* * this atom is not to be deleted */ /* * copy over atom information */ strcpy(atomName[j], oldstate->atomName[i]); charges[j] = oldstate->charges[i]; masses[j] = oldstate->masses[i]; /* * check to see if we are in the same residue or not * and copy relevant information */ if (ires == -1 || ires != curres) { jres++; strcpy(residueName[jres], oldstate->residueName[curres]); ipres[jres] = j+1; ires = curres; } /* * deal with the molecule information */ if (oldstate->molecules) { curmol = atomToMolecule(i+1, oldstate->molecules, oldstate->moleculeInfo)-1; if (imol == -1 || imol != curmol) { jmol++; moleculeInfo[jmol] = oldstate->moleculeInfo[curmol]; for (k=1; k < oldstate->moleculeInfo[curmol]; k++) if (mask[i+k] == (strip ? 1 : 0)) moleculeInfo[jmol]--; imol = curmol; } } /* * deal with the solvent information */ if (oldstate->solventMolecules) { cursol = atomToSolventMolecule(i+1, oldstate->solventMolecules, oldstate->solventMoleculeStart, oldstate->solventMoleculeStop)-1; if (cursol >= 0 && (isol == -1 || isol != cursol)) { jsol++; startsol[jsol] = j; stopsol[jsol] = j; for (k=i; k < oldstate->solventMoleculeStop[cursol]; k++) if (mask[k] == (strip ? 0 : 1)) stopsol[jsol]++; isol = cursol; } } /* * increment the new atom counter */ j++; } } /* * fix up IPRES */ ipres[jres+1] = j+1; /* * set up the newstate using the date placed into the temporary arrays; * free up unneeded space as we go... */ newstate->atoms = j; newstate->masses = (double *) safe_malloc(sizeof(double) * newstate->atoms); newstate->charges = (double *) safe_malloc(sizeof(double) * newstate->atoms); newstate->atomName = (Name *) safe_malloc(sizeof(Name) * newstate->atoms); for (i=0; i < newstate->atoms; i++) { newstate->masses[i] = masses[i]; newstate->charges[i] = charges[i]; strcpy(newstate->atomName[i], atomName[i]); } safe_free(masses); safe_free(charges); safe_free(atomName); newstate->residues = jres+1; newstate->ipres = (int *) safe_malloc(sizeof(int) * (newstate->residues+1)); newstate->residueName = (Name *) safe_malloc(sizeof(Name) * newstate->residues); for (i=0; i < newstate->residues; i++) { newstate->ipres[i] = ipres[i]; strcpy(newstate->residueName[i], residueName[i]); } newstate->ipres[i] = ipres[i]; safe_free(ipres); safe_free(residueName); newstate->IFBOX = oldstate->IFBOX; if (jsol < 0) { jsol = 0; newstate->solventMolecules = jsol; newstate->solventMoleculeStart = NULL; newstate->solventMoleculeStop = NULL; newstate->solventMask = NULL; } else { newstate->solventMolecules = jsol+1; newstate->solventMoleculeStart = (int *) safe_malloc(sizeof(int) * (jsol+1)); newstate->solventMoleculeStop = (int *) safe_malloc(sizeof(int) * (jsol+1)); for (i=0; i < newstate->solventMolecules; i++) { newstate->solventMoleculeStart[i] = startsol[i]; newstate->solventMoleculeStop[i] = stopsol[i]; } newstate->solventMask = (int *) safe_malloc(sizeof(int) * newstate->atoms); newstate->solventAtoms = 0; cursol = 0; for (i=0; i < newstate->atoms; i++) { newstate->solventMask[i] = 0; if (i == newstate->solventMoleculeStart[cursol]) { for (j=i; j < newstate->solventMoleculeStop[cursol]; j++) { newstate->solventAtoms++; newstate->solventMask[j] = 1; } i = newstate->solventMoleculeStop[cursol]-1; cursol++; } } safe_free(startsol); safe_free(stopsol); } newstate->maxFrames = oldstate->maxFrames; if (oldstate->molecules) { newstate->molecules = jmol+1; newstate->moleculeInfo = (int *) safe_malloc(sizeof(int) * newstate->molecules); for (i=0; imolecules; i++) newstate->moleculeInfo[i] = moleculeInfo[i]; safe_free(moleculeInfo); } for (i=0; i<6; i++) newstate->box[i] = oldstate->box[i]; *newstatep = newstate; } void printAtomMask(int *mask, ptrajState *state) { printAtomMaskDetailed(mask, state->atoms, state->residues, state->atomName, state->residueName, state->ipres); } void printHBondMask(int *mask, int *maskH1, int *maskH2, int *maskH3, ptrajState *state) { int i; fprintf(stdout, " Atom# Residue# Name -- Atom# Residue# Name\n"); for (i=0; i < state->atoms; i++) { if (mask[i] != 0) { fprintf(stdout, " %5i %4i %s ", i+1, atomToResidue(i+1, state->residues, state->ipres), state->atomName[i]); if (maskH1 != NULL) { if (maskH1[i] >= 0) { fprintf(stdout, " -- %5i %4i %s ", maskH1[i]+1, atomToResidue(maskH1[i]+1, state->residues, state->ipres), state->atomName[maskH1[i]]); } if (maskH2[i] >= 0) { fprintf(stdout, "\n %5i %4i %s ", i+1, atomToResidue(i+1, state->residues, state->ipres), state->atomName[i]); fprintf(stdout, " -- %5i %4i %s ", maskH2[i]+1, atomToResidue(maskH2[i]+1, state->residues, state->ipres), state->atomName[maskH2[i]]); } if (maskH3[i] >= 0) { fprintf(stdout, "\n %5i %4i %s ", i+1, atomToResidue(i+1, state->residues, state->ipres), state->atomName[i]); fprintf(stdout, " -- %5i %4i %s ", maskH3[i]+1, atomToResidue(maskH3[i]+1, state->residues, state->ipres), state->atomName[maskH3[i]]); } } fprintf(stdout, "\n"); } } } void printHBondInfo(int numdonor, int *donor, int numacceptor, int *acceptor, int *acceptorH, ptrajState *state) { int i; if (numdonor > 0 && donor != NULL) { fprintf(stdout, " HBOND DONOR LIST:\n"); fprintf(stdout, " Atom# Residue# Name\n"); for (i=0; i < numdonor; i++) { fprintf(stdout, " %5i %4i %s\n", donor[i]+1, atomToResidue(donor[i]+1, state->residues, state->ipres), state->atomName[donor[i]]); } } if (numacceptor > 0 && acceptor != NULL) { fprintf(stdout, " HBOND ACCEPTOR LIST:\n"); fprintf(stdout, " Atom# Residue# Name -- Atom# Residue# Name\n"); for (i=0; i < numacceptor; i++) { fprintf(stdout, " %5i %4i %s ", acceptor[i]+1, atomToResidue(acceptor[i]+1, state->residues, state->ipres), state->atomName[acceptor[i]]); fprintf(stdout, " -- %5i %4i %s\n", acceptorH[i]+1, atomToResidue(acceptorH[i]+1, state->residues, state->ipres), state->atomName[acceptorH[i]]); } } } int * processAtomMaskDetailed( char *maskString, int atoms, int residues, Name *atomName, Name *residueName, int *ipres) { int actualAtoms; int not = 0; char *maskp; char *tmp; int *residueMask; int residueMaskActive = 0; int *atomMask; int atomMaskActive = 0; int res1, res2; int atom1, atom2; int continuation = 0; int i, j; stackType *residueStack = NULL; stackType *atomStack = NULL; parseEntry *pp; Name name; maskp = maskString; skipWhitespace(maskp); /* * allocate mask strings */ atomMask = (int *) safe_malloc(sizeof(int) * atoms); residueMask = (int *) safe_malloc(sizeof(int) * residues); bzero(atomMask, sizeof(int) * atoms); bzero(residueMask, sizeof(int) * residues); /* * special case, choose ALL atoms */ if ( maskp[0] == (char) 0 || maskp[0] == '*' ) { for (i=0; i < atoms; i++) atomMask[i] = 1; goto clean_return; } /* * get rid of all NOT characters "~"; only one is significant * and set NOT status if one is found... */ while ( (tmp = strchr(maskString, '~' )) != NULL ) { not = 1; tmp[0] = ' '; } /* * check for error */ if (strchr(maskp, ':') == NULL && strchr(maskp, '@') == NULL) { fprintf(stdout, "WARNING: Error in mask string, no \"@\" or \":\" present (%s)\n", maskString); safe_free( atomMask ); safe_free( residueMask ); return NULL; } /* * do the main "parsing" */ skipWhitespace(maskp); while ( maskp[0] != (char) 0 ) { if ( maskp[0] == ':' ) { maskp++; for (;;) { skipWhitespace(maskp); pp = parseToken(&maskp, PARSE_RESIDUE); pushBottomStack(&residueStack, (void *) pp); if (maskp[0] == (char) 0 || maskp[0] == '@' || maskp[0] == ':') break; } } if ( maskp[0] == '@' ) { maskp++; for (;;) { skipWhitespace(maskp); pp = parseToken(&maskp, PARSE_ATOM); pushBottomStack(&atomStack, (void *) pp); if ( maskp[0] == (char) 0 || maskp[0] == ':' ) break; if ( maskp[0] == '@' ) maskp++; } } /* * now process the atomStack and residueStack */ if ( not ) { for (i=0; i < atoms; i++) atomMask[i] = 1; for (i=0; i < residues; i++) residueMask[i] = 1; } while ( residueStack != NULL ) { if ( continuation ) { res1 = res2; res2 = -1; } pp = (parseEntry *) popStack( &residueStack ); if ( pp->isnumber ) { if ( sscanf(pp->token, "%i", &res2) != 1 ) { fprintf(stdout, "WARNING: error parsing atom mask\n"); safe_free( atomMask ); safe_free( residueMask ); return NULL; } res2--; if (continuation) { continuation = 0; if (res1 < 0) res1 = 0; if (res2 >= residues) res2 = residues-1; if (res2 < res1) res2 = res1; for (i = res1; i <= res2; i++) { residueMask[i] = (not ? 0 : 1); residueMaskActive = 1; } } else { residueMask[res2] = (not ? 0 : 1); residueMaskActive = 1; } if ( pp->operator == PARSE_CONTINUATION ) continuation = 1; } else { continuation = 0; strcpy(name, " "); for (i=0; i < strlen(pp->token) && i < NAME_SIZE; i++) { name[i] = pp->token[i]; } name[NAME_SIZE-1] = (char) 0; for (i=0; i < residues; i++) if ( isMatch(name, residueName[i]) ) { residueMask[i] = (not ? 0 : 1); residueMaskActive = 1; } } safe_free(pp->token); pp->token=NULL; safe_free(pp); } while ( atomStack != NULL ) { if ( continuation ) { atom1 = atom2; atom2 = -1; } pp = (parseEntry *) popStack( &atomStack ); if ( pp->isnumber ) { if ( sscanf(pp->token, "%i", &atom2) != 1 ) { fprintf(stdout, "WARNING: error parsing atom mask\n"); safe_free( atomMask ); safe_free( residueMask ); return NULL; } atom2--; if (continuation) { continuation = 0; if (atom1 < 0) atom1 = 0; if (atom2 > atoms) atom2 = atoms; if (atom2 < atom1) atom2 = atom1; for (i = atom1; i < atom2; i++) { atomMask[i] = (not ? 0 : 1); atomMaskActive = 1; } } else { atomMask[atom2] = (not ? 0 : 1); atomMaskActive = 1; } if ( pp->operator == PARSE_CONTINUATION ) continuation = 1; } else { continuation = 0; strcpy(name, " "); for (i=0; i < strlen(pp->token) && i < NAME_SIZE; i++) { name[i] = pp->token[i]; } name[NAME_SIZE-1] = (char) 0; for (i=0; i < atoms; i++) if ( isMatch(name, atomName[i]) ) { atomMask[i] = (not ? 0 : 1); atomMaskActive = 1; } } safe_free(pp->token); pp->token = NULL; safe_free(pp); } if ( atomMaskActive && residueMaskActive ) { for (i=0; i < residues; i++) { if ( residueMask[i] == 0 ) { for (j = ipres[i]-1; j < ipres[i+1]-1; j++) { atomMask[j] = 0; } } } } else if ( residueMaskActive ) { for (i=0; i < residues; i++) { for (j = ipres[i]-1; j < ipres[i+1] - 1; j++) { if ( residueMask[i] ) atomMask[j] = 1; else if (not) atomMask[j] = 0; } } } atomMaskActive = 0; residueMaskActive = 0; } clean_return: actualAtoms = 0; for (i=0; i < atoms; i++) if ( atomMask[i] ) actualAtoms++; if ( tmp = strchr(maskString, '\n') ) tmp[0] = (char) 0; if (actualAtoms > 0) { fprintf(stdout, "Mask %s%s] represents %i atoms\n", (not ? "[~" : "["), maskString, actualAtoms); } else { fprintf(stdout, "Mask %s%s] represents %i atoms ", (not ? "[~" : "["), maskString, actualAtoms); fprintf(stdout, "!!!NO ATOMS DETECTED!!!\n"); safe_free(atomMask); atomMask = NULL; } safe_free( residueMask ); return ( atomMask ); } int * processAtomMask( char *maskString, ptrajState *state ) { return (processAtomMaskDetailed(maskString, state->atoms, state->residues, state->atomName, state->residueName, state->ipres)); } int * processAtomMaskWrapper(char *buffer, int printit, int returnit) { int *mask; int i; ptrajState *state; state = (ptrajState *) safe_malloc( sizeof(ptrajState) ); INITIALIZE_ptrajState(state); state->atoms = parm->NTOTAT; state->atomName = (Name *) safe_malloc( sizeof(Name) * state->atoms); state->masses = (double *) safe_malloc(sizeof(double) * state->atoms); state->charges = NULL; for (i=0; i < state->atoms; i++) { strcpy(state->atomName[i], parm->atom[i].igraph); state->masses[i] = parm->atom[i].amass; } state->residues = parm->NTOTRS; state->residueName = (Name *) safe_malloc (sizeof(Name) * state->residues); for (i=0; i < state->residues; i++) { strcpy(state->residueName[i], parm->residue[i].labres); } state->ipres = (int *) safe_malloc(sizeof(int) * (state->residues+1)); for (i=0; i <= state->residues; i++) { state->ipres[i] = parm->residue[i].ipres; } mask = processAtomMaskDetailed(buffer, state->atoms, state->residues, state->atomName, state->residueName, state->ipres); if (printit) { printAtomMaskDetailed(mask, state->atoms, state->residues, state->atomName, state->residueName, state->ipres); } if ( state != NULL ) { safe_free( state->atomName ); safe_free( state->residueName ); safe_free( state->ipres ); safe_free( state->masses ); safe_free( state ); } if (returnit) { return(mask); } else { safe_free(mask); return(NULL); } } void checkAtomMask(char *buffer) { processAtomMaskWrapper(buffer, 1, 0); } int * returnAtomMask(char *buffer) { return( processAtomMaskWrapper(buffer, 0, 1) ); } void atomMaskIsActive(int *mask, ptrajState *state, int *activep, int *firstp) { int i, active, first; if (mask == NULL) { *activep = 0; return; } active = 0; for (i=0; i < state->atoms; i++) { if (mask[i] != 0) { active++; if (active == 1) first = i; } } *activep = active; *firstp = first; } scalarInfo * scalarStackGetName(stackType **scalarStackP, char *name) { stackType *s; scalarInfo *info, *match; match = NULL; for (s = *scalarStackP; s != NULL; s = s->next) { info = (scalarInfo *) s->entry; if ( strcmp(info->name, name) == 0 ) match = info; } return (match); } transformHBondInfo * hbondInfoStackGetName(stackType **scalarStackP, char *name) { stackType *s; transformHBondInfo *info, *match; match = NULL; for (s = *scalarStackP; s != NULL; s = s->next) { info = (transformHBondInfo *) s->entry; if ( strcmp(info->name, name) == 0 ) match = info; } return (match); } void boxToRecip(double box[6], double ucell[9], double recip[9]) { double u12x,u12y,u12z; double u23x,u23y,u23z; double u31x,u31y,u31z; double volume; ucell[0] = box[0]; ucell[1] = 0.0; ucell[2] = 0.0; ucell[3] = box[1]*cos(DEGRAD*box[5]); ucell[4] = box[1]*sin(DEGRAD*box[5]); ucell[5] = 0.0; ucell[6] = box[2]*cos(DEGRAD*box[4]); ucell[7] = (box[1]*box[2]*cos(DEGRAD*box[3]) - ucell[6]*ucell[3]) / ucell[4]; ucell[8] = sqrt(box[2]*box[2] - ucell[6]*ucell[6] - ucell[7]*ucell[7]); u23x = ucell[4]*ucell[8] - ucell[5]*ucell[7]; u23y = ucell[5]*ucell[6] - ucell[3]*ucell[8]; u23z = ucell[3]*ucell[7] - ucell[4]*ucell[6]; u31x = ucell[7]*ucell[2] - ucell[8]*ucell[1]; u31y = ucell[8]*ucell[0] - ucell[6]*ucell[2]; u31z = ucell[6]*ucell[1] - ucell[7]*ucell[0]; u12x = ucell[1]*ucell[5] - ucell[2]*ucell[4]; u12y = ucell[2]*ucell[3] - ucell[0]*ucell[5]; u12z = ucell[0]*ucell[4] - ucell[1]*ucell[3]; volume=ucell[0]*u23x + ucell[1]*u23y + ucell[2]*u23z; recip[0] = u23x/volume; recip[1] = u23y/volume; recip[2] = u23z/volume; recip[3] = u31x/volume; recip[4] = u31y/volume; recip[5] = u31z/volume; recip[6] = u12x/volume; recip[7] = u12y/volume; recip[8] = u12z/volume; } /* * Calculate the minimum possible distance between periodic images. * This routine assumes that the ucell and recip information have * previously been set via a call to boxToRecip */ double calculateMinImagedDistance2(double *box, double *ucell, double *recip, double x1, double y1, double z1, double x2, double y2, double z2, int *rix, int *riy, int *riz) { double X, Y, Z; double min, dist; double fx, fy, fz; int ix, iy, iz, i; if (box == NULL) return(0.0); min = 100.0 * (box[0]*box[0]+box[1]*box[1]+box[2]*box[2]); if (prnlev > 6) { fprintf(stdout, "ATOM 0 XXX A1 1 %7.3f %7.3f %7.3f\n", x1, y1, z1); fprintf(stdout, "ATOM 1 XXX A2 1 %7.3f %7.3f %7.3f\n", x2, y2, z2); } fx = x1*recip[0] + y1*recip[1] + z1*recip[2]; fy = x1*recip[3] + y1*recip[4] + z1*recip[5]; fz = x1*recip[6] + y1*recip[7] + z1*recip[8]; if (prnlev > 6) { i = 2; for (ix = -1; ix <= 1; ix++) { for (iy = -1; iy <= 1; iy++) { for (iz = -1; iz <= 1; iz++) { X = (fx+ix)*ucell[0] + (fy+iy)*ucell[3] + (fz+iz)*ucell[6]; Y = (fx+ix)*ucell[1] + (fy+iy)*ucell[4] + (fz+iz)*ucell[7]; Z = (fx+ix)*ucell[2] + (fy+iy)*ucell[5] + (fz+iz)*ucell[8]; fprintf(stdout, "ATOM %3i XXX B%-2i 1 %7.3f %7.3f %7.3f\n", i, i, X, Y, Z); i++; } } } } i = 1; for (ix = -1; ix <= 1; ix++) { for (iy = -1; iy <= 1; iy++) { for (iz = -1; iz <= 1; iz++) { X = ((fx+ix)*ucell[0] + (fy+iy)*ucell[3] + (fz+iz)*ucell[6]) - x2; Y = ((fx+ix)*ucell[1] + (fy+iy)*ucell[4] + (fz+iz)*ucell[7]) - y2; Z = ((fx+ix)*ucell[2] + (fy+iy)*ucell[5] + (fz+iz)*ucell[8]) - z2; dist = X*X + Y*Y + Z*Z; if (prnlev > 6) { fprintf(stdout, " IMAGE FAMILIAR distance %3i: %6.3f (%5i %5i %5i)\n", i++, dist, ix, iy, iz); } if (dist < min) { min = dist; *rix = ix; *riy = iy; *riz = iz; } } } } if (prnlev > 4) { fprintf(stdout, " IMAGE FAMILIAR, min distance is %6.3f (%5i %5i %5i)\n", min, *rix, *riy, *riz); } return(min); } /* * This routine will calculate a distance squared performing imaging, * such that the minimum imaged distance (squared) is returned. If * the box is orthorhomic, simply subtract off multiple of the box lengthes. * For non-orthorhomic, this is a little more tricky since this procedure * is no longer applicable. In this case, the ucell, recip and closest2 * information is used. "closest2" represents a cutoff distance, such that if * the calculate distance**2 is less than this return it (without calculating * all possible image distances). All possible images in each direction are * investigated currently. This should be made "smarter" to only search in the * appropriate octant based on the angle values. */ double calculateDistance2(int i, int j, double *x, double *y, double *z, double *box, double *ucell, double *recip, double closest2) /* i -- index for first atom j -- index for second atom x -- coordinate arrays y z box -- box coordinates ucell recip closest2 -- minimum distance between a periodic image */ { double X, Y, Z, W; double fx, fy, fz, min, dist; int ix, iy, iz; X = x[i] - x[j]; Y = y[i] - y[j]; Z = z[i] - z[j]; if (box == NULL || box[0] == 0.0) return(X*X + Y*Y + Z*Z); if (box[3] == 90.0 && box[4] == 90.0 && box[5] == 90.0) { /* * DO ORTHORHOMIBIC IMAGING (this is faster!) */ /* * rid sign information */ if ( X < 0 ) X = -X; if ( Y < 0 ) Y = -Y; if ( Z < 0 ) Z = -Z; /* * rid multiples of the box length and image */ while ( X > box[0] ) X = X - box[0]; while ( Y > box[1] ) Y = Y - box[1]; while ( Z > box[2] ) Z = Z - box[2]; /* * find shortest distance in the periodic reference */ W = box[0] - X; if ( W < X ) X = W; W = box[1] - Y; if ( W < Y ) Y = W; W = box[2] - Z; if ( W < Z ) Z = W; return ( X*X + Y*Y + Z*Z ); } else { /* * NON-ORTHORHOMBIC CASE: find shortest distance in periodic reference * This is a brute force check requiring up to 26 distance evaluations. * It has been adapted to be smarter by returning the first distance that * is shorter than the minimum possible distance between images. */ fx = x[j]*recip[0] + y[j]*recip[1] + z[j]*recip[2]; fy = x[j]*recip[3] + y[j]*recip[4] + z[j]*recip[5]; fz = x[j]*recip[6] + y[j]*recip[7] + z[j]*recip[8]; X = fx*ucell[0] + fy*ucell[3] + fz*ucell[6] - x[i]; Y = fx*ucell[1] + fy*ucell[4] + fz*ucell[7] - y[i]; Z = fx*ucell[2] + fy*ucell[5] + fz*ucell[8] - z[i]; min = X*X + Y*Y + Z*Z; if (prnlev > 3) { if (min > 0.0) printf("DISTANCE: NON IMAGED is %8.3f\n", sqrt(min)); printf("DISTANCE BOX: is %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n", box[0],box[1],box[2],box[3],box[4],box[5]); } if (closest2 != 0.0 && min < closest2) return (min); for (ix = -1; ix <= 1; ix++) { for (iy = -1; iy <= 1; iy++) { for (iz = -1; iz <= 1; iz++) { if (! (ix == 0 && iy == 0 && iz == 0) ) { X = ((fx+ix)*ucell[0] + (fy+iy)*ucell[3] + (fz+iz)*ucell[6]) - x[i]; Y = ((fx+ix)*ucell[1] + (fy+iy)*ucell[4] + (fz+iz)*ucell[7]) - y[i]; Z = ((fx+ix)*ucell[2] + (fy+iy)*ucell[5] + (fz+iz)*ucell[8]) - z[i]; dist = X*X + Y*Y + Z*Z; if (prnlev > 3) printf("DISTANCE: %3i %3i %3i is %8.3f\n", ix, iy, iz, sqrt(dist)); if (dist < min) { min = dist; if (closest2 != 0.0 && min < closest2) return(min); } } } } } return(min); } } ptrajState ** ptrajCurrentState() { return(&ptrajCurrentStatePointer); } ptrajState * ptrajCopyState(ptrajState **stateinp) { ptrajState *state, *statein; int i; /* * Make a copy of the current state */ statein = *stateinp; state = (ptrajState *) safe_malloc( sizeof(ptrajState) ); INITIALIZE_ptrajState(state); state->atoms = statein->atoms; state->atomName = (Name *) safe_malloc( sizeof(Name) * state->atoms); state->masses = (double *) safe_malloc(sizeof(double) * state->atoms); state->charges = (double *) safe_malloc(sizeof(double) * state->atoms); for (i=0; i < state->atoms; i++) { strcpy(state->atomName[i], statein->atomName[i]); state->masses[i] = statein->masses[i]; state->charges[i] = statein->charges[i]; } state->residues = statein->residues; state->residueName = (Name *) safe_malloc (sizeof(Name) * state->residues); for (i=0; i < state->residues; i++) { strcpy(state->residueName[i], statein->residueName[i]); } state->ipres = (int *) safe_malloc(sizeof(int) * (state->residues+1)); for (i=0; i <= state->residues; i++) { state->ipres[i] = statein->ipres[i]; } state->solventMolecules = statein->solventMolecules; state->solventAtoms = statein->solventAtoms; if (statein->solventMolecules) { state->solventMask = (int *) safe_malloc(sizeof(int) * statein->atoms); for (i=0; i < state->atoms; i++) state->solventMask[i] = statein->solventMask[i]; state->solventMoleculeStart = (int *) safe_malloc(sizeof(int) * statein->solventMolecules); state->solventMoleculeStop = (int *) safe_malloc(sizeof(int) * statein->solventMolecules); for (i=0; i < state->solventMolecules; i++) { state->solventMoleculeStart[i] = statein->solventMoleculeStart[i]; state->solventMoleculeStop[i] = statein->solventMoleculeStop[i]; } } else { state->solventMoleculeStart = NULL; state->solventMoleculeStop = NULL; state->solventMask = NULL; } state->IFBOX = statein->IFBOX; for (i=0; i < 6; i++) state->box[i] = statein->box[i]; if ( statein->molecules > 0 ) { state->molecules = statein->molecules; state->moleculeInfo = (int *) safe_malloc( sizeof(int) * state->molecules ); for (i=0; i < state->molecules; i++) { state->moleculeInfo[i] = statein->moleculeInfo[i]; } } return(state); } void ptrajPrintState(ptrajState *state) { int i,j; int curres; printf("Dumping state information...\n"); printf(" atoms: %i\n", state->atoms); printf(" residues: %i\n", state->residues); printf(" box length: %8.3f %8.3f %8.3f\n", state->box[0], state->box[1], state->box[2]); printf(" box angles: %8.3f %8.3f %8.3f\n", state->box[3], state->box[4], state->box[5]); printf(" molecules: %i\n", state->molecules); printf(" max frames: %i\n", state->maxFrames); if (prnlev > 1) { printf(" ATOM NAME RESIDUE NAME CHARGE MASS\n"); curres = 0; for (i=0; i < state->atoms; i++) { if (i == state->ipres[curres+1]-1) curres++; printf(" %8i %s %8i %s %8.3f %8.3f\n", i+1, state->atomName[i], curres+1, state->residueName[curres], state->charges[i], state->masses[i]); } printf(" Molecule information (atoms in each molecule):\n"); for (i=0; i < state->molecules; i++) { printf(" %5i", state->moleculeInfo[i]); if (i!=0 && i%10==0) printf("\n"); } printf("\n"); } if (state->solventMolecules > 0) { printf(" solvent molecules: %i (%i atoms)\n", state->solventMolecules, state->solventAtoms); printf(" solvent mask is: "); printAtomMask(state->solventMask, state); if (prnlev > 1) { for (i=0; i < state->solventMolecules; i++) { printf(" SOLVENT %4i: ", i+1); for (j=state->solventMoleculeStart[i]; j < state->solventMoleculeStop[i]; j++) { printf(" (%s %5i)", state->atomName[j], j+1); } printf("\n"); } } } } void ptrajClearState(ptrajState **stateinp) { ptrajState *state; state = *stateinp; if (state != NULL) { safe_free( state->atomName ); safe_free( state->residueName ); safe_free( state->ipres ); safe_free( state->masses ); safe_free( state->charges ); if (state->solventMolecules) { safe_free( state->solventMoleculeStart ); safe_free( state->solventMoleculeStop ); safe_free( state->solventMask ); } safe_free( state->moleculeInfo ); INITIALIZE_ptrajState(state); safe_free( state ); } } /* * checkCoordinates(): a routine to determine what type of coordinate * file specified by "filename" is and how many coordinate frames the * file represents. The file is opened up, checked, then the file * is closed. It is the responsibility of later routines to reopen * the file. [This requirement is so the file limit will not be blown * in the case of a user processing a boatload of files; in this way, * sequential file access is maintained.] * * On success, a "coordinateInfo *" will be returned with the filename, * format and start/stop represented. * * On failure, NULL is returned. */ coordinateInfo * checkCoordinates(char *filename, int totalAtoms) { FILE *fp; pdb_record r; coordType type = COORD_UNKNOWN; char buffer1[BUFFER_SIZE], buffer2[BUFFER_SIZE]; float junk1, junk2, junk3, junk4, junk5, junk6, junk7, junk8, junk9; float *binposScratch; int i, actualAtoms, binposeof; int isBox = 0; int isVelocity = 0; int lines_per_set; int start = 1; int stop = 1; double *x = NULL; Restart *restrt; coordinateInfo *returnValue; charmmTrajectoryInfo **charmmTrajectoryp = NULL; charmmTrajectoryInfo *charmmTrajectory = NULL; /* * Open up the file specified by "filename" */ if ( openFile(&fp, filename, "r") == 0 ) { fprintf(stdout, "WARNING in checkCoordinates(): Could not open file (%s)\n", filename); return NULL; } /* * Try to determine what "type" of file the input file is... * NOTE: if new coordType's are implemented, this section may * need to be modified. * * (a) read two lines of text * (b) check for pdb [this requires reading both lines of * text since the comment from an AMBER restrt or * trajectory may have text which "fools" the pdb routines...] * (c) look for the BINPOS magic string * (d) look for CHARMM trajectory magic string * (e) if not a pdb or BINPOS, check the second line to see if it * has three numbers which implies it is an AMBER trajectory file. */ if (fgets(buffer1, BUFFER_SIZE, fp) == NULL || fgets(buffer2, BUFFER_SIZE, fp) == NULL ) { fprintf(stdout, "WARNING in checkCoordinates(): EOF encountered prematurely (%s)\n", filename); return NULL; } /* * Is this file a PDB file? */ r = pdb_read_string(buffer1); if (r.record_type != PDB_UNKNOWN) { r = pdb_read_string(buffer2); if (r.record_type != PDB_UNKNOWN) type = COORD_PDB; } else if (strncmp(buffer1, "REMARK", 6) == 0) /* * this is a hack to allow CHARMM pdb files to be read or other PDB's * which have malformed remarks */ type = COORD_PDB; /* * Does this file have the magic header for the Scripps binary format? */ if (buffer1[0] == 'f' && buffer1[1] == 'x' && buffer1[2] == 'y' && buffer1[3] == 'z') { type = COORD_BINPOS; } /* * Does file file have the magic header for a CHARMM trajectory??? */ if (buffer1[4] == 'C' && buffer1[5] == 'O' && buffer1[6] == 'R' && buffer1[7] == 'D') { type = COORD_CHARMM_TRAJECTORY; } /* * Can we scan in three numbers from the first line, if so it is * an AMBER trajectory otherwise assume it is an AMBER restrt file */ if ( type == COORD_UNKNOWN ) { if ( (i = sscanf(buffer2, "%f%f%f", &junk1, &junk2, &junk3)) == 3 ) type = COORD_AMBER_TRAJECTORY; else type = COORD_AMBER_RESTART; } if ( type == COORD_UNKNOWN ) return NULL; /* * REOPEN THE FILE */ if ( (fp = safe_freopen(fp)) == NULL ) return NULL; /* * Read through the file, depending on type, and make sure * that they have the same number of atoms and casually check * to make sure the file isn't corrupted. In the case of a * trajectory, check to see how many frames are present in the * file... */ actualAtoms = 0; switch ( type ) { case COORD_PDB: /* * read PDB records until the EOF is encountered counting * the number of atoms encountered along the way... * NOTE: this code will NOT properly handle pdbs that are * strung together into a single file. */ for (;;) { r = pdb_read_record(fp); if ( r.record_type == PDB_ATOM || r.record_type == PDB_HETATM ) actualAtoms++; if ( r.record_type == PDB_END ) break; } break; case COORD_AMBER_RESTART: restrt = readAmberRestart(totalAtoms, filename, fp); if (restrt == NULL) return NULL; actualAtoms = restrt->natoms; isBox = restrt->isbox; isVelocity = restrt->restart; break; case COORD_AMBER_TRAJECTORY: /* * there is no perfect, easy, way to make sure that the * trajectory file has the "correct" number of atoms * in a general way. If the number of atoms * 3 is not * a multiple of 10, it is possible; however we do not check * this currently. */ actualAtoms = totalAtoms; if (fgets(buffer1, BUFFER_SIZE, fp) == NULL) { fprintf(stdout, "WARNING in checkCoordinates(): EOF encountered during reading of\n"); fprintf(stdout, " title from (%s)\n", filename); return NULL; } lines_per_set = (int) (actualAtoms * 3) / 10; if ( (actualAtoms * 3) % 10 ) lines_per_set += 1; for (i=0; fgets(buffer1, BUFFER_SIZE, fp) != NULL; i++) { if ( strchr(buffer1, (int) '*') != NULL ) { fprintf(stdout, "WARNING in checkCoordinates(): File is corrupted...\n"); fprintf(stdout, " '*' detected in AMBER trajectory (%s)\n", filename); break; } if (i == lines_per_set) { if (actualAtoms > 2 && sscanf(buffer1, "%f%f%f%f%f%f%f%f%f", &junk1, &junk2, &junk3, &junk4, &junk5, &junk6, &junk7, &junk8, &junk9) < 9) isBox = 1; } } stop = i / (lines_per_set + isBox); break; case COORD_BINPOS: if (openbinpos(fp) < 0) break; binposScratch = safe_malloc(sizeof(float) * totalAtoms * 3); binposeof = 0; for (i=0; binposeof == 0; i++) { if (readbinpos(fp, &actualAtoms, binposScratch, &binposeof) < 0) binposeof = 1; } stop = i-1; safe_free(binposScratch); break; case COORD_CHARMM_TRAJECTORY: /* * pre-process */ stop = -1; x = NULL; charmmTrajectoryp = (charmmTrajectoryInfo **) safe_malloc(sizeof(charmmTrajectoryInfo *)); *charmmTrajectoryp = NULL; readCharmmTrajectory(fp, charmmTrajectoryp, x, x, x, x, stop); charmmTrajectory = *charmmTrajectoryp; x = (double *) safe_malloc(sizeof(double) * charmmTrajectory->natrec); /* * read in sets until there are no more */ stop = 0; while ( readCharmmTrajectory(fp, charmmTrajectoryp, x, x, x, x, stop+1) ) stop++; actualAtoms = charmmTrajectory->natrec; isVelocity = 0; isBox = charmmTrajectory->icntrl[10]; if (charmmTrajectory->icntrl[0] != stop) { fprintf(stdout, "NOTE: this charmm trajectory contains %i sets, expecting %i\n", stop, charmmTrajectory->icntrl[0]); } safe_free(x); x = NULL; } if (actualAtoms != totalAtoms) { fprintf(stdout, "WARNING in checkCoordinates(): The actual number of atoms (%d)\n", actualAtoms ); fprintf(stdout, " does not match the expected number of atoms (%d) in (%s)\n", totalAtoms, filename); } returnValue = (coordinateInfo *) safe_malloc(sizeof(coordinateInfo)); INITIALIZE_coordinateInfo(returnValue); returnValue->filename = (char *) safe_malloc( sizeof(char) * (strlen(filename)+1)); strcpy(returnValue->filename, filename); returnValue->start = start; returnValue->stop = stop; returnValue->offset = 1; returnValue->isBox = isBox; returnValue->isVelocity = isVelocity; returnValue->type = type; if (type == COORD_CHARMM_TRAJECTORY) { returnValue->info = (void *) charmmTrajectory; } safe_fclose(fp); return ( returnValue ); } void checkCoordinatesWrap(char *filename) { coordinateInfo *info; charmmTrajectoryInfo *charmmTraj; byte u; info = checkCoordinates(filename, parm->NTOTAT); if (info == NULL) { fprintf(stdout, "WARNING in checkCoordinatesWrap. Encountered a problem analyzing\n"); fprintf(stdout, "coordinates from file %s\n", filename); return; } switch ( info->type ) { case COORD_PDB: fprintf(stdout, "File (%s) is a PDB file\n", info->filename); break; case COORD_AMBER_TRAJECTORY: fprintf(stdout, "File (%s) is an AMBER trajectory%s", info->filename, (info->isBox ? " with box coordinates" : "")); if (info->stop > 0) fprintf(stdout, " representing %i sets\n", info->stop); else fprintf(stdout, "\n"); break; case COORD_CHARMM_TRAJECTORY: charmmTraj = (charmmTrajectoryInfo *) info->info; fprintf(stdout, "File (%s) is a CHARMM trajectory in %s endian binary format %s", info->filename, (charmmTraj->byteorder ? "little" : "big"), (info->isBox ? "with box coordinates" : "")); if (info->stop > 0) fprintf(stdout, "representing %i sets\n", info->stop); else fprintf(stdout, "\n"); if (prnlev > 2) { printf(" NFILE = %i\n", charmmTraj->icntrl[0]); /* number of coordinate sets in file */ printf(" ISTEP = %i\n", charmmTraj->icntrl[1]); /* number of previous dynamics steps */ printf(" NINTV = %i\n", charmmTraj->icntrl[2]); /* frequency for saving coordinates */ printf(" NSAVC = %i\n", charmmTraj->icntrl[3]); /* number of steps for creation runs */ printf(" NSAVV = %i\n", charmmTraj->icntrl[4]); printf(" NDEGF = %i\n", charmmTraj->icntrl[7]); printf(" NFREA = %i\n", charmmTraj->icntrl[8]); u.i = charmmTraj->icntrl[9]; printf(" DELTA = %f\n", u.f); printf(" QCRYS = %i\n", charmmTraj->icntrl[10]); printf(" QDIM4 = %i\n", charmmTraj->icntrl[11]); printf(" VERNU = %i\n", charmmTraj->icntrl[19]); printf(" Dumping the title:\n"); printStack(&charmmTraj->titleStack, printString, NULL); } break; case COORD_AMBER_RESTART: fprintf(stdout, "File (%s) is an AMBER restart file ", info->filename); if ( info->isBox && info->isVelocity ) fprintf(stdout, "with box and velocity information\n"); else if (info->isBox) fprintf(stdout, "with box information\n"); else if (info->isVelocity) fprintf(stdout, "with velocity information\n"); else fprintf(stdout, "\n"); break; } } ptrajState * loadCharmmPSF(FILE *fd, int skipheader) { char *bufferp, *buffer, *pbuffer; int natom; double *charges, *masses; float charge, mass; int *tmp_ipres, *tmp_molinfo; Name *atomName; Name *tmp_resName; Name segid, cursegid, vdwt; int atom; int residue; int fixed; int vdwp; int i,j,k,len; int curres, curresidx, curmol, numcurmol; int solventMolecules, *solventMoleculeStart, *solventMoleculeStop; int *solventMask, solventAtoms; int *tmp_solventMoleculeStart, *tmp_solventMoleculeStop; int ntitle; ptrajState *state; buffer = (char *) safe_malloc(sizeof(char) * BUFFER_SIZE); pbuffer = (char *) safe_malloc(sizeof(char) * BUFFER_SIZE); /* * Open up the PSF file */ if (fd == NULL) { error("loadCharmmPSF()", "Couldn't open file psf\n"); } /* * grab header and make sure it matches "PSF" */ if (skipheader != 1) { bufferp = fgets(buffer, BUFFER_SIZE, fd); if (bufferp == NULL || strncmp(buffer, "PSF", 3) != 0) { error("loadCharmmPSF()", "File is not a CHARMM PSF file!\n"); } } fprintf(stdout, "Reading in CHARMM PSF file\n"); /* * skip the next line, read the title and the blank line that follows */ bufferp = fgets(buffer, BUFFER_SIZE, fd); bufferp = fgets(buffer, BUFFER_SIZE, fd); if (strstr(bufferp, "!NTITLE")) { if ( sscanf(bufferp, "%i", &ntitle) != 1 ) ntitle = 0; } fprintf(stdout, "Reading in the title...\n\n"); if (ntitle) { for (i=0; i<=ntitle; i++) { bufferp = fgets(buffer, BUFFER_SIZE, fd); fprintf(stdout, "%s", buffer); } } else { bufferp = fgets(buffer, BUFFER_SIZE, fd); while (buffer[0] == '*') bufferp = fgets(buffer, BUFFER_SIZE, fd); fprintf(stdout, "%s", buffer); } /* * get the total number of atoms */ bufferp = fgets(buffer, BUFFER_SIZE, fd); if (sscanf(buffer, "%i", &natom) != 1) error("loadCharmmPSF()", "Scanning natoms\n"); fprintf(stdout, "Total number is atoms is %d\n", natom); /* * read in the atom segid, residue, name and charge/mass info */ charges = (double *) safe_malloc(sizeof(double) * natom); masses = (double *) safe_malloc(sizeof(double) * natom); tmp_ipres = (int *) safe_malloc(sizeof(double) * natom); tmp_molinfo = (int *) safe_malloc(sizeof(double) * natom); atomName = (Name *) safe_malloc(sizeof(Name) * natom); tmp_resName = (Name *) safe_malloc(sizeof(Name) * natom); tmp_solventMoleculeStart = (int *) safe_malloc(sizeof(int) * natom); tmp_solventMoleculeStop = (int *) safe_malloc(sizeof(int) * natom); solventMask = (int *) safe_malloc(sizeof(int) * natom); bzero(solventMask, sizeof(int) * natom); bzero(tmp_solventMoleculeStart, sizeof(int) * natom); bzero(tmp_solventMoleculeStop, sizeof(int) * natom); curres = 1; curresidx = 1; tmp_ipres[0] = 1; solventMolecules = 0; solventAtoms = 0; fprintf(stdout, "Reading in the atom information...\n"); for (i=0; i < natom; i++) { bufferp = fgets(buffer, BUFFER_SIZE, fd); /* * atom seg res# res atom vdwp charge mass fixed? * 1 WAT 1 WAT OH 2 -0.820000 15.9994 0 */ if (sscanf(buffer, "%8i %4s %4i %4s %4s %4i %f %f%8i", &atom, (char *) &segid, &residue, (char *) &tmp_resName[i], (char *) &atomName[i], &vdwp, &charge, &mass, &fixed) < 9) { /* * assume this is an XPLOR psf file w/o SEGID information... */ if (sscanf(buffer, "%8i %8i %4s %4s %4s %f %f%8i", &atom, &residue, (char *) &tmp_resName[i], (char *) &atomName[i], (char *) &vdwt, &charge, &mass, &fixed) < 8) { error("loadCharmmPSF()", "Scanning atom information\n"); } } atomName[i][NAME_SIZE-1] = (char) 0; tmp_resName[i][NAME_SIZE-1] = (char) 0; for (j=NAME_SIZE-2; atomName[i][j] == (char) 0; j--) { atomName[i][j] = ' '; } for (j=NAME_SIZE-2; tmp_resName[i][j] == (char) 0; j--) { tmp_resName[i][j] = ' '; } if (i==0) { strcpy(cursegid, segid); curmol = 0; numcurmol = 1; } else { if (strcmp(cursegid, segid) == 0) { numcurmol++; } else { tmp_molinfo[curmol] = numcurmol; curmol++; numcurmol = 1; strcpy(cursegid, segid); } } charges[i] = charge; masses[i] = mass; if (prnlev > 1) printf("Atom %4i (%s) seg (%4s) res %4i (%s) charge: %7.3f mass: %7.3f fixed? %s %i\n", atom, atomName[i], segid, residue, tmp_resName[i], charge, mass, (fixed ? "yes" : "no"), vdwp); if (residue != curres) { tmp_ipres[curresidx] = i+1; curresidx++; curres = residue; } } tmp_ipres[curresidx] = i+1; tmp_molinfo[curmol] = numcurmol; for (i=0; i < curresidx; i++) { j = tmp_ipres[i]-1; if ( strcmp(tmp_resName[j], "WAT ") == 0 || strcmp(tmp_resName[j], "IP3 ") == 0 || strcmp(tmp_resName[j], "TIP3") == 0) { tmp_solventMoleculeStart[solventMolecules] = j; tmp_solventMoleculeStop[solventMolecules] = j+3; solventMolecules++; } } solventMoleculeStart = (int *) safe_malloc(sizeof(int) * (solventMolecules+1)); solventMoleculeStop = (int *) safe_malloc(sizeof(int) * (solventMolecules+1)); for (i = 0; i < solventMolecules; i++) { solventMoleculeStart[i] = tmp_solventMoleculeStart[i]; solventMoleculeStop[i] = tmp_solventMoleculeStop[i]; for (j=solventMoleculeStart[i]; j < solventMoleculeStop[i]; j++) { solventMask[j] = 1; solventAtoms++; } } if (prnlev > 1) { printf("The total number of molecules is %i\n", curmol+1); printf("The total number of residues is %i\n", curresidx); } fprintf(stdout, "Dumping out residue names:\n"); j = 0; k = 0; len = 0; pbuffer[0] = (char) 0; for (i=0; i < curresidx; i++) { k = tmp_ipres[i]-1; sprintf(buffer+j, "%5s", tmp_resName[k]); j += 5; if ( i == curresidx - 1 || ( i != 0 && (i+1) % 10 == 0 ) ) { buffer[j] = (char) 0; if ( strcmp(pbuffer, buffer) == 0 ) { if ( ! len ) { fprintf(stdout, " ...\n"); len = 1; } } else { fprintf(stdout, "%s\n", buffer); len = 0; strcpy(pbuffer, buffer); } j = 0; } } state = (ptrajState *) safe_malloc(sizeof(ptrajState)); INITIALIZE_ptrajState(state); state->box[0] = 0.0; state->box[1] = 0.0; state->box[2] = 0.0; state->box[3] = 90.0; state->box[4] = 90.0; state->box[5] = 90.0; state->atoms = natom; state->residues = curresidx; state->charges = charges; state->masses = masses; state->ipres = (int *) safe_malloc(sizeof(int) * (state->residues+1)); state->residueName = (Name *) safe_malloc(sizeof(Name) * state->residues); for (i=0; i <= state->residues; i++) state->ipres[i] = tmp_ipres[i]; for (i=0; i < state->atoms; i++) { curresidx = atomToResidue(i+1, state->residues, state->ipres)-1; strcpy(state->residueName[curresidx], tmp_resName[i]); i = state->ipres[curresidx+1]-2; } state->IFBOX = 0; state->molecules = curmol+1; state->moleculeInfo = (int *) safe_malloc(sizeof(int) * state->molecules); for (i=0; i < state->molecules; i++) { state->moleculeInfo[i] = tmp_molinfo[i]; } state->solventMask = solventMask; state->solventAtoms = solventAtoms; state->solventMolecules = solventMolecules; state->solventMoleculeStart = solventMoleculeStart; state->solventMoleculeStop = solventMoleculeStop; state->atomName = atomName; safe_free(buffer); safe_free(tmp_ipres); safe_free(tmp_molinfo); safe_free(tmp_resName); safe_free(tmp_solventMoleculeStart); safe_free(tmp_solventMoleculeStop); if (prnlev >= 0) ptrajPrintState(state); return state; } void ptrajInitializeState(ptrajState **statep, char *filename) { FILE *fp; char *buffer; ptrajState *state; int i, j, k, *start, *stop; /* * open up the filename (if it exists, otherwise prompt the user) * and determine whether it is a AMBER prmtop or CHARMM PSF file */ if (filename == NULL) filename = promptToOpenFile(&fp, "", "r", "Input the name of an AMBER prmtop or CHARMM PSF: "); else if ( openFile(&fp, filename, "r") == 0 ) error("ptrajInitializeState()", "Attempting to open parameter/topology file %s", filename); buffer = (char *) safe_malloc(sizeof(char)* BUFFER_SIZE); if (fgets(buffer, BUFFER_SIZE, fp) == NULL) error("ptrajInitializeState()", "Attempting to read parameter/topology file"); if (strcmp(buffer, "PSF \n") == 0) state = loadCharmmPSF(fp, 1); else { safe_freopen(fp); if ( parm != NULL ) clearParm( parm ); parm = safe_malloc( sizeof( Parm ) ); initializeParm(parm); parm->filename = filename; parm->fp = fp; readParm(); state = (ptrajState *) safe_malloc( sizeof(ptrajState) ); INITIALIZE_ptrajState(state); state->atoms = parm->NTOTAT; state->atomName = (Name *) safe_malloc( sizeof(Name) * state->atoms); state->masses = (double *) safe_malloc(sizeof(double) * state->atoms); state->charges = (double *) safe_malloc(sizeof(double) * state->atoms); for (i=0; i < state->atoms; i++) { strcpy(state->atomName[i], parm->atom[i].igraph); state->masses[i] = parm->atom[i].amass; state->charges[i] = parm->atom[i].chrg / CHARGE_TO_KCALS; } state->residues = parm->NTOTRS; state->residueName = (Name *) safe_malloc (sizeof(Name) * state->residues); for (i=0; i < state->residues; i++) { strcpy(state->residueName[i], parm->residue[i].labres); } state->ipres = (int *) safe_malloc(sizeof(int) * (state->residues+1)); for (i=0; i <= state->residues; i++) { state->ipres[i] = parm->residue[i].ipres; } state->box[0] = 0.0; state->box[1] = 0.0; state->box[2] = 0.0; state->box[3] = 90.0; state->box[4] = 90.0; state->box[5] = 90.0; state->maxFrames = 0; state->IFBOX = parm->IFBOX; if ( state->IFBOX ) { state->molecules = parm->box->nspm; state->moleculeInfo = (int *) safe_malloc( sizeof(int) * state->molecules ); for (i=0; i < state->molecules; i++) { state->moleculeInfo[i] = parm->box->nsp[i]; } state->box[0] = parm->box->box[0]; state->box[1] = parm->box->box[1]; state->box[2] = parm->box->box[2]; if (parm->box->beta != 90.0) { fprintf(stdout, "\nInitializing state: detected PBC beta angle in prmtop that is not 90.0!\n"); if (parm->box->beta > 109.47 && parm->box->beta < 109.48) { state->box[3] = 2.0*acos(1.0/sqrt(3.0))*RADDEG; state->box[4] = state->box[3]; state->box[5] = state->box[3]; fprintf(stdout, "Assuming this box is a truncated octahedron, angle is %f\n", state->box[3]); } else if (parm->box->beta == 60.0) { fprintf(stdout, "Assuming this box is a rhombic dodecahedron, i.e. alpha=gamma=60.0, beta=90.0\n"); state->box[3] = 60.0; state->box[4] = 90.0; state->box[5] = 60.0; } } } state->solventMolecules = 0; state->solventAtoms = 0; state->solventMask = NULL; state->solventMoleculeStart = NULL; state->solventMoleculeStop = NULL; start = (int *) safe_malloc(sizeof(int) * state->atoms); stop = (int *) safe_malloc(sizeof(int) * state->atoms); state->solventMask = (int *) safe_malloc(sizeof(int) * state->atoms); for (i=0; i < state->atoms; i++) state->solventMask[i] = 0; if ( ! parm->IFBOX) { /* * treat all the molecules starting with parm->box->nspsol * as solvent molecules IF there is box information */ j = 0; for (i=0; i < state->molecules; i++) { if (i+1 >= parm->box->nspsol) { /* * add this molecule to the solvent list */ state->solventAtoms += state->moleculeInfo[i]; for (k = j; k < j+state->moleculeInfo[i]; k++) state->solventMask[k] = 1; start[state->solventMolecules] = j; stop[ state->solventMolecules] = j+state->moleculeInfo[i]; state->solventMolecules++; } j += state->moleculeInfo[i]; } } else { /* * treat all the residues named "WAT " as solvent molecules */ for (i=0; i < state->residues; i++) { if (strcmp("WAT ", state->residueName[i]) == 0) { /* * add this residue to the list of solvent */ j = state->ipres[i+1]-state->ipres[i]; state->solventAtoms += j; start[state->solventMolecules] = state->ipres[i]-1; stop[ state->solventMolecules] = state->ipres[i+1]-1; state->solventMolecules++; for (k=state->ipres[i]-1; k < state->ipres[i+1]-1; k++) state->solventMask[k] = 1; } } } /* state->solventMoleculeStart = (int *) safe_malloc(sizeof(int) * state->solventMolecules); state->solventMoleculeStop = (int *) safe_malloc(sizeof(int) * state->solventMolecules); for (i=0; i < state->solventMolecules; i++) { state->solventMoleculeStart[i] = start[i]; state->solventMoleculeStop[i] = stop[i]; } safe_free(stop); safe_free(start); */ state->solventMoleculeStart = start; state->solventMoleculeStop = stop; } *statep = state; } void printCoordinateInfo( void *entry ) { coordinateInfo *p; charmmTrajectoryInfo *charmmTraj; byte u; int i; p = (coordinateInfo *) entry; if (p == NULL) { fprintf(stdout, " NULL entry\n"); return; } switch ( p->type ) { case COORD_BINPOS: fprintf(stdout, " File (%s) is a BINPOS file ", p->filename); if (p->stop > 0) { i = (p->stop - p->start)/p->offset+1; if (i != p->stop) fprintf(stdout, " with %i sets (processing only %i)\n", p->stop, i); else fprintf(stdout, " with %i sets\n", i); } else { fprintf(stdout, "\n"); } break; case COORD_PDB: fprintf(stdout, " File (%s) is a PDB file", p->filename); if (p->option1 == 1) fprintf(stdout, ": AMBER charges and radii to occupancy and temp factor columns"); else if (p->option1 == 2) fprintf(stdout, ": AMBER charges and PARSE radii to occupancy and temp factor columns"); fprintf(stdout, "\n"); break; case COORD_AMBER_TRAJECTORY: fprintf(stdout, " File (%s) is an AMBER trajectory%s", p->filename, (p->isBox ? " (with box info)" : "")); if (p->stop > 0) { i = (p->stop - p->start)/p->offset+1; if (i != p->stop) fprintf(stdout, " with %i sets (processing only %i)\n", p->stop, i); else fprintf(stdout, " with %i sets\n", i); } else { fprintf(stdout, "\n"); } break; case COORD_CHARMM_TRAJECTORY: charmmTraj = (charmmTrajectoryInfo *) p->info; fprintf(stdout, " File (%s) is a CHARMM trajectory in %s endian binary format\n", p->filename, (charmmTraj->byteorder ? "little" : "big")); fprintf(stdout, "%s", (p->isBox ? " with box information " : " ")); if (p->stop > 0) { i = (p->stop - p->start)/p->offset+1; if (i != p->stop) fprintf(stdout, "representing %i sets (processing only %i)\n", p->stop, i); else fprintf(stdout, "representing %i sets\n", p->stop); } else fprintf(stdout, "\n"); /* * NFILE -- number of coordinate sets in file * ISTEP -- number of previous dynamics steps * NINTV -- frequency for saving coordinates * NSAVC -- number of steps for creation run */ printf(" NFILE = %8i ISTEP = %8i NINTV = %8i NSAVC = %8i NSAVV = %8i\n", charmmTraj->icntrl[0], charmmTraj->icntrl[1], charmmTraj->icntrl[2], charmmTraj->icntrl[3], charmmTraj->icntrl[4]); u.i = charmmTraj->icntrl[9]; printf(" NDEGF = %8i NFREA = %8i DELTA = %8.4f QCRYS = %8i QDIM4 = %8i\n", charmmTraj->icntrl[7], charmmTraj->icntrl[8], u.f, charmmTraj->icntrl[10], charmmTraj->icntrl[11]); printf(" VERNU = %8i NATREC= %8i NFREAT = %8i\n", charmmTraj->icntrl[19], charmmTraj->natrec, charmmTraj->nfreat); printf(" Dumping the title:\n"); printStack(&charmmTraj->titleStack, printString, NULL); break; case COORD_AMBER_RESTART: fprintf(stdout, " File (%s) is an AMBER restart file ", p->filename); if ( p->isBox && p->isVelocity ) fprintf(stdout, "with box and velocity information\n"); else if (p->isBox) fprintf(stdout, "with box information\n"); else if (p->isVelocity) fprintf(stdout, "with velocity information\n"); else fprintf(stdout, "\n"); break; } if (prnlev > 1) { if ( p->file == NULL ) { fprintf(stdout, " [the FILE is currently closed]\n"); } else if (p->file == stdin) { fprintf(stdout, " [the FILE is standard input]\n"); } else if (p->file == stdout) { fprintf(stdout, " [the FILE is standard output]\n"); } else if (p->file == stderr) { fprintf(stdout, " [the FILE is standard error]\n"); } } if ( p->mask ) fprintf(stdout, " [The file has an active atom mask]\n"); } void ptrajCleanup() { actionInformation *a; coordinateInfo *f; int *mask; /* * Clean up the INPUT files, transformFileStack */ while (transformFileStack != NULL && (f = (coordinateInfo *) popStack(&transformFileStack)) != NULL ) { safe_free(f->filename); safe_free(f->info); safe_free(f->mask); safe_free(f->x); safe_free(f->y); safe_free(f->z); INITIALIZE_coordinateInfo(f); safe_free(f); } transformFileStack = NULL; /* * Clean up the ACTIONS */ while (transformActionStack != NULL && (a = (actionInformation *) popStack(&transformActionStack)) != NULL) { /* * Free any associated mask */ if (a->mask) safe_free(a->mask); /* * Free any associated state information */ if (a->state != NULL) { ptrajClearState(&a->state); a->state = NULL; } /* * Clean up any of the complex arguments as necessary; this * is done by the associated action function in the PTRAJ_CLEANUP * mode */ if (a->type != TRANSFORM_TRANSFORM && a->type != TRANSFORM_NOOP) { if (a->fxn != NULL) { a->fxn(a, NULL, NULL, NULL, NULL, PTRAJ_CLEANUP); } INITIALIZE_actionInformation(a); } /* * Free up the action */ safe_free(a); } transformActionStack = NULL; /* * Free up outInfo */ if (outInfo != NULL) { safe_free(outInfo->filename); outInfo->filename = NULL; safe_free(outInfo->mask); outInfo->mask = NULL; safe_free(outInfo->info); outInfo->info = NULL; safe_free(outInfo); outInfo = NULL; } /* * Free up referenceInfo */ while (transformReferenceStack != NULL && (f = (coordinateInfo *) popStack(&transformReferenceStack)) != NULL ) { safe_free(f->filename); safe_free(f->info); safe_free(f->mask); safe_free(f->x); safe_free(f->y); safe_free(f->z); INITIALIZE_coordinateInfo(f); safe_free(f); } transformReferenceStack = NULL; referenceInfo = NULL; /* * Clean up hbond donor/acceptor information */ while (hbondDonorStack != NULL && (mask = (int *) popStack(&hbondDonorStack)) != NULL) { safe_free(mask); } while (hbondAcceptorStack != NULL && (mask = (int *) popStack(&hbondAcceptorStack)) != NULL) { safe_free(mask); } while (hbondAcceptorH1Stack != NULL && (mask = (int *) popStack(&hbondAcceptorH1Stack)) != NULL) { safe_free(mask); } while (hbondAcceptorH2Stack != NULL && (mask = (int *) popStack(&hbondAcceptorH2Stack)) != NULL) { safe_free(mask); } while (hbondAcceptorH3Stack != NULL && (mask = (int *) popStack(&hbondAcceptorH3Stack)) != NULL) { safe_free(mask); } } int ptrajPreprocessInputCoordinates(coordinateInfo *currentCoordinateInfo) { charmmTrajectoryInfo **charmmTrajp; charmmTrajectoryInfo *charmmTraj; char buffer[BUFFER_SIZE]; /* * open up the file */ if ( ! openFile(¤tCoordinateInfo->file, currentCoordinateInfo->filename, "r") ) { fprintf(stdout, "WARNING in ptrajPreprocessInputCoordinates(): Error on opening\n"); fprintf(stdout, "input coordinate file (%s)\n", currentCoordinateInfo->filename); return 1; } /* * preprocess the input coordinates as necessary * (i.e. to remove titles, etc.) */ switch(currentCoordinateInfo->type) { case COORD_PDB: fprintf(stdout, "\nProcessing PDB file %s\n", currentCoordinateInfo->filename); break; case COORD_BINPOS: fprintf(stdout, "\nProcessing BINPOS file %s\n", currentCoordinateInfo->filename); break; case COORD_AMBER_RESTART: fprintf(stdout, "\nProcessing AMBER restart file %s\n", currentCoordinateInfo->filename); break; case COORD_AMBER_TRAJECTORY: fprintf(stdout, "\nProcessing AMBER trajectory file %s\n", currentCoordinateInfo->filename); break; case COORD_CHARMM_TRAJECTORY: fprintf(stdout, "\nProcessing CHARMM trajectory file %s\n", currentCoordinateInfo->filename); break; } switch ( currentCoordinateInfo->type ) { case COORD_AMBER_TRAJECTORY: /* * we need to read the title line to set up for processing * by readAmberTrajectory() */ if ( fgets(buffer, BUFFER_SIZE, currentCoordinateInfo->file) == NULL ) { fprintf(stdout, "WARNING: Error on processing the title from the AMBER\n"); fprintf(stdout, "trajectory (%s)\n", currentCoordinateInfo->filename); return 1; } break; case COORD_BINPOS: if ( openbinpos( currentCoordinateInfo->file) < 0 ){ fprintf(stdout, "Error on opening BINPOS file %s\n", currentCoordinateInfo->filename); return 1; } break; case COORD_CHARMM_TRAJECTORY: /* * read in all the header information; this is done by calling with a * negative value for the current set... */ charmmTraj = (charmmTrajectoryInfo *) currentCoordinateInfo->info; charmmTrajp = (charmmTrajectoryInfo **) safe_malloc(sizeof(charmmTrajectoryInfo *)); *charmmTrajp = NULL; readCharmmTrajectory(currentCoordinateInfo->file, charmmTrajp, NULL, NULL, NULL, NULL, -1); charmmTraj = *charmmTrajp; break; case COORD_UNKNOWN: fprintf(stdout, "WARNING: Attempting to process a coordinate file of unknown type\n"); fprintf(stdout, "in ptraj(): ptrajPreprocessInputCoordinates(), ignoring...\n"); return 1; } return 0; } void ptrajProcessInputCoordinates(coordinateInfo *currentCoordinateInfo, ptrajState *state, double *X, double *Y, double *Z, double *box, int set, int *readCoordinates, int *processCoordinates) { pdb_record *pdb; float *binpos; int i,j,n_atoms,eof; charmmTrajectoryInfo **charmmTrajp; charmmTrajectoryInfo *charmmTraj; /* * NOTE: it is assumed that the box information is set to valid values on entry */ /* * READ IN THE CURRENT FILE, ONE SET AT A TIME */ switch( currentCoordinateInfo->type ) { case COORD_PDB: i = loadPdb(currentCoordinateInfo->file, &pdb); j = getCoordinatesFromPdb(pdb, X, Y, Z); if (i != state->atoms || j != state->atoms) { /* * on error, print warning but do not stop processing of this set. This * allows the specification of multiple reference sets */ fprintf(stdout, "\nWARNING in ptrajProcessInputCoordinates(): Unexpected number of atoms\n"); fprintf(stdout, " encountered when reading PDB!\n"); fprintf(stdout, " loadPdb returned %i\n", i); fprintf(stdout, " getCoordinatesFromPdb returned %i\n", j); fprintf(stdout, " actual number of atoms is %i\n", state->atoms); } /* * We assume that a PDB only contains ONE set of coordinates!!! */ safe_fclose(currentCoordinateInfo->file); currentCoordinateInfo->file = NULL; *readCoordinates = 0; *processCoordinates = 1; break; case COORD_AMBER_RESTART: if ( getCoordinatesFromRestart(currentCoordinateInfo->file, X, Y, Z, &box[0], &box[1], &box[2], &box[3], &box[4], &box[5]) != state->atoms ) { fprintf(stdout, "WARNING in ptrajProcessInputCoordinates(): Unexpected\n"); fprintf(stdout, "number of atoms encountered reading restrt file\n"); ptrajCleanup(); *readCoordinates = 0; *processCoordinates = 0; return; } /* * get default box information if necessary */ if ( state->IFBOX && (box[0] == 0.0 && box[1] == 0.0 && box[2] == 0.0) ) { box[0] = state->box[0]; box[1] = state->box[1]; box[2] = state->box[2]; } if ( state->IFBOX && (box[3] == 0.0 && box[4] == 0.0 && box[5] == 0.0) ) { box[3] = state->box[3]; box[4] = state->box[4]; box[5] = state->box[5]; } /* * AMBER restrt files only contain a single set of coordinates */ safe_fclose(currentCoordinateInfo->file); currentCoordinateInfo->file = NULL; *readCoordinates = 0; *processCoordinates = 1; break; case COORD_AMBER_TRAJECTORY: *readCoordinates = readAmberTrajectory(currentCoordinateInfo->file, state->atoms, X, Y, Z, box, set, currentCoordinateInfo->isBox); /* * get default box information if necessary */ if ( state->IFBOX && (box[0] == 0.0 && box[1] == 0.0 && box[2] == 0.0) ) { box[0] = state->box[0]; box[1] = state->box[1]; box[2] = state->box[2]; } if ( state->IFBOX && (box[3] == 0.0 && box[4] == 0.0 && box[5] == 0.0) ) { box[3] = state->box[3]; box[4] = state->box[4]; box[5] = state->box[5]; } /* * check to see if we've already loaded enough of the current * trajectory file */ if ( set > currentCoordinateInfo->stop ) *readCoordinates = 0; if ( *readCoordinates == 0 ) { *processCoordinates = 0; safe_fclose(currentCoordinateInfo->file); currentCoordinateInfo->file = NULL; } else { *processCoordinates = 1; } break; case COORD_CHARMM_TRAJECTORY: charmmTrajp = (charmmTrajectoryInfo ** ) safe_malloc(sizeof(charmmTrajectoryInfo *)); *charmmTrajp = (charmmTrajectoryInfo *) currentCoordinateInfo->info; *readCoordinates = readCharmmTrajectory(currentCoordinateInfo->file, charmmTrajp, X, Y, Z, box, set); break; case COORD_BINPOS: n_atoms = state->atoms; binpos = safe_malloc(sizeof(float) * n_atoms * 3); *readCoordinates = readbinpos( currentCoordinateInfo->file, &n_atoms, binpos, &eof); j=0; for (i=0; i currentCoordinateInfo->stop) ) { *readCoordinates = 0; } if ( *readCoordinates == 0 ) { *processCoordinates = 0; safe_fclose(currentCoordinateInfo->file); currentCoordinateInfo->file = NULL; } else { *processCoordinates = 1; } break; default: *readCoordinates = 0; *processCoordinates = 0; fprintf(stdout, "WARNING in ptrajProcessInputCoordinates(): Attempting to process\n"); fprintf(stdout, "a file of unknown type (%s)", currentCoordinateInfo->filename); fprintf(stdout, "Ignoring this file...\n"); } } void ptrajOutputCoordinates(ptrajState *state, int set, int append, int first, int atoms, double *X, double *Y, double *Z, double *box) { char buffer[BUFFER_SIZE]; pdb_record *pdb; switch( outInfo->type ) { case COORD_PDB: if (append) sprintf(buffer, "%s.%i", outInfo->filename, set); else sprintf(buffer, "%s", outInfo->filename); if ( openFile(&outInfo->file, buffer, "w") ) { pdb = ptrajStateToPdb(state, NULL, outInfo->option2); putCoordinatesInPdb(pdb, atoms, X, Y, Z); if (outInfo->option1 > 0) { /* * include charges and radii and dump out in higher precision */ putQandRInPdb(pdb, outInfo->option1, outInfo->option2, 0); savePdbHigh(outInfo->file, pdb); } else { savePdb(outInfo->file, pdb); } safe_fclose(outInfo->file); outInfo->file = NULL; safe_free( (void *) pdb ); } else { fprintf(stdout, "WARNING in ptrajOutputCoordinates(): Could not open output\n"); fprintf(stdout, "coordinate file %s, NOT dumping output to file.\n", buffer); } break; case COORD_AMBER_TRAJECTORY: /* * preprocessing to open file, if necessary */ if (first) { if ( ! openFile(&outInfo->file, outInfo->filename, "w") ) { fprintf(stdout, "WARNING in ptrajOutputCoordinates(): Error opening\n"); fprintf(stdout, "output coordinate file (%s)\n", outInfo->filename); } fprintf(outInfo->file, "rdparm transformed trajectory\n"); } if (outInfo->isBox) { dumpAmberTrajectory(outInfo->file, atoms, X, Y, Z, box); } else { dumpAmberTrajectory(outInfo->file, atoms, X, Y, Z, NULL); } break; case COORD_CHARMM_TRAJECTORY: /* * preprocessing to open file and write header */ if (first) { if ( ! openFile(&outInfo->file, outInfo->filename, "w") ) { fprintf(stdout, "WARNING in ptrajOutputCoordinates(): Error opening\n"); fprintf(stdout, "output coordinate file (%s)\n", outInfo->filename); } dumpCharmmTrajectory(outInfo->file, (charmmTrajectoryInfo *) outInfo->info, atoms, NULL, NULL, NULL, NULL, -1); } dumpCharmmTrajectory(outInfo->file, (charmmTrajectoryInfo *) outInfo->info, atoms, X, Y, Z, box, set); break; case COORD_AMBER_RESTART: if (append) sprintf(buffer, "%s.%i", outInfo->filename, set); else sprintf(buffer, "%s", outInfo->filename); if ( openFile(&outInfo->file, buffer, "w") ) { dumpAmberRestart(outInfo->file, atoms, X, Y, Z, NULL, NULL, NULL, (outInfo->isBox ? box : NULL)); } else { fprintf(stdout, "WARNING in ptrajOutputCoordinates(): Could not open\n"); fprintf(stdout, "output coordinate file %s, not dumping to output file.\n", buffer); } break; case COORD_BINPOS: /* * preprocessing; open input file and write magic header */ if (first) { if ( ! openFile(&outInfo->file, outInfo->filename, "w") ) { fprintf(stdout, "WARNING in ptrajOutputCoordinates(): Error opening\n"); fprintf(stdout, "output coordinate file (%s)\n", outInfo->filename); } fwrite( "fxyz", 4, 1, outInfo->file ); } writebinpos(outInfo->file, atoms, X, Y, Z); break; default: fprintf(stdout, "WARNING in ptrajOutputCoordinates(): Attempting to output an unknown\n"); fprintf(stdout, "trajectory type. Ignoring.\n"); ptrajCleanup(); return; } } void parseHBondDonor(stackType **argumentStackPointer, ptrajState *state, int *hbondDonor) { char *buffer, *buffer1; int *mask, i, j; Name atom, res; buffer = argumentStackKeyToString( argumentStackPointer, "mask", NULL ); /* * mask */ if (buffer != NULL) { mask = processAtomMask(buffer, state); atomMaskIsActive(mask, state, &i, &j); if (i==0) { fprintf(stdout, "WARNING in ptraj, donor: No atoms selected (%s), ignoring...\n", buffer); safe_free(buffer); safe_free(mask); return; } for (i=0; i < state->atoms; i++) { if (mask[i] != 0) { hbondDonor[i] = 1; if (prnlev > 0) fprintf(stdout, " DONOR: adding atom %5i, residue %4i, atom name %s\n", i+1, atomToResidue(i+1, state->residues, state->ipres), state->atomName[i]); } } safe_free(buffer); safe_free(mask); return; } /* * resname atomname */ buffer = getArgumentString( argumentStackPointer, NULL ); buffer1 = getArgumentString( argumentStackPointer, NULL ); if (buffer != NULL && buffer1 != NULL) { /* * copy in atom and residue names and pad with spaces as appropriate */ for (i=0; i < 4; i++) { res[i] = (char) 0; atom[i] = (char) 0; } strncpy(res, buffer, 5); strncpy(atom, buffer1, 5); res[4] = (char) 0; atom[4] = (char) 0; i = 3; while (res[i] == (char) 0) { res[i] = ' '; i--; } i = 3; while (atom[i] == (char) 0) { atom[i] = ' '; i--; } for (i=0; i < state->atoms; i++) { j = atomToResidue(i+1, state->residues, state->ipres)-1; if (strcmp(atom, state->atomName[i]) == 0 && strcmp(res, state->residueName[j]) == 0) { hbondDonor[i] = 1; if (prnlev > 0) fprintf(stdout, " DONOR: adding atom %5i, residue %4i, atom name %s\n", i+1, j+1, state->atomName[i]); } } safe_free(buffer); safe_free(buffer1); } else { fprintf(stdout, "WARNING in ptraj, donor: either the residue or atom name specification is blank!\n"); } return; } void parseHBondAcceptor(stackType **argumentStackPointer, ptrajState *state, int *hbondAcceptor, int *hbondAcceptorH1, int *hbondAcceptorH2, int *hbondAcceptorH3) { char *buffer, *buffer1; int *mask, *mask1, i, j, i1, j1, k, stop; Name atom, atom1, res; buffer1 = NULL; buffer = argumentStackKeyToString( argumentStackPointer, "mask", NULL ); if (buffer != NULL) { buffer1 = getArgumentString( argumentStackPointer, NULL ); } /* * mask */ if (buffer != NULL && buffer1 != NULL) { mask = processAtomMask(buffer, state); atomMaskIsActive(mask, state, &i, &i1); mask1 = processAtomMask(buffer1, state); atomMaskIsActive(mask1, state, &j, &j1); stop = 0; if (i==0) { fprintf(stdout, "WARNING in ptraj, acceptor: No heavy atom was selected (%s), ignoring...\n", buffer); stop = 1; } if (j==0) { fprintf(stdout, "WARNING in ptraj, acceptor: No hydrogen atom was selected (%s), ignoring...\n", buffer1); stop = 1; } if (i != j) { fprintf(stdout, "WARNING in ptraj, acceptor: There is not a 1-1 correspondence between the\n"); fprintf(stdout, "atom selection in the two masks %s and %s which contain %i and %i\n", buffer, buffer1, i, j); fprintf(stdout, "atoms respectively. Ignoring...\n"); stop = 1; } if (stop) { safe_free(buffer); safe_free(buffer1); safe_free(mask); safe_free(mask1); return; } stop = i; for (i=0; i < stop; i++) { hbondAcceptor[i1] = 1; if (hbondAcceptorH1[i1] >= 0) { if (hbondAcceptorH2[i1] >= 0) { if (hbondAcceptorH3[i1] >= 0) { fprintf(stdout, "WARNING in ptraj, acceptor: More than three hydrogens have been selected\n"); fprintf(stdout, "as acceptors for the atom %5i. Ignoring the latest...\n", i1+1); } else hbondAcceptorH3[i1] = j1; } else hbondAcceptorH2[i1] = j1; } else hbondAcceptorH1[i1] = j1; if (prnlev > 0) { fprintf(stdout, " ACCEPTOR: adding atom %5i, res %4i, atom name %s -- ", i1+1, atomToResidue(i1+1, state->residues, state->ipres), state->atomName[i1]); fprintf(stdout, "atom %5i, res %4i, atom name %s\n", j1+1, atomToResidue(j1+1, state->residues, state->ipres), state->atomName[j1]); } mask[i1] = 0; mask1[j1] = 0; while (mask[i1] == 0 && i1 < state->atoms) i1++; while (mask1[j1] == 0 && j1 < state->atoms) j1++; if (i1 == state->atoms || j1 == state->atoms) return; } safe_free(mask); safe_free(mask1); safe_free(buffer); safe_free(buffer1); return; } else { if (buffer1 != NULL || buffer != NULL) { fprintf(stdout, "WARNING in ptraj, acceptor: Error in mask specification. Ignoring...\n"); safe_free(buffer); safe_free(buffer1); return; } } /* * resname atomname atomname */ stop = 0; buffer = getArgumentString( argumentStackPointer, NULL ); if (stop || buffer == NULL) stop = 1; else { for (i=0; i < 4; i++) { res[i] = (char) 0; } strncpy(res, buffer, 5); res[4] = (char) 0; i = 3; while (res[i] == (char) 0) { res[i] = ' '; i--; } safe_free(buffer); } buffer = getArgumentString( argumentStackPointer, NULL ); if (stop || buffer == NULL) stop = 1; else { for (i=0; i < 4; i++) { atom[i] = (char) 0; } strncpy(atom, buffer, 5); atom[4] = (char) 0; i = 3; while (atom[i] == (char) 0) { atom[i] = ' '; i--; } safe_free(buffer); } buffer = getArgumentString( argumentStackPointer, NULL ); if (stop || buffer == NULL) stop = 1; else { for (i=0; i < 4; i++) { atom1[i] = (char) 0; } strncpy(atom1, buffer, 5); atom1[4] = (char) 0; i = 3; while (atom1[i] == (char) 0) { atom1[i] = ' '; i--; } safe_free(buffer); } if (stop == 0) { for (i=0; i < state->atoms; i++) { j = atomToResidue(i+1, state->residues, state->ipres)-1; if (strcmp(atom, state->atomName[i]) == 0 && strcmp(res, state->residueName[j]) == 0) { for (k=state->ipres[j]-1; k < state->ipres[j+1]-1; k++) { if (strcmp(atom1, state->atomName[k]) == 0) { hbondAcceptor[i] = 1; if (hbondAcceptorH1[i] >= 0) if (hbondAcceptorH2[i] >= 0) if (hbondAcceptorH3[i] >= 0) { fprintf(stdout, "WARNING in ptraj, acceptor: More than three hydrogens "); fprintf(stdout, "have been selected\nas acceptors for atom %i. ", i+1); fprintf(stdout, "Ignoring the latest...\n"); } else hbondAcceptorH3[i] = k; else hbondAcceptorH2[i] = k; else hbondAcceptorH1[i] = k; if (prnlev > 0) { fprintf(stdout, " ACCEPTOR: adding atom %5i, res %4i, atom name %s -- ", i+1, atomToResidue(i+1, state->residues, state->ipres), state->atomName[i]); fprintf(stdout, "atom %5i, res %4i, atom name %s\n", k+1, atomToResidue(k+1, state->residues, state->ipres), state->atomName[k]); } } } } } } else fprintf(stdout, "WARNING in ptraj, acceptor: error in specification of res/atom selection\n"); } /* * ptrajSetupIO(): This is called upon receipt of the trigger strings * that specify input and output files, reference structures and * global variables that may be used by various actions (such as solvent * information or hydrogen bonding donor/acceptor atoms). * For I/O, the transformFileStack list of input files is set up as well * as the output trajectory file. */ void ptrajSetupIO(stackType *argumentStack, actionType type) { char *filename; char *buffer; char *title, *titlenew; int i, j, k, start, stop, offset; coordinateInfo *info; ptrajState *state, **statep; int *mask, byres, bytype, curres, *startsol, *stopsol; Name res; charmmTrajectoryInfo *charmmTraj, *ctrj; stackType *sp; coordinateInfo *infiles; double boxtmp; statep = ptrajCurrentState(); state = *statep; /* * Input/Output setup for ptraj. The following command/arguments are parsed: * * TRANSFORM_BOX * * box [x value] [y value] [z value] [alpha value] [beta value] [gamma value] * [fixx] [fixy] [fixz] [fixalpha] [fixbeta] [fixgamma] * * TRANSFORM_PRNLEV * * prnlev * * TRANSFORM_REFERENCE * * reference filename * * TRANSFORM_SOLVENT * * solvent mask1 [mask2] ... [maskN] [byres | bytype | byname] * * TRANSFORM_TRAJIN * * trajin filename [start] [stop] [delta] * * TRANSFORM_TRAJOUT * * trajout filename [nobox] [ PDB | RESTART | BINPOS ] * * TRANSFORM_DONOR * * donor {print | clear | | mask } * * TRANSFORM_ACCEPTOR * * acceptor print | clear | * { } | * { mask } * */ switch ( type ) { case TRANSFORM_PRNLEV: prnlev = getArgumentInteger(&argumentStack, 0.0); fprintf(stdout, " PRNLEV: value is now %i\n", prnlev); return; case TRANSFORM_BOX: boxtmp = argumentStackKeyToDouble(&argumentStack, "x", 0.0); if (boxtmp != 0.0) { state->IFBOX = 1; state->box[0] = boxtmp; fprintf(stdout, " BOX X = %f\n", boxtmp); } boxtmp = argumentStackKeyToDouble(&argumentStack, "y", 0.0); if (boxtmp != 0.0) { state->IFBOX = 1; state->box[1] = boxtmp; fprintf(stdout, " BOX Y = %f\n", boxtmp); } boxtmp = argumentStackKeyToDouble(&argumentStack, "z", 0.0); if (boxtmp != 0.0) { state->IFBOX = 1; state->box[2] = boxtmp; fprintf(stdout, " BOX Z = %f\n", boxtmp); } boxtmp = argumentStackKeyToDouble(&argumentStack, "alpha", 0.0); if (boxtmp != 0.0) { state->IFBOX = 1; state->box[3] = boxtmp; fprintf(stdout, " BOX ALPHA = %f\n", boxtmp); } boxtmp = argumentStackKeyToDouble(&argumentStack, "beta", 0.0); if (boxtmp != 0.0) { state->IFBOX = 1; state->box[4] = boxtmp; fprintf(stdout, " BOX BETA = %f\n", boxtmp); } boxtmp = argumentStackKeyToDouble(&argumentStack, "gamma", 0.0); if (boxtmp != 0.0) { state->IFBOX = 1; state->box[5] = boxtmp; fprintf(stdout, " BOX GAMMA = %f\n", boxtmp); } if (argumentStackContains( &argumentStack, "fixx" )) state->boxfixed[0] = 1; if (argumentStackContains( &argumentStack, "fixy" )) state->boxfixed[1] = 1; if (argumentStackContains( &argumentStack, "fixz" )) state->boxfixed[2] = 1; if (argumentStackContains( &argumentStack, "fixalpha" )) state->boxfixed[3] = 1; if (argumentStackContains( &argumentStack, "fixbeta" )) state->boxfixed[4] = 1; if (argumentStackContains( &argumentStack, "fixgamma" )) state->boxfixed[5] = 1; return; case TRANSFORM_DONOR: /* * print */ if (argumentStackContains( &argumentStack, "print" )) { if (hbondDonor != NULL) { printHBondMask(hbondDonor, NULL, NULL, NULL, state); } return; } /* * clear */ if (argumentStackContains( &argumentStack, "clear" )) { if (hbondDonor != NULL) { pushBottomStack(&hbondDonorStack, (void *) hbondDonor); hbondDonor = NULL; } return; } /* * allocate space for the hbondDonor if necessary... */ if (hbondDonor == NULL) { hbondDonor = (int *) safe_malloc(sizeof(int) * state->atoms); for (k=0; k < state->atoms; k++) hbondDonor[k] = 0; } parseHBondDonor(&argumentStack, state, hbondDonor); /* * make sure at least one donor has been chosen or free the memory */ atomMaskIsActive(hbondDonor, state, &i, &j); if (i == 0) { safe_free(hbondDonor); hbondDonor = NULL; } return; case TRANSFORM_ACCEPTOR: /* * print */ if (argumentStackContains( &argumentStack, "print" )) { if (hbondAcceptor != NULL) { printHBondMask(hbondAcceptor, hbondAcceptorH1, hbondAcceptorH2, hbondAcceptorH3, state); } return; } /* * clear */ if (argumentStackContains( &argumentStack, "clear" )) { if (hbondAcceptor != NULL) { pushBottomStack(&hbondAcceptorStack, (void *) hbondAcceptor); pushBottomStack(&hbondAcceptorH1Stack, (void *) hbondAcceptorH1); pushBottomStack(&hbondAcceptorH2Stack, (void *) hbondAcceptorH2); pushBottomStack(&hbondAcceptorH3Stack, (void *) hbondAcceptorH3); hbondAcceptor = NULL; hbondAcceptorH1 = NULL; hbondAcceptorH2 = NULL; hbondAcceptorH3 = NULL; } return; } /* * allocate space for the hbondAcceptor lists if necessary... */ if (hbondAcceptor == NULL) { hbondAcceptor = (int *) safe_malloc(sizeof(int) * state->atoms); hbondAcceptorH1 = (int *) safe_malloc(sizeof(int) * state->atoms); hbondAcceptorH2 = (int *) safe_malloc(sizeof(int) * state->atoms); hbondAcceptorH3 = (int *) safe_malloc(sizeof(int) * state->atoms); for (k=0; k < state->atoms; k++) { hbondAcceptor[k] = 0; hbondAcceptorH1[k] = -1; hbondAcceptorH2[k] = -1; hbondAcceptorH3[k] = -1; } } parseHBondAcceptor(&argumentStack, state, hbondAcceptor, hbondAcceptorH1, hbondAcceptorH2, hbondAcceptorH3); /* * make sure at least one acceptor has been chosen or free the memory */ atomMaskIsActive(hbondAcceptor, state, &i, &j); if (i == 0) { safe_free(hbondAcceptor); safe_free(hbondAcceptorH1); safe_free(hbondAcceptorH2); safe_free(hbondAcceptorH3); hbondAcceptor = NULL; hbondAcceptorH1 = NULL; hbondAcceptorH2 = NULL; hbondAcceptorH3 = NULL; } return; case TRANSFORM_REFERENCE: filename = getArgumentString( &argumentStack, NULL ); if (filename == NULL) { fprintf(stdout, "WARNING in ptrajSetupIO(): reference command lacks a filename!\n"); return; } info = checkCoordinates(filename, state->atoms); safe_free(filename); if ( info == NULL ) return; info->mask = NULL; info->offset = 1; info->start = 1; info->stop = -1; info->option1 = 0; info->option2 = 0; if ( ptrajPreprocessInputCoordinates(info) ) { safe_free(info->filename); safe_free(info->info); return; } info->x = (double *) safe_malloc(sizeof(double) * state->atoms); info->y = (double *) safe_malloc(sizeof(double) * state->atoms); info->z = (double *) safe_malloc(sizeof(double) * state->atoms); /* * note, we initialize the memory here since the checks to disable processing * of truncated sets have been removed, i.e. it is possible to read in a truncated * PDB reference structure noting that the truncated part will be replaced by zeros. * If the set is truncated, the user will be warned... */ for (i=0; i < state->atoms; i++) { info->x[i] = 0.0; info->y[i] = 0.0; info->z[i] = 0.0; } ptrajProcessInputCoordinates(info, state, info->x, info->y, info->z, state->box, 1, &i, &j); if (j == 0) { safe_free(info->x); safe_free(info->y); safe_free(info->z); safe_free(info->filename); safe_free(info->info); return; } referenceInfo = info; pushBottomStack( &transformReferenceStack, (void *) info ); return; case TRANSFORM_SOLVENT: /* * solvent [byres | bytype | byname] mask1 [mask2] [mask3] ... */ byres = argumentStackContains( &argumentStack, "byres"); bytype = argumentStackContains( &argumentStack, "bytype"); if ( argumentStackContains( &argumentStack, "byname") == 1) { byres = 0; bytype = 0; } if (bytype == 1) { fprintf(stderr, "WARNING in ptrajSetupIO(): solvent \"bytype\" option has not yet\n"); fprintf(stderr, "been implemented. Defaulting to byname...\n"); bytype = 0; } buffer = (char *) argumentStack->entry; /* * if we've exhausted the argument stack, assume this is a call to * solvent with no arguments; in this case we do not alter the solvent * information, we just print a summary of it... */ if (buffer[0] == (char) 0) { ptrajPrintState(state); return; } /* * re-initialize the solvent information */ state->solventMolecules = 0; state->solventAtoms = 0; if (state->solventMask != NULL) safe_free(state->solventMask); state->solventMask = NULL; if (state->solventMoleculeStart != NULL) safe_free(state->solventMoleculeStart); state->solventMoleculeStart = NULL; if (state->solventMoleculeStop != NULL) safe_free(state->solventMoleculeStop); state->solventMoleculeStop = NULL; state->solventMask = (int *) safe_malloc(sizeof(int) * state->atoms); for (i=0; i < state->atoms; i++) state->solventMask[i] = 0; startsol = (int *) safe_malloc(sizeof(int) * state->atoms); stopsol = (int *) safe_malloc(sizeof(int) * state->atoms); buffer = NULL; while ( (buffer = getArgumentString(&argumentStack, NULL)) != NULL ) { if (byres) { /* * load up solvent by residue using the residues specified in * the input masks... */ mask = processAtomMask(buffer, state); fprintf(stdout, " Searching for solvent by mask "); printAtomMask(mask, state); for (i=0; i < state->atoms; i++) { if (mask[i]) { curres = atomToResidue(i+1, state->residues, state->ipres)-1; j = isActiveResidue(curres, mask, state); if (j > 0) { /* * all the atoms in "curres" are active, hence add this * solvent molecule to the list of solvent. Note that * isActiveResidue returns the number of active atoms * found (if not zero)... */ state->solventAtoms += j; startsol[state->solventMolecules] = i; stopsol[ state->solventMolecules] = i+j; state->solventMolecules++; for (k=i; k < i+j; k++) state->solventMask[k] = 1; i += j-1; } } } } else if (bytype) { /* * search for solvent using the mask to define * what a representative solvent molecule is... * * CURRENTLY NOT IMPLEMENTED */ mask = processAtomMask(buffer, state); } else { /* * search for solvent by residue name */ /* * copy the residue name from the buffer into "res" and pad * with spaces to conform to standard atom/residue naming */ for (i=0; i < 4; i++) res[i] = (char) 0; strncpy(res, buffer, 5); res[4] = (char) 0; i = 3; while (res[i] == (char) 0) { res[i] = ' '; i--; } fprintf(stdout, " Searching for solvent by residue name %s\n", res); /* * loop over all residues and check for matches */ for (i=0; i < state->residues; i++) { if (strcmp(res, state->residueName[i]) == 0) { /* * add this residue to the list of solvent */ j = state->ipres[i+1]-state->ipres[i]; state->solventAtoms += j; startsol[state->solventMolecules] = state->ipres[i]-1; stopsol[ state->solventMolecules] = state->ipres[i+1]-1; state->solventMolecules++; for (k=state->ipres[i]-1; k < state->ipres[i+1]-1; k++) state->solventMask[k] = 1; } } } safe_free(mask); safe_free(buffer); } /* * update the solvent information and free any unnecessary memory */ state->solventMoleculeStart = (int *) safe_malloc(sizeof(int) * state->solventMolecules); state->solventMoleculeStop = (int *) safe_malloc(sizeof(int) * state->solventMolecules); for (i=0; i < state->solventMolecules; i++) { state->solventMoleculeStart[i] = startsol[i]; state->solventMoleculeStop[i] = stopsol[i]; } safe_free(buffer); safe_free(stopsol); safe_free(startsol); if (prnlev > 0) ptrajPrintState(state); return; case TRANSFORM_TRAJIN: filename = getArgumentString( &argumentStack, NULL ); if (filename == NULL) { fprintf(stdout, "WARNING in ptrajSetupIO: trajin command lacks a filename!\n"); return; } /* * check the coordinates and push the information to the bottom * of the transformFileStack */ fprintf(stderr, "Checking coordinates: %s\n", filename); info = checkCoordinates(filename, state->atoms); if (info == (coordinateInfo *) NULL) { fprintf(stdout, "WARNING in ptrajSetupIO(): trajin %s, cannot open file...\n", filename); safe_free(filename); return; } safe_free(filename); /* * set start, stop and offset if they were specified and relevent */ start = getArgumentInteger( &argumentStack, 1 ); stop = getArgumentInteger( &argumentStack, -1 ); offset = getArgumentInteger( &argumentStack, 1 ); if (stop > 0) { if (stop > info->stop) { fprintf(stdout, "FYI ptraj(): trajin stop value (%i) is greater than the number of sets read in\n", stop); fprintf(stdout, "FYI ptraj(): Setting stop to the maximum value (%i)\n", info->stop); } else { info->stop = stop; } } if (start > 1) { if (info->stop > 0 && start > info->stop) { fprintf(stdout, "WARNING in ptraj(): trajectory start is > stop; no\n"); fprintf(stdout, "configurations will be processed\n"); } info->start = start; } info->offset = offset; if (info->offset != 1 && info->offset > info->stop - info->start) { fprintf(stdout, "WARNING in ptraj(): Offset is so large that only 1 set\n"); fprintf(stdout, " will be processed...\n"); } pushBottomStack( &transformFileStack, (void *) info ); return; case TRANSFORM_TRAJOUT: filename = getArgumentString( &argumentStack, NULL ); if (filename == NULL) { fprintf(stdout, "WARNING in ptraj(): trajout command lacks a filename!\n"); safe_free(filename); return; } /* * do not actually open the filename yet; wait until we * know whether the output is a trajectory style (i.e. multiple * sets) or a single frame */ info = (coordinateInfo *) safe_malloc(sizeof(coordinateInfo)); INITIALIZE_coordinateInfo(info); info->filename = filename; info->isBox = 1; /* * check to see if we do not want box coordinates dumped */ if (argumentStackContains( &argumentStack, "nobox" )) info->isBox = 0; /* * check to see if a format other than amber trajectory is wanted */ if (argumentStackContains( &argumentStack, "pdb" )) info->type = COORD_PDB; else if (argumentStackContains( &argumentStack, "rest" )) info->type = COORD_AMBER_RESTART; else if (argumentStackContains( &argumentStack, "binpos" )) info->type = COORD_BINPOS; else if (argumentStackContains( &argumentStack, "charmm" )) info->type = COORD_CHARMM_TRAJECTORY; else info->type = COORD_AMBER_TRAJECTORY; /* * check to see if we want charges/radii dumped to pdb */ if (info->type == COORD_PDB) { if (argumentStackContains( &argumentStack, "dumpq" )) info->option1 = 1; if (argumentStackContains( &argumentStack, "parse" )) info->option1 = 2; if (argumentStackContains( &argumentStack, "nowrap" )) info->option2 = 1; } /* * check to see if other CHARMM related information is present */ if (info->type == COORD_CHARMM_TRAJECTORY) { /* * search through the list of input files to setup the CHARMM information * structure defaults; if none is present, make it up and/or modify according * to what the user specifies... */ charmmTraj = (charmmTrajectoryInfo *) safe_malloc(sizeof(charmmTrajectoryInfo)); INITIALIZE_charmmTrajectoryInfo(charmmTraj); ctrj = NULL; for (sp = transformFileStack; ctrj == NULL && sp != NULL; sp = sp->next) { infiles = (coordinateInfo *) sp->entry; if (infiles->type == COORD_CHARMM_TRAJECTORY) ctrj = (charmmTrajectoryInfo *) infiles->info; } if (ctrj != NULL) { charmmTraj->byteorder = ctrj->byteorder; charmmTraj->magic = ctrj->magic; for (i=0; i < 20; i++) charmmTraj->icntrl[i] = ctrj->icntrl[i]; charmmTraj->ntitle = ctrj->ntitle; for (sp = ctrj->titleStack; sp != NULL; sp = sp->next) { title = (char *) sp->entry; titlenew = (char *) safe_malloc(sizeof(char) * (strlen(title) + 1)); strcpy(titlenew, title); pushBottomStack(&charmmTraj->titleStack, (void *) title); } charmmTraj->natrec = ctrj->natrec; charmmTraj->nfreat = ctrj->nfreat; if (ctrj->nfreat != ctrj->natrec) { charmmTraj->freeat = (int *) safe_malloc(sizeof(int) * ctrj->nfreat); for (i=0; i < ctrj->nfreat; i++) charmmTraj->freeat[i] = ctrj->freeat[i]; } for (i=0; i<6; i++) charmmTraj->xtlabc[i] = ctrj->xtlabc[i]; } else { charmmTraj->byteorder = 0; charmmTraj->magic.c[0] = 'C'; charmmTraj->magic.c[1] = 'O'; charmmTraj->magic.c[2] = 'R'; charmmTraj->magic.c[3] = 'D'; for (i=0; i<20; i++) charmmTraj->icntrl[i] = 0; charmmTraj->icntrl[19] = 26; if (state->IFBOX) { charmmTraj->icntrl[10] = 1; /* QCRYS */ } } if (info->isBox == 0) charmmTraj->icntrl[10] = 0; if (argumentStackContains( &argumentStack, "big" )) charmmTraj->byteorder = 0; else if (argumentStackContains( &argumentStack, "little" )) charmmTraj->byteorder = 1; /*!!!!!!!!!!!!!!!*/ info->info = (void *) charmmTraj; } outInfo = info; return; } } /* * ptrajSetup(): This routine is called for every trigger that is related * to coordinate processing (i.e. not those commands that are I/O * related, such as trajin, trajout or reference and not those that * are involved with postprocessing any acculated data). This creates * the transformActionStack stack of "actions" to be performed and * is called by dispatchToken() upon receipt of the appropriate trigger. * Most of the actual setup of the action function is performed by the * actual action function itself in the PTRAJ_SETUP mode. See the * detailed comments in actions.c for more information. */ void ptrajSetup(stackType *argumentStack, actionType type) { actionInformation *action; int ierr; /* * Allocate and initialize a actionInformation structure */ action = (actionInformation *) safe_malloc(sizeof(actionInformation)); INITIALIZE_actionInformation(action); /* * Place a copy of the current state into the action */ action->state = ptrajCopyState(ptrajCurrentState()); /* * Set the action type */ action->type = type; /* * Set the action function */ switch ( type ) { case TRANSFORM_ANGLE: action->type = TRANSFORM_ANGLE; action->fxn = (actionFunction) transformAngle; break; case TRANSFORM_ATOMICFLUCT: action->type = TRANSFORM_ATOMICFLUCT; action->fxn = (actionFunction) transformAtomicFluct; break; case TRANSFORM_AVERAGE: action->type = TRANSFORM_AVERAGE; action->fxn = (actionFunction) transformAverage; break; case TRANSFORM_CENTER: action->type = TRANSFORM_CENTER; action->fxn = (actionFunction) transformCenter; break; case TRANSFORM_CHECKOVERLAP: action->type = TRANSFORM_CHECKOVERLAP; action->fxn = (actionFunction) transformCheckOverlap; break; case TRANSFORM_CLOSESTWATERS: action->type = TRANSFORM_CLOSESTWATERS; action->fxn = (actionFunction) transformClosestWaters; break; case TRANSFORM_CORRELATION: action->type = TRANSFORM_CORRELATION; action->fxn = (actionFunction) transformCorr; break; case TRANSFORM_DIHEDRAL: action->type = TRANSFORM_DIHEDRAL; action->fxn = (actionFunction) transformDihedral; break; case TRANSFORM_DIFFUSION: action->type = TRANSFORM_DIFFUSION; action->fxn = (actionFunction) transformDiffusion; break; case TRANSFORM_DIPOLE: action->type = TRANSFORM_DIPOLE; action->fxn = (actionFunction) transformDipole; break; case TRANSFORM_DISTANCE: action->type = TRANSFORM_DISTANCE; action->fxn = (actionFunction) transformDistance; break; case TRANSFORM_GRID: action->type = TRANSFORM_GRID; action->fxn = (actionFunction) transformGrid; break; case TRANSFORM_HBOND: action->type = TRANSFORM_HBOND; action->fxn = (actionFunction) transformHBond; break; case TRANSFORM_IMAGE: action->type = TRANSFORM_IMAGE; action->fxn = (actionFunction) transformImage; break; case TRANSFORM_PRINCIPAL: action->type = TRANSFORM_PRINCIPAL; action->fxn = (actionFunction) transformPrincipal; break; case TRANSFORM_PUCKER: action->type = TRANSFORM_PUCKER; action->fxn = (actionFunction) transformPucker; break; case TRANSFORM_RADIAL: action->type = TRANSFORM_RADIAL; action->fxn = (actionFunction) transformRadial; break; case TRANSFORM_RMS: action->type = TRANSFORM_RMS; action->fxn = (actionFunction) transformRMS; break; case TRANSFORM_RUNNINGAVERAGE: action->type = TRANSFORM_RUNNINGAVERAGE; action->fxn = (actionFunction) transformRunningAverage; break; case TRANSFORM_SCALE: action->type = TRANSFORM_SCALE; action->fxn = (actionFunction) transformScale; break; case TRANSFORM_STRIP: action->type = TRANSFORM_STRIP; action->fxn = (actionFunction) transformStrip; break; case TRANSFORM_TRANSLATE: action->type = TRANSFORM_TRANSLATE; action->fxn = (actionFunction) transformTranslate; break; case TRANSFORM_TRUNCOCT: action->type = TRANSFORM_TRUNCOCT; action->fxn = (actionFunction) transformTruncOct; break; case TRANSFORM_VECTOR: action->type = TRANSFORM_VECTOR; action->fxn = (actionFunction) transformVector; break; case TRANSFORM_WATERSHELL: action->type = TRANSFORM_WATERSHELL; action->fxn = (actionFunction) transformWatershell; break; case TRANSFORM_2DRMS: action->type = TRANSFORM_2DRMS; action->fxn = (actionFunction) transform2dRMS; break; case TRANSFORM_TRANSFORM: case TRANSFORM_NOOP: action->type = type; action->fxn = NULL; break; default: fprintf(stdout, "WARNING in ptrajSetup(): unknown action type %i\n", type); error("ptraj()", "terminating...\n"); } /* * Parse the arguments. This is done by the action->fxn in the * PTRAJ_SETUP mode. The argumentStack is placed into the * complex argument 1 slot for this. */ ierr = 0; if (action->type != TRANSFORM_TRANSFORM && action->type != TRANSFORM_NOOP) { action->carg1 = (void *) &argumentStack; ierr = action->fxn(action, NULL, NULL, NULL, NULL, PTRAJ_SETUP); } /* * If the setup fails, -1 is returned and therefore this action * should not be placed on the action stack and the associated * memory should be freed */ if (ierr < 0) { safe_free(action->state); action->state= NULL; safe_free(action); action = NULL; } else { /* * Place the now setup action structure onto the transformActionStack */ pushBottomStack( &transformActionStack, (void *) action ); } } void ptrajSetupAnalyze(stackType *argumentStack, actionType type) { analyzeInformation *analyze; char *buffer; int ierr; /* * Make sure that this is indeed an "analyze" action (which really isn't * necessary) */ if (type != TRANSFORM_ANALYZE) { fprintf(stdout, "Error in ptrajSetupAnalyze(): Called with the wrong type!\n"); fprintf(stdout, "Ignoring this command...\n"); return; } /* * Grab the first argument off the argument stack. This is the "trigger" for * the analyze function */ buffer = getArgumentStringLower(&argumentStack, NULL); if (buffer == NULL) { fprintf(stdout, "ptrajSetupAnalyze(): No command passed to analyze, ignoring...\n"); return; } /* * Allocate and initialize a analyzeInformation structure and set the type */ analyze = (analyzeInformation *) safe_malloc(sizeof(analyzeInformation)); INITIALIZE_analyzeInformation(analyze); /* * search for a match to the trigger (stored in buffer) */ if (strncmp(buffer, "correlationcoe", 14) == 0) { analyze->type = ANALYZE_CORRELATIONCOEFFICIENT; analyze->fxn = (analyzeFunction) analyzeCorrelationCoefficient; } else if (strncmp(buffer, "hbond", 5) == 0) { analyze->type = ANALYZE_HBOND; analyze->fxn = (analyzeFunction) analyzeHBond; } else if (strncmp(buffer, "set", 3) == 0) { analyze->type = ANALYZE_SET; analyze->fxn = (analyzeFunction) analyzeSet; } else if (strncmp(buffer, "stat", 4) == 0) { analyze->type = ANALYZE_STATISTICS; analyze->fxn = (analyzeFunction) analyzeStatistics; } else if (strncmp(buffer, "test", 4) == 0) { analyze->type = ANALYZE_TEST; analyze->fxn = (analyzeFunction) analyzeTest; } else { fprintf(stdout, "WARNING in ptrajSetupAnalyze(): unknown analyze type %i\n", type); safe_free(analyze); return; } /* * Parse the arguments. This is done by the analyze->fxn itself in the * PTRAJ_SETUP mode. The argumentStack is placed into the * complex argument 1 slot for this. */ ierr = 0; if (analyze->type != ANALYZE_NOOP) { analyze->carg1 = (void *) &argumentStack; ierr = analyze->fxn(analyze, scalarStack, PTRAJ_SETUP); } /* * If the setup fails, -1 is returned and therefore this action * should not be placed on the action stack and the associated * memory should be freed */ if (ierr < 0) { safe_free(analyze); analyze = NULL; } else { /* * Place the now setup analyze structure onto the transformAnalyzeStack */ pushBottomStack( &transformAnalyzeStack, (void *) analyze ); } } void ptraj(char *filenamep) { FILE *infile; char buffer[BUFFER_SIZE]; char *bufferp; coordinateInfo *currentCoordinateInfo; double *X, *Y, *Z; double box[6], boxnew[6]; int boxfixed[6]; actionInformation *action; analyzeInformation *analyze; int set; int local_set; int processed; int i; int suppressProcessing; stackType *actionStackTemp = NULL; int outputTrajectory = 0; stackType *sp, *argumentStack; ptrajState *startingState, *currentState, **statep; char *continuation; int readCoordinates; int processCoordinates; int firstOutput; /* * --------------- INPUT FILE PROCESSING -------------------- * * if an input file was specified, open it up, else use * standard input... */ fprintf(stdout, "\nPTRAJ: Processing input file...\n"); if ( filenamep == NULL || strcmp(filenamep, "") == 0 || strcmp(filenamep, "stdin") == 0 ) { fprintf(stdout, " Input is from standard input\n"); infile = stdin; } else if ( openFile(&infile, filenamep, "r") == 0 ) { fprintf(stdout, "WARNING in ptraj(): Could not open input file (%s), exiting\n", filenamep); ptrajCleanup(); return; } else { fprintf(stdout, " Input is from file %s\n", filenamep); } /* * ----------------- SETUP INITIAL STATE -------------------- * * this gives a snapshot of the current "state" based on the * appropriate parameter/topology information (i.e. AMBER prmtop) * upon entry. Note this assumes that the GLOBAL ptrajState * information was previously set in the main routine (main.c) * via a call to ptrajInitializeState(). */ statep = ptrajCurrentState(); startingState = *statep; currentState = startingState; /* * --------------- PROCESS THE INPUT FILE ------------------- * * Read the input file, line by line, using the "dispatchToken()" * routine (dispatch.c) to find a match for each command, ignoring * comments in the input file. Each "command" typed by the user * has an associated routine that is called for setup. Currently, * this is ptrajSetup() for "actions" and ptrajSetupIO for input/ * output functions. Note that ptrajSetup() will actually call the * "action" routine to perform the setup and parse arguments. Note * also that the global state pointer may be altered!!! * * Lines of text are processed until EOF is encountered. */ argumentStack = NULL; while ( (bufferp = fgets(buffer, BUFFER_SIZE, infile)) != NULL && strcmp(bufferp, "go\n") != 0 ) { continuation = bufferp; while ( continuation != NULL ) { if (strlen(buffer) >= BUFFER_SIZE) continuation = NULL; else { continuation = strrchr(buffer, '\\'); if (continuation) continuation = fgets(continuation, (BUFFER_SIZE - strlen(buffer) - 1), infile); } } skipWhitespace(bufferp); /* * skip blank lines and/or comments denoted by "#" or "%" */ if (bufferp[0] != (char) 0 && bufferp[0] != '#' && bufferp[0] != '%') { fprintf(stdout, "\nPTRAJ: %s", buffer); dispatchToken((Token *) &ptrajTokenlist, argumentStack, (char *) buffer); } } if (infile != stdin) safe_fclose(infile); /* * ----------------- ERROR CHECKING --------------------- * */ if (transformFileStack == NULL) { fprintf(stdout, "WARNING in ptraj(): No input trajectories specified (trajin), aborting...\n"); ptrajCleanup(); return; } if (outInfo == NULL) { fprintf(stdout, "FYI: No output trajectory specified (trajout), none will be saved.\n"); } else { outputTrajectory = 1; } /* * set default box information. This will allow, at setup time, specification * of the current box sizes (if none were specified such as can happen when * reading CHARMM PSF files) and also the option to FIX box sizes. */ if (currentState->IFBOX) { for (i=0; i<6; i++) { box[i] = currentState->box[i]; boxfixed[i] = currentState->boxfixed[i]; } } else { box[0] = 0.0; box[1] = 0.0; box[2] = 0.0; box[3] = 90.0; box[4] = 90.0; box[5] = 90.0; boxfixed[0] = 0; boxfixed[1] = 0; boxfixed[2] = 0; boxfixed[3] = 0; boxfixed[4] = 0; boxfixed[5] = 0; } /* * --------- CHECK HOW MANY FRAMES WILL BE PROCESSED ------- * * How many frames will be used? */ startingState->maxFrames = 0; for (sp = transformFileStack; sp != NULL; sp = sp->next) { currentCoordinateInfo = (coordinateInfo *) sp->entry; startingState->maxFrames += (currentCoordinateInfo->stop - currentCoordinateInfo->start) / currentCoordinateInfo->offset + 1; } fprintf(stdout, "\nPTRAJ: Successfully read the input file.\n"); fprintf(stdout, " Coordinate processing will occur on %i frames.\n", startingState->maxFrames); fprintf(stdout, " Summary of I/O and actions follows:\n\n"); /* * ----- PRINT A SUMMARY OF THE FILES IN/OUT AND ACTIONS ------ */ fprintf(stdout, "INPUT COORDINATE FILES\n"); printStack( &transformFileStack, printCoordinateInfo, NULL ); fprintf(stdout, "\nOUTPUT COORDINATE FILE\n"); printCoordinateInfo( (void *) outInfo ); if (transformReferenceStack != NULL) { fprintf(stdout, "\nREFERENCE FILE\n"); printStack( &transformReferenceStack, printCoordinateInfo, NULL ); } fprintf(stdout, "\nACTIONS\n"); if ( transformActionStack == NULL ) { fprintf(stdout, " Stack is NULL\n"); } else { i = 1; for (actionStackTemp = transformActionStack; actionStackTemp != NULL; actionStackTemp = actionStackTemp->next) { action = (actionInformation *) actionStackTemp->entry; /* * With each action's state variable local, it is necessary * to update each with the current states maxFrames value!!! */ action->state->maxFrames = startingState->maxFrames; /* * call the action function with the mode PTRAJ_STATUS */ if (action->type != TRANSFORM_NOOP && action->type != TRANSFORM_TRANSFORM) { fprintf(stdout, "%3i>", i); for (i=0; i < 6; i++) boxnew[i] = box[i]; /* protect the current box info */ action->fxn(action, X, Y, Z, boxnew, PTRAJ_STATUS); i++; } } fprintf(stdout, "\n"); } if (transformAnalyzeStack != NULL) { fprintf(stdout, "\nANALYZE\n"); i = 1; for (actionStackTemp = transformAnalyzeStack; actionStackTemp != NULL; actionStackTemp = actionStackTemp->next) { analyze = (analyzeInformation *) actionStackTemp->entry; if (analyze->type != ANALYZE_NOOP) { fprintf(stdout, "%3i>", i); analyze->fxn(analyze, scalarStack, PTRAJ_STATUS); i++; } } fprintf(stdout, "\n"); } /* * ---------- ALLOCATE SPACE FOR COORDINATES ------------ * * perform and initial setup necessary prior to reading in * the coordinates */ X = (double *) safe_malloc(sizeof(double) * startingState->atoms); Y = (double *) safe_malloc(sizeof(double) * startingState->atoms); Z = (double *) safe_malloc(sizeof(double) * startingState->atoms); /* * --------- MAIN LOOP FOR COORDINATE PROCESSING ------------- * * loop over each of the files representing the coordinates. * * "set" -- the global counter over all sets, all files * "local_set" -- the counter over sets in each individual file */ processed = 0; firstOutput = 1; set = 1; for (sp = transformFileStack; sp != NULL; sp = sp->next) { currentCoordinateInfo = (coordinateInfo *) sp->entry; /* * ------------- PREPROCESS COORDINATE FILES----------------- */ /* * open up the file and preprocess the coordinates */ if ( ptrajPreprocessInputCoordinates(currentCoordinateInfo) ) { ptrajCleanup(); return; } local_set = 1; /* * -------- READ IN, PROCESS, AND OUTPUT COORDINATES -------- * * The following two variables control whether to continue * reading configurations from this file or not and also * whether the coordinates should be processed after this * read... * * "readCoordinates" > 0 if there are more coordinate sets * to read in this file * * "processCoordinates" > 0 if this coordinate set should be * processed */ readCoordinates = 1; processCoordinates = 1; while ( readCoordinates ) { /* * read in the current file of coordinates, a single set at a time. */ for (i=0; i<6; i++) boxnew[i] = box[i]; ptrajProcessInputCoordinates(currentCoordinateInfo, startingState, X, Y, Z, boxnew, local_set, &readCoordinates, &processCoordinates); for (i=0; i<6; i++) if (boxfixed[i] == 0) box[i] = boxnew[i]; /* * process coordinates if necessary */ if ( processCoordinates ) { /* * check to see if this snapshot is within bounds/offset */ if ((local_set >= currentCoordinateInfo->start && local_set <= currentCoordinateInfo->stop) && ((currentCoordinateInfo->offset == 1) || ((local_set - currentCoordinateInfo->start) % currentCoordinateInfo->offset == 0)) ) { /* * TRAVERSE THE ACTION STACK to perform each * action on each coordinate set. Note that a * particular action can suppress processing of * further actions on the stack and prevent output * by setting the suppressProcessing flag in the * actionInformation structure. This is useful * when various sets are to be accumulated (for * example when calculating running averages) prior * to output and further processing... */ suppressProcessing = 0; for (actionStackTemp = transformActionStack; actionStackTemp != NULL; actionStackTemp = actionStackTemp->next) { action = (actionInformation *) actionStackTemp->entry; if (action->type != TRANSFORM_NOOP && action->type != TRANSFORM_TRANSFORM && suppressProcessing == 0) { for (i=0; i<6; i++) boxnew[i] = box[i]; /* protect box coordinates */ action->fxn(action, X, Y, Z, boxnew, PTRAJ_ACTION); for (i=0; i<6; i++) if (boxfixed[i] == 0) box[i] = boxnew[i]; /* * update the current state */ currentState = action->state; /* * check if any of the actions have suppressed output */ if (suppressProcessing == 0) suppressProcessing = action->suppressProcessing; } } /* * perform output as necessary */ if ( outputTrajectory && suppressProcessing == 0) { ptrajOutputCoordinates(currentState, set, 1, firstOutput, currentState->atoms, X, Y, Z, box); firstOutput = 0; } /* * print out a dot for each snapshot processed to show progress */ processed++; if ( local_set % 50 == 0 || local_set == 1 ) fprintf(stdout, "\nSet %6i ", local_set); if ( local_set % 50 != 0 ) fprintf(stdout, "."); } else { /* * for coordinates that were not processed, print a space */ if ( local_set % 50 == 0 || local_set == 1 ) fprintf(stdout, "\nSet %6i ", local_set); if ( local_set % 50 != 0 ) fprintf(stdout, " "); } fflush(stdout); set += 1; local_set += 1; } /* IF (processCoordinates) */ } /* WHILE (readCoordinates) */ } /* FOR over transformFileStack */ /* * -------------- FINAL POSTPROCESSING ----------------- */ if (outputTrajectory) { switch (outInfo->type) { case COORD_CHARMM_TRAJECTORY: case COORD_AMBER_TRAJECTORY: case COORD_BINPOS: safe_fclose(outInfo->file); outInfo->file = NULL; break; } } /* * -------------- DUMP ACCUMULATED DATA ------------------- */ fprintf(stdout, "\n\nPTRAJ: Successfully read in %i sets and processed %i sets.\n", set-1, processed); fprintf(stdout, " Dumping accumulated results (if any)\n\n"); for (actionStackTemp = transformActionStack; actionStackTemp != NULL; actionStackTemp = actionStackTemp->next) { action = (actionInformation *) actionStackTemp->entry; if (action->type != TRANSFORM_NOOP && action->type != TRANSFORM_TRANSFORM) { for (i=0; i<6; i++) boxnew[i] = box[i]; action->fxn(action, X, Y, Z, boxnew, PTRAJ_PRINT); } } /* * ----------- PERFORM ANY REQUESTED ANALYSIS ------------- */ if (transformAnalyzeStack != NULL) { fprintf(stdout, "\nPTRAJ: Analyzing accumulated data\n"); for (actionStackTemp = transformAnalyzeStack; actionStackTemp != NULL; actionStackTemp = actionStackTemp->next) { analyze = (analyzeInformation *) actionStackTemp->entry; if (analyze->type) analyze->fxn(analyze, scalarStack, PTRAJ_ACTION); } } ptrajCleanup(); ptrajClearState(&startingState); safe_free(X); safe_free(Y); safe_free(Z); X = NULL; Y = NULL; Z = NULL; }