/*{{{--probe.c general comments */ /* name: probe.c */ /* author: J. Michael Word */ /* date written: 2/26/96 */ /* purpose: compute intersecton surfaces */ /*****************************************************************/ /* NOTICE: This is free software and the source code is freely */ /* available. You are free to redistribute or modify under the */ /* conditions that (1) this notice is not removed or modified */ /* in any way and (2) any modified versions of the program are */ /* also available for free. */ /* ** Absolutely no Warranty ** */ /* Copyright (C) 1996-2013 J. Michael Word */ /*****************************************************************/ /*Updated to work with the remediated PDB names rmi 070727 */ /* Essentially ennumerated additional possible atom names and added */ /* the residue name info to the routine that assigns element type */ /*Jan 2006 DCR compromise comments to enable jEdit type folding {{{ }}} */ /* without messing up vi type beginning-end bracket identification shift % */ /* So jEdit folds have to stay within code regions, esp subroutines 060129*/ /*Oct 2004 DCR modifications to version 2.10.031014 starting with dcr041007 */ /*dcr: sorry Mike, but changing to Berkeley-Altman brackets really helps me!*/ /*i.e. some routines changed: opening and closing brackets vertically in line*/ /*also substituted spaces for tabs which are not robust across text editors*/ /*also tend to shorten lines to <=80 col, standard card-image window size*/ /*which many text editors use as a default width*/ /*050111: in trying to track behavior of autobondrot mode, the lack of */ /*global mode state is a severe limitation. */ /*probe trys to be so object oriented as to not have global modal states */ /*but given the existance of a whole slew of probe.c global parameters, */ /*a few more would only help clear up dependencies and program flow! */ /*hence logical Lautobondrot as a probe.c global*/ /* modifications: see dump_changes() */ /*}}}--probe.c general comments */ /*{{{--includes and declarations (globals for probe.c) */ #include #include #include #include #include "probe.h" #define INLINE_FOR_SPEED 1 static char* versionString = "probe: version 2.26.021123, Copyright 1996-2016, J. Michael Word; 2021-2023 Richardson Lab"; static char* shortVersionStr = "probe.2.26.021123"; static char *referenceString = "Word, et. al. (1999) J. Mol. Biol. 285, 1711-1733."; static char *electronicReference = "http://kinemage.biochem.duke.edu/"; static int LweakHbonds = FALSE; /*111215dcr global Logical to group weak H-bonds*/ static int LworseOverlap = FALSE; /* 04/09/2015 SJ - to separate the bad and worse overlaps if set to true*/ static int LOneDotEach = FALSE; /*111013 global One Dot for Each neighbor */ static int Lnearest = FALSE; /*111021dcr selection only on nearest neighbor*/ static int Lfilternearest = FALSE;/*111022dcr final filter nearest by extraPat*/ static int Lfilterocc1 = FALSE; /*20120120dcr final filter use only OCC==1 */ static int Lfilternoalt = FALSE; /*20120120dcr final filter exclude any alts */ static int Lgapbins = FALSE; /*111020dcr gapbin ptmasters*/ static int LXHvector = FALSE; /*20120120dcr for H: vector parent--H--Target angle*/ static int Ldotdump = FALSE; /* global flag for Ldotdump messages in examineDots*/ static int Lhonlyocclude = FALSE; /*130427dcr global flag: only Hatoms can occlude*/ static pattern *extraPat = NULL;/*111022dcr finally filter nearest by extraPat*/ static pattern *causePat = NULL;/*111022dcr initial select nearest by targPat*/ static pattern *sourcePat= NULL;/*111021dcr select src bases on srcPat */ static int LMasterName = FALSE; /*global extra master={name} on lists 060129*/ static int LMergeContacts = TRUE; /*global combine wide & close contacts 060129*/ /*static int Lautobondrot = FALSE;*/ /* global flag for AUTOBONDROT mode 050111*/ static int ImplicitH = FALSE; /* global controling dot radii */ static int UsePolarH = TRUE; /* global controling VDW radii of polar Hs */ static int Verbose = TRUE; /* global flag for output messages */ static int DoMcMc = FALSE; /* global controling mainchain->mainchain dots */ static int DoHet = TRUE; /* global controling non-Water HET dots */ static int DoH2O = TRUE; /* global controling Water dots */ static int DoWatWat = FALSE; /* global controling water->water dots */ static int LimitDots = TRUE; /* global controling chopping around bumps*/ static int LensDots = FALSE; /* global controling lens keyword in kin file*/ static int ColorGap = TRUE; /* global controling coloration of Gaps */ static int ByNABaseColor= FALSE; /* global controling coloration of dots */ static int DumpNewHO = FALSE; /*global controling dump of water newH atom data*/ static int Maxbonded = 4; /* global controling exclusion bond count */ static int UseStdBond= FALSE; /* global flag: if TRUE std AA and NA bonding assumed */ static int UseHParent= TRUE; /* global flag: bonding based on parent names on std H atoms */ static int ShowTicks = TRUE; /* global controling display of residue name ticker */ static int OldStyleU = FALSE; /* global controling -u output (true gives kiss edge stuff) */ static int SplitBySegID=FALSE;/* global controling the splitting of residues (-SEGID) */ static int HB2aromFace=TRUE; /* global controling aromatic face Hbonding */ static int AtomMasters=FALSE; /* global controling use of masters */ static int PermitCHXHB=FALSE; /* global controling use CHX (CHO) hbonds */ static int DebugLevel= 0; /* global controling debugging information */ int NuclearRadii = FALSE; /* global controling Hydrogen vdW radii, cannot be static */ #ifdef JACKnANDREA static int writeHbonds = FALSE; /* global flag: if true, write out Hbonds every time they are detected */ #endif static int OutputFmtType = 0; /* global selecting output format */ /* 0 ==> kinemage format */ /* 1 ==> O format */ /* 2 ==> XtalView format */ /* 3 ==> oneline multi:count:dcr041101*/ static int OutputHBs = TRUE; /* global controling display of contact categories */ static int OutputClashes = TRUE; /* global controling display of contact categories */ static int OutputVDWs = TRUE; /* global controling display of contact categories */ static int OnlyBadOut = FALSE; /* global controling display of contact categories dcr041010*/ static int ContactSUMMARY = FALSE; /* global controling output of contact summaries dcr041101*/ static float Min_regular_hb_cutoff=0.6f; /* globals controling hbond cutoff */ static float Min_charged_hb_cutoff=0.8f; /* defaults set in processCommandline() */ static float RadScaleFactor =1.0; /* global VDW radius scale Factor r*f */ static float RadScaleOffset =0.0; /* global VDW radius scale Offset r+o */ static float CORadScale =(float)(1.65/ATOMC_EXPLICIT_VDW); /* global VDW radius scale Factor for C=O */ static float GAPweight =0.25;/* global raw GAP score weight */ static float BUMPweight =10.0;/* global raw BUMP score scale Factor */ static float HBweight = 4.0;/* global raw HBond score scale Factor */ static float CHOHBfactor = 0.5;/* global CH..O HBond score scale Factor */ static float LowGoodCut =-0.4f;/* global cutoff for good bin */ static float HighGoodCut = 0.25;/* global cutoff for good bin */ static float WorseBumpCut =-0.5; /* SJ 04/10/2015 - global cutoff for worse overlap*/ static float OccupancyCutoff =0.02f; /* global occupancy below which atom has presence but no clashes */ /*unfortunately, one needs intervening atoms to avoid clashes between atoms*/ /*that are separated by <= Maxbonded, which thus needs to pass through any*/ /*atom no matter what it's own occupancy happens to be */ /*050118 atom->occ < OccupancyCutoff neither clash nor transfer bonded info*/ /* so flanking atoms show horrible spurious clashes (e.g. 1BKR lys 5 lys 78*/ /*mechanism seems to be atom->flags | IGNORE_FLAG */ /*050119 atoms irrespective of occ now put into neighbor bins for bonding*/ /* but atom->occ < OccupancyCutoff are not allowed to show clashes*/ static int OccupancyCriteria = TRUE; /*global: use occupancy criteria 060831*/ /*explicit default: unset by -noocc flag, 060831*/ static char *OutPointColor = "gray"; static pointSet COdots; /* adjustible radius dots for C=O carbon */ /*Old: dcr041020 NODEWIDTH defined 5 probe.h, hard code 6 to match initialization*/ /*Old: 0 wide contact,1 close contact,2 small overlap,3 bad overlap,4 H-bonds*/ /*Old: 5 running sum which can be tested for existance of any of these spikes*/ /* New 04/15/2015 SJ: 0 wide contact,1 close contact,2 weak H bonds, 3 small overlap, 4 bad overlap, 5 worse overlap, 6 H-bonds*/ /* changing the size of the array to 8*/ static long mcmccont[2][8] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; /*dcr041017,dcr041020, SJ 04/15/2015*/ static long scsccont[2][8] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; /*dcr041017,dcr041020, SJ 04/15/2015*/ static long mcsccont[2][8] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; /*dcr041017,dcr041020, SJ 04/15/2015*/ static long othrcont[2][8] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; /*dcr041017,dcr041020, SJ 04/15/2015*/ static long summcont[2][8] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; /*dcr041017,dcr041020, SJ 04/15/2015*/ /* comment: the above number 8 should not be hardcoded, maybe initialize in the main function? - SJ 04/09/2015*/ static char inputfilename[256]; /*dcr041023 global, subroutines can report it*/ static int modelNumber[100]; /*global record of incoming model # 041114*/ static int modelCount = 0; /*global count of such records, all files 041114*/ static int modelLimit = 0; /*global loop's flag for multiple models 041114*/ /*model counting simple-minded only good for one input file cases like SELF*/ static int modelSrc = 0; /*global model specified for source 041114*/ static int modelTarg = 0; /*global model specified for target 041114*/ static int modelToProcess = 0; /*global model specified for processing 041114*/ static char* dumpFileName = 0; /* global name of file to dump atom info to (default none) */ #define ATOMS_IN_THE_SAME_RES(aa, bb) IS_THE_SAME_RES((aa)->r, (bb)->r) /*041112 IS_THE_SAME_RES() cannot avoid modeltest, i.e. implied Lmodeltest==1*/ #define IS_THE_SAME_RES(ra, rb) \ ( ((ra) == (rb)) \ ||( (ra)->resid == (rb)->resid \ && (ra)->resInsCode == (rb)->resInsCode \ && strcmp((ra)->chain, (rb)->chain) == 0 \ && (ra)->model == (rb)->model ) \ && ((SplitBySegID == FALSE) || ( strcmp((ra)->segid, (rb)->segid) == 0) ) ) #define TRACEATOM(where, txt, atptr, endtxt) \ (fprintf((where), "%s{%s%c %s%d%c}[%d]%s", (txt),\ (atptr)->atomname, (atptr)->altConf,\ (atptr)->r->resname, (atptr)->r->resid,\ (atptr)->r->resInsCode, (atptr)->r->rescnt, (endtxt))) /*}}}--includes and declarations */ /*{{{main()***************************** MAIN ********************************/ int main(int argc, char **argv) { int rc = 0; FILE *outf = stdout; initalizeAtomTbl(); /*atomprops.c*/ initStdConnTable(); /*stdconntable.c*/ inputfilename[0] = '\0'; /*global dcr041023*/ rc = mainProbeProc(argc, argv, outf); return rc; } /*}}}main()__________________________________________________________________*/ /*{{{mainProbeProc()******** called from main() ******************************/ int mainProbeProc(int argc, char **argv, FILE *outf) {/*mainProbeProc()*/ int method; int keepUnselected; region boundingBoxA; atom *a, *allMainAtoms = NULL; atomBins *abins = NULL; pointSet dots[NUMATOMTYPES]; char *srcArg=NULL, *targArg=NULL, *extraArg=NULL, *ignoreArg=NULL; /* extraArg dcr 111022*/ char *groupLabel=NULL; pattern *srcPat = NULL, *targPat = NULL, *ignorePat = NULL; float density, probeRad, spikelen; int drawSpike, countDots, rawOutput, sayGroup, addKinToFile, conFlag; //conFlag by SJ - 10/07/2011 char message[200]; movingAtomBuildInfo mabis; residue *resLst = NULL; FILE* dumpFile = NULL; /*mabis moving atom build info structure */ /*buffer to pass data to newMovingAtom for autobondrot*/ allMainAtoms = processCommandline(argc, argv, &method, &boundingBoxA, &density, &probeRad, &drawSpike, &spikelen, &countDots, &keepUnselected, &srcArg, &targArg, &extraArg, &ignoreArg, &groupLabel, &rawOutput, &conFlag, &sayGroup, &addKinToFile, &mabis, &resLst);//conFlag by SJ -10/07/2011 /*called with address of mabis, i.e. mabip moving atom build info ptr */ /*autobondrot mode: not read in any mobile atoms yet*/ /*though may have static set of atoms already in. */ /* dumpRes(resLst); */ if ((allMainAtoms == NULL) && (mabis.inf == NULL)) { /*static atoms read into allMainAtoms linked list, */ /* and/or autobondrot file name for mobile atoms held in mabis.inf */ sprintf(message, "No atom data in input."); note(message); } else {/*there are atoms for probe to work on...*/ if (Verbose) {/*Verbose: lots of stuff to screen...*/ note(versionString); sprintf(message, "density: %g, probe radius: %g, maxBonded: %d", density, probeRad, Maxbonded); note(message); sprintf(message, "regular HB cutoff: %g, charged HB cutoff: %g", Min_regular_hb_cutoff, Min_charged_hb_cutoff); note(message); sprintf(message, "Dot gap bins: low to %g, %g to 0, 0 to %g, %g to high, & Hbonds", LowGoodCut, LowGoodCut, HighGoodCut, HighGoodCut); note(message); if (RadScaleOffset < -0.0001 || RadScaleOffset > 0.0001) { sprintf(message, "Add %g to VDW Rradius", RadScaleOffset); note(message); } if (RadScaleFactor < 0.9999 || RadScaleFactor > 1.0001) { sprintf(message, "Scale VDW Rradius by %g", RadScaleFactor); note(message); } sprintf(message, "C=O carbon VDW scaled by %.3f to a radius of %g A", CORadScale, CORadScale* ATOMC_EXPLICIT_VDW); note(message); if (rawOutput || countDots) { sprintf(message, "Score Weights: gapWt=%g, bumpWt=%g, HBWt=%g", GAPweight, BUMPweight, HBweight); note(message); } if (drawSpike) { sprintf(message, "draw spike, len: %g", spikelen); note(message); } if (ImplicitH) { note("implicit hydrogens"); } if (PermitCHXHB) { sprintf(message, "CH..O type hbonds recognized (scale factor %g)", CHOHBfactor); note(message); } if (!DoMcMc) { note("Excluding mainchain-mainchain interactions"); } if (!DoHet && !DoH2O) { note("Excluding both HET groups and waters"); } else { if (!DoHet) { note("Excluding non-water HET groups"); } if (!DoH2O) { note("Excluding waters"); } else if (!DoWatWat) { note("Excluding water-water interactions"); } } if (LimitDots) { note("Limiting dots from bumps to max non-bump cap"); } if (rawOutput) { /* nothing required for raw format */ } else if (OutputFmtType == 1) { note("Formatting for display in O"); } else if (OutputFmtType == 2) { note("Formatting for display in XtalView"); } else if (OutputFmtType == 3) { /*dcr041101*/ note("outputing one line colon:deliminated:counts:by:severity:type:"); } if (ByNABaseColor) { if (ColorGap) { note("Coloring dots by gap (& list color by nucleic acid base)"); } else { note("Coloring dots by nucleic acid base"); } } else if (ColorGap) { note("Coloring dots by gap (& list color by atom)"); } if (!UsePolarH) { note("Using old (long) radii for polar H atoms"); } if (!OutputHBs) { note("Attention: ***Output not generated for H-Bond contacts***"); } if (!OutputClashes) { note("Attention: ***Output not generated for clashes***"); } if (!OutputVDWs) { note("Attention: ***Output not generated for van der Waals contacts***"); } if (OnlyBadOut) { note("Attention: ***Output only generated for Bad Clashes***");/*dcr041010*/ } if (!keepUnselected) { note("Attention: ***Unselected atoms dropped (ignored)***"); } }/*Verbose: lots of stuff to screen...*/ /*now do some work...*/ if( addKinToFile && (OutputFmtType == 0) && (! countDots) && (! rawOutput) ) { fprintf(outf, "@kinemage 1\n"); } srcPat = getPat(srcArg, "pattern 1", Verbose); sourcePat = getPat(srcArg, "pattern 1", Verbose);/*111021dcr*/ modelSrc = modelInPattern(srcPat); /*parse.c 041114*/ if(Verbose && modelSrc > 0) { sprintf(message, "specified src model==%d", modelSrc); note(message); } targPat = getPat(targArg, "pattern 2", Verbose); causePat = getPat(targArg, "pattern 2", Verbose); /*111021dcr nearest*/ modelTarg = modelInPattern(targPat); /*parse.c 041114*/ if(Verbose && modelSrc > 0) { sprintf(message, "specified targ model==%d", modelTarg); note(message); } extraPat = getPat(extraArg,"extra", Verbose); /*111023dcr*/ ignorePat = getPat(ignoreArg,"ignore", Verbose); loadDotSpheres(dots, density); if (allMainAtoms) /*atoms usually read in from processCommandline()*/ { /*build the bins for all static atoms*/ selectSource(allMainAtoms, srcPat, SET1, targPat, SET2, ignorePat); abins = binAtoms(allMainAtoms, &boundingBoxA, 'a', probeRad, keepUnselected, SET1|SET2); } else /*but autobondrot mode may not yet have read in any atoms*/ {/*minimal boundingBoxes setup*/ boundingBoxA.min.x = boundingBoxA.min.y = boundingBoxA.min.z = 0.0; boundingBoxA.max.x = boundingBoxA.max.y = boundingBoxA.max.z = 0.1; abins = initBins('a', &boundingBoxA, 1); /* dummy bbox */ } if (mabis.inf == NULL) /*pseudo modal flag for autobondrot*/ {/*there is NOT an autobondrot atom info input file */ if (! ImplicitH) { /*preliminary step acts on all atoms*/ updateHydrogenInfo(outf, allMainAtoms, abins, NULL, NULL, SET1|SET2, FALSE); /*(returned param used on separate call for autobondrot processing)*/ /*e.g. creates: newH->elem = atomHOd , dummy Hydrogen to make an H-Bond*/ } /* Dump the atom information if we've been asked to. */ if (dumpFileName) { dumpFile = fopen(dumpFileName, "wb"); for (a = allMainAtoms; a; a = a->next) { /* Don't dump Phantom Hydrogen information. */ if (a->elem != atomHOd) { fprintf(dumpFile, "%s %s %3d %-4s %c %7.3f %7.3f %7.3f %5.2f %s %s %s\n", a->r->chain, a->r->resname, a->r->resid, a->atomname, a->altConf == ' ' ? '-' : a->altConf, a->loc.x, a->loc.y, a->loc.z, a->radius, a->props & ACCEPTOR_PROP ? "isAcceptor" : "noAcceptor", a->props & DONOR_PROP ? "isDonor" : "noDonor", a->props & METAL_PROP ? "isMetallic" : "noMetallic"); } } fclose(dumpFile); } /*now do the real work: */ doCommand(outf, method, allMainAtoms, abins, NULL, NULL, dots, probeRad, density, spikelen, countDots, rawOutput, conFlag, "", 0.0, drawSpike, sayGroup, groupLabel, argc, argv, message); } else {/*setup to call autobondrot*/ xformDatabase* xdb = NULL; movingCommandInfo mcis; mcis.firstPass = TRUE; mcis.keepUnselected = keepUnselected; mcis.outf = outf; mcis.method = method; mcis.allMainAtoms = allMainAtoms; mcis.waterClones= NULL; mcis.abins = abins; mcis.dots = dots; mcis.probeRad = probeRad; mcis.density = density; mcis.spikelen = spikelen; mcis.countDots = countDots; mcis.rawOutput = rawOutput; mcis.drawSpike = drawSpike; mcis.sayGroup = sayGroup; mcis.groupLabel = groupLabel; mcis.argc = argc; mcis.argv = argv; mcis.message = message; mabis.srcPat = srcPat; mabis.targPat = targPat; #define RAW_HEADER_COMMENT "#" descrCommand(outf, RAW_HEADER_COMMENT, RAW_HEADER_COMMENT, argc, argv); /*now actually read in both the mobile atoms and transformation info*/ xdb = readTransformationDatabase(mabis.inf, outf, newMovingAtom, &mabis, movingAtomListProcessing, NULL, RAW_HEADER_COMMENT); /* mabis.inf == FILE *inf */ /* outf == FILE *outf */ /* newMovingAtom == abrMkAtomProc mkAtom */ /* &mabis == void *atomstuff */ /* movingAtomListProcessing == abrAtomListProc inputListProc */ /* NULL == void *liststuff */ /* RAW_HEADER_COMMENT == char *cch */ /*autobondrot/readTransformationDatabase() */ /* also calls autobondrot/describeXformDB() */ /* which writes the header-comments to the .map file! */ if (mabis.close) { fclose(mabis.inf); } autobondrot(stderr, xdb, movingDoCommand, &mcis, deleteMovingAtom, &mabis, Verbose); /*autobondrot.c/autobondrot() is the call to do the autobondrot stuff*/ /*movingDoCommand is the name of the probeProc() called from there */ /*probe.c/movingDoCommand() in turn calls probe.c/doCommand()*/ discardXformDB(xdb, deleteMovingAtom, &mabis); }/*setup to call autobondrot*/ /*after finish, then drop through to here to release memory, etc.*/ disposeBins(abins); abins = NULL; freeDotSphere(&COdots); unloadDotSpheres(dots); freePattern(ignorePat); ignorePat = NULL; freePattern(targPat); targPat = NULL; freePattern(srcPat); srcPat = NULL; freePattern(sourcePat); sourcePat = NULL; /*111023dcr*/ freePattern(causePat); causePat = NULL; /*111023dcr*/ freePattern(extraPat); extraPat = NULL; /*111023dcr*/ if (Verbose) { note("If you publish work which uses probe, please cite:"); note(referenceString); sprintf(message, "For more information see %s", electronicReference); note(message); } }/*there are atoms for probe to work on...*/ disposeListOfResidues(resLst); resLst = NULL; disposeListOfAtoms(allMainAtoms); allMainAtoms = NULL; return 0; /*to probe.c/main() */ }/*mainProbeProc()*/ /*}}}mainProbeProc()_________________________________________________________*/ /*{{{newMovingAtom()**********************************************************/ atom * newMovingAtom(char *rec, void *userdata) { atom* a = NULL; movingAtomBuildInfo *m = (movingAtomBuildInfo *)userdata; if (rec == NULL) { return NULL; } if (m->scratchRes == NULL) { m->scratchRes = newResidueData(); } /*atom * newAtom(char *rec, int file, int model, residue * resDataBuf)*/ a = newAtom(rec, m->filenum, 0, m->scratchRes); if (a) { a->next = NULL; if (ImplicitH && isHatom(a->elem)) { deleteAtom(a); /* filter out implicit hydrogens */ return NULL; } selectSource(a, m->srcPat, SET1, m->targPat, SET2, NULL); /* note: new moving atom can not be ignored using -IGNORE flag */ if (*(m->reslstptr) == NULL) {/* first residue goes on residue list */ *(m->reslstptr) = m->scratchRes; m->scratchRes = NULL; } else { /* otherwise we have to compare residue data blocks */ if (resDiffersFromPrev(*(m->reslstptr), m->scratchRes)) { m->scratchRes->nextRes = *(m->reslstptr); *(m->reslstptr) = m->scratchRes; m->scratchRes = NULL; } else { a->r = *(m->reslstptr); /* same, so point to prior block */ a->r->a = a; /* makes this atom the head of atom list for this res */ } } } return a; } /*}}}newMovingAtom()_________________________________________________________*/ /*{{{deleteMovingAtom()*******************************************************/ void deleteMovingAtom(atom *a, void *userdata) { movingAtomBuildInfo *m = (movingAtomBuildInfo *)userdata; if (m && m->scratchRes) { deleteResidueData(m->scratchRes); /* clean up extra memory */ m->scratchRes = NULL; } deleteAtom(a); } /*}}}deleteMovingAtom()______________________________________________________*/ /*{{{movingDoCommand()********************************************************/ void movingDoCommand(char* orientationName, double scoreBias, atom *allMovingAtoms, void *userdata) {/*movingDoCommand() for autobondrot*/ /*orientationName holds the angle values as a character string*/ /* NOTE: allMainAtoms/abins & allMovingAtoms/bbins must be disjoint */ /* sets of atoms (none in common) or allMovingAtoms/bbins can */ /* be NULL. */ /* allMovingAtoms refers to autobondrot set of atoms */ movingCommandInfo *m = (movingCommandInfo *)userdata; region boundingBoxB, nonMovingBB; atom *a = NULL; atomBins *bbins = NULL; /*autobondrot.c: alst == atom *allMovingAtoms existance implies autobondrot*/ if (allMovingAtoms) { boundingBoxB.min = allMovingAtoms->loc; boundingBoxB.max = allMovingAtoms->loc; for(a = allMovingAtoms; a; a = a->next) { updateBoundingBox(&(a->loc), &boundingBoxB); } /* Expand bounding box because of extra water hydrogens */ /* added in updateHydrogenInfo() */ nonMovingBB.min = m->abins->min; nonMovingBB.max = m->abins->max; addBBox2BBox(&nonMovingBB, &boundingBoxB); bbins = binAtoms(allMovingAtoms, &boundingBoxB, 'b', m->probeRad, m->keepUnselected, SET1|SET2); if (m->abins == NULL || bbins == NULL) { halt("no atoms for autobondrot processing"); } if (! ImplicitH) { if (m->firstPass) { m->waterClones = updateHydrogenInfo(m->outf, m->allMainAtoms, m->abins, allMovingAtoms, bbins, SET1|SET2, TRUE); m->firstPass = FALSE; } else { /*the DONOR prop tells updateHydrogenInfo() to build H? atoms*/ for(a = m->waterClones; a; a = a->next) { if ((a->props & AMBIGWATER_PROP) && (a->elem == atomO)) { a->props |= DONOR_PROP; /* turn water donor property back on*/ } } /* since each pass through updateHydrogenInfo turns it off */ for(a = allMovingAtoms; a; a = a->next) { if ((a->props & AMBIGWATER_PROP) && (a->elem == atomO)) { a->props |= DONOR_PROP; /* turn water donor property back on*/ } } /* since each pass through updateHydrogenInfo turns it off */ updateHydrogenInfo(m->outf, m->waterClones, m->abins, allMovingAtoms, bbins, SET1|SET2, FALSE); } } /*autobondrot at this stage has m->method == INTERSECTONCE == 1 */ /* or m->method == SELFINTERSECT == 3 (seems as given in probe command)*/ /* and m->countDots==1, m->rawOutput==1, m->drawSpike==1 */ /* orientationName==rawname== char* holding angle values */ doCommand(m->outf, m->method, m->allMainAtoms, m->abins, allMovingAtoms, bbins, m->dots, m->probeRad, m->density, m->spikelen, m->countDots, m->rawOutput, FALSE, orientationName, scoreBias, m->drawSpike, m->sayGroup, m->groupLabel, m->argc, m->argv, m->message); // SJ -10/07/2011 by default no condensed output for this /*allMovingAtoms is closest thing to a flag for autobondrot mode*/ disposeBins(bbins); } }/*movingDoCommand() for autobondrot*/ /*}}}movingDoCommand()_______________________________________________________*/ /*{{{doCommand()******* called from mainProbeProc() ***calls writeOutput()****/ void doCommand(FILE *outf, int method, atom *allMainAtoms, atomBins *abins, atom *allMovingAtoms, atomBins *bbins, pointSet dots[], float probeRad, float density, float spikelen, int countDots, int rawOutput, int conFlag, char* rawname, double scoreBias, int drawSpike, int sayGroup, char* groupLabel, int argc, char **argv, char message[]) { /* NOTE: allMainAtoms/abins & allMovingAtoms/bbins must be disjoint */ /* sets of atoms (none in common) or allMovingAtoms/bbins can */ /* be NULL. */ /* allMovingAtoms refers to autobondrot set of atoms */ /*doCommand is called from mainProbeProc() with rawname=="" */ /* and called from movingDoCommand() with rawname==orientationName */ /* which holds the string of the current autobondrot cycle's angle values*/ /*for autobondrot: method is set by probe command, e.g. */ /* INTERSECTONCE==1 as used in example of mobile sidechain in a protein */ /* SELFINTERSECT==3 as used by alaphitaupsi ramachandran calculation*/ /* autobondrot option defines: countDots==1, rawOutput==1*/ /* autobondrot seems to get here with: drawSpike==1*/ dotNode *rslts[NUMATOMTYPES][NODEWIDTH]; int nsel = 0, numSkinDots = 0; int usesMovingAtoms = FALSE; int j=0; char extrastr[32]; /*060129 for extra master to control orig vs fitted dots*/ /*allMovingAtoms is closest thing to a flag for autobondrot mode*/ /*so if it exists and bbins exist (dcr?: mobile atoms near enough?) */ /*then usesMovingAtoms becomes the flag for autobondrot mode */ /*in effect probe is modal but too object oriented to admit such globalness*/ usesMovingAtoms = ((allMovingAtoms != NULL) && (bbins != NULL)); if (! usesMovingAtoms) { allMovingAtoms = NULL; bbins = NULL; } initResults(rslts); if (!countDots && rawOutput && Verbose) { /*write rawOutput col headers*/ if (OldStyleU) { if(conFlag) note(">>name:pat:type:srcAtom:targAtom:dot-count:min-gap:gap:" "kissEdge2BullsEye:dot2BE:dot2SC:spike:score:" "stype:ttype:x:y:z:sBval:tBval"); else note(">>name:pat:type:srcAtom:targAtom:min-gap:gap:" "kissEdge2BullsEye:dot2BE:dot2SC:spike:score:" "stype:ttype:x:y:z:sBval:tBval"); } else { if(conFlag) note(">>name:pat:type:srcAtom:targAtom:dot-count:min-gap:gap:" "spX:spY:spZ:spikeLen:score:" "stype:ttype:x:y:z:sBval:tBval"); else note(">>name:pat:type:srcAtom:targAtom:min-gap:gap:" "spX:spY:spZ:spikeLen:score:" "stype:ttype:x:y:z:sBval:tBval"); } }/*write rawOutput col headers*/ if (method == SELFINTERSECT) {/*{{{(method == SELFINTERSECT)___expected default method___****/ if (Verbose && !(countDots && rawOutput)) note("SelfIntersect"); /*SELFINTERSECT case, using one input file (and where srcPat == targPat)*/ /*is where unwashed NMR files with multiple models can be processed */ /*when there is not a single model specified in the pattern.*/ /*This uses a loop through all the model numbers actually in the file.*/ if(modelSrc == 0 && modelTarg == 0 && modelCount > 1) /*041114*/ /*multiple models in file but no model specified in pattern*/ modelLimit = modelCount; else modelLimit = 1; for(j = 1; j <= modelLimit; j++) /*041114*/ {/*loop over models*/ if(modelLimit > 1) {/*defines the multiple model case*/ modelToProcess = modelNumber[j]; /*models can be in any order 041114*/ if(Verbose) fprintf(stderr,"processing modelNumber== %d\n",modelNumber[j]); } if(j > 1) { freeResults(rslts); initResults(rslts); } genDotIntersect(allMainAtoms, abins, allMovingAtoms, bbins,dots, probeRad, spikelen,SET1, SET1, rslts); /*does this for all atoms... calls examineDots()*/ if (countDots) {/*autobondrot sets this, commandline can set this*/ if (!rawOutput) {descrCommand(outf, "program:", "command:", argc, argv);} numSkinDots = enumDotSkin(allMainAtoms, abins, allMovingAtoms, bbins, dots, SET1); /*numSkinDots used to normalize output score*/ if (!rawOutput) { fprintf(outf, "selection: self\nname: %s\n", groupLabel?groupLabel:"dots"); fprintf(outf, "density: %.1f dots per A^2\nprobeRad: %.3f A\nVDWrad: (r * %.3f) + %.3f A\n", density, probeRad, RadScaleFactor, RadScaleOffset); fprintf(outf, "score weights: gapWt=%g, bumpWt=%g, HBWt=%g\n", GAPweight, BUMPweight, HBweight); } if (usesMovingAtoms) nsel = countSelected(allMainAtoms, SET1) + countSelected(allMovingAtoms, SET1); else nsel = countSelected(allMainAtoms, SET1); if (rawOutput) {/*autobondrot sets this*/ rawEnumerate(outf, "", rslts, method, nsel, drawSpike, FALSE, numSkinDots, density, groupLabel?groupLabel:"", rawname, scoreBias); } else { enumerate(outf, "self dots", rslts, probeRad, method, nsel, drawSpike, FALSE, numSkinDots, density); } }/*countDots*/ else { if (rawOutput) { writeRaw(outf, "1->1", rslts, probeRad,groupLabel?groupLabel:"", density,conFlag); } else if (OutputFmtType == 1) writeAltFmtO(outf, TRUE, TRUE, "self_dots", rslts, drawSpike); else if (OutputFmtType == 2) { descrCommand(outf, "# software:", "# command:", argc, argv); writeAltFmtXV(outf, TRUE, TRUE, "self_dots", rslts, drawSpike); } else if (OutputFmtType == 3) { /*dcr041101 ONELINE :count:summaries:*/ countsummary(outf,"SelfIntersect", 1, 1); /*dcr041101 Lines,Pass*/ } else {/*kinemage output*/ if(ContactSUMMARY) /*dcr041101 Lines,Pass*/ countsummary(outf,"SelfIntersect", 9, 1); descrCommand(outf, "@caption", " command:", argc, argv); if (sayGroup) { if(modelLimit > 1) {/*doing jth of multiple models of an ensemble*/ fprintf(outf, "@group dominant {%s M%d} animate\n",groupLabel?groupLabel:"dots",j); } else fprintf(outf, "@group dominant {%s}\n",groupLabel?groupLabel:"dots"); } sprintf(extrastr,"%s",groupLabel?groupLabel:"dots"); /*060129*/ writeOutput(outf, "self dots", rslts, drawSpike, method, extrastr, probeRad); /*add probeRad 20111220dcr*/ }/*kinemage output*/ } }/*loop over models*/ }/*}}}(method == SELFINTERSECT)________________________________*/ else if (method == INTERSECTONCE) {/*{{{(method == INTERSECTONCE)*********************************/ if (Verbose && !(countDots && rawOutput)) note("IntersectOnce"); genDotIntersect(allMainAtoms, abins, allMovingAtoms, bbins,dots, probeRad, spikelen, SET1, SET2, rslts); if (countDots) {/*autobondrot sets this*/ if (!rawOutput) descrCommand(outf, "program:", "command:", argc, argv); numSkinDots = enumDotSkin(allMainAtoms, abins, allMovingAtoms, bbins, dots, SET1); /*numSkinDots used to normalize output score*/ if (!rawOutput) { fprintf(outf, "selection: once\nname: %s\n", groupLabel?groupLabel:"dots"); fprintf(outf, "density: %.1f dots per A^2\nprobeRad: %.3f A\nVDWrad: (r * %.3f) + %.3f A\n", density, probeRad, RadScaleFactor, RadScaleOffset); fprintf(outf, "score weights: gapWt=%g, bumpWt=%g, HBWt=%g\n", GAPweight, BUMPweight, HBweight); } if (usesMovingAtoms) nsel = countSelected(allMainAtoms, SET1) + countSelected(allMovingAtoms, SET1); else nsel = countSelected(allMainAtoms, SET1); if (rawOutput) {/*autobondrot sets this*/ rawEnumerate(outf,"", rslts, method, nsel, drawSpike, FALSE, numSkinDots, density, groupLabel?groupLabel:"", rawname, scoreBias); } else { enumerate(outf,"once dots", rslts, probeRad, method, nsel, drawSpike, FALSE, numSkinDots, density); } } else { if (rawOutput) { writeRaw(outf, "1->2", rslts, probeRad,groupLabel?groupLabel:"", density,conFlag); } else if (OutputFmtType == 1) writeAltFmtO(outf, TRUE, TRUE, "once_dots", rslts, drawSpike); else if (OutputFmtType == 2) { descrCommand(outf, "# software:", "# command:", argc, argv); writeAltFmtXV(outf, TRUE, TRUE, "once_dots", rslts, drawSpike); } else if (OutputFmtType == 3) { /*dcr041101 ONELINE :count:summaries:*/ countsummary(outf,"IntersectOnce", 1, 1); /*dcr041101 Lines,Pass*/ } else {/*write kinemage*/ if(ContactSUMMARY) /*dcr041101 Lines,Pass*/ countsummary(outf,"IntersectOnce", 9, 1); descrCommand(outf, "@caption", " command:", argc, argv); if (sayGroup) fprintf(outf, "@group dominant {%s}\n", groupLabel?groupLabel:"dots"); sprintf(extrastr,"%s",groupLabel?groupLabel:"dots"); /*060129*/ writeOutput(outf, "once dots", rslts, drawSpike, method, extrastr, probeRad); /*add probeRad 20111220dcr*/ } } }/*}}}(method == INTERSECTONCE)________________________________*/ else if (method == INTERSECTBOTHWAYS) {/*{{{(method == INTERSECTBOTHWAYS)*****************************/ if (Verbose && !(countDots && rawOutput)) note("IntersectBothWays"); if (countDots) { if (!rawOutput) { descrCommand(outf, "program:", "command:", argc, argv); fprintf(outf, "selection: both\nname: %s\n", groupLabel?groupLabel:"dots"); fprintf(outf, "density: %.1f dots per A^2\nprobeRad: %.3f A\nVDWrad: (r * %.3f) + %.3f A\n", density, probeRad, RadScaleFactor, RadScaleOffset); fprintf(outf, "score weights: gapWt=%g, bumpWt=%g, HBWt=%g\n", GAPweight, BUMPweight, HBweight); } } else { if (rawOutput) {} /* do nothing on purpose */ else if (OutputFmtType == 1) {} else if (OutputFmtType == 2) descrCommand(outf, "# software:", "# command:", argc, argv); else {/*kinemage: keywords before double pass*/ descrCommand(outf, "@caption", " command:", argc, argv); if (sayGroup) fprintf(outf, "@group {%s}\n",groupLabel?groupLabel:"dots"); } } genDotIntersect(allMainAtoms, abins, allMovingAtoms, bbins, dots, probeRad, spikelen, SET1, SET2, rslts); if (countDots) { numSkinDots = enumDotSkin(allMainAtoms, abins, allMovingAtoms, bbins, dots, SET1); /*numSkinDots used to normalize output score*/ if (usesMovingAtoms) nsel = countSelected(allMainAtoms, SET1) + countSelected(allMovingAtoms, SET1); else nsel = countSelected(allMainAtoms, SET1); if (rawOutput) { rawEnumerate(outf, "1->2", rslts, method, nsel, drawSpike, FALSE, numSkinDots, density, groupLabel?groupLabel:"", rawname, scoreBias); } else { enumerate(outf, "1->2", rslts, probeRad, method, nsel, drawSpike, FALSE, numSkinDots, density); } } else { if (rawOutput) { writeRaw(outf, "1->2", rslts, probeRad, groupLabel?groupLabel:"", density,conFlag); } else if (OutputFmtType == 1) writeAltFmtO(outf, TRUE, !sayGroup, "1->2", rslts, drawSpike); else if (OutputFmtType == 2) writeAltFmtXV(outf, TRUE, !sayGroup, "1->2", rslts, drawSpike); else if (OutputFmtType == 3) { /*dcr041101 ONELINE :count:summaries:*/ countsummary(outf,"IntersectBothWays 1->2", 1, 0); /*dcr041101 Lines,Pass==0 no output for this pass*/ } else { sprintf(extrastr,"%s",groupLabel?groupLabel:"dots"); /*060129*/ writeOutput(outf, "1->2", rslts, drawSpike, method, extrastr,probeRad); /*add probeRad 20111220dcr*/ if(ContactSUMMARY) /*dcr041101 Lines,Pass*/ countsummary(outf,"IntersectBothWays 1->2", 9, 0); /*dcr041101 Lines,Pass==0 no output for this pass*/ } } freeResults(rslts); initResults(rslts); genDotIntersect(allMainAtoms, abins, allMovingAtoms, bbins, dots, probeRad, spikelen, SET2, SET1, rslts); if (countDots) { numSkinDots = enumDotSkin(allMainAtoms, abins, allMovingAtoms, bbins, dots, SET2); /*numSkinDots used to normalize output score*/ if (usesMovingAtoms) nsel = countSelected(allMainAtoms, SET2) + countSelected(allMovingAtoms, SET2); else nsel = countSelected(allMainAtoms, SET2); if (rawOutput) { rawEnumerate(outf, "2->1", rslts, method, nsel, drawSpike, FALSE, numSkinDots, density, groupLabel?groupLabel:"", rawname, scoreBias); } else { enumerate(outf, "2->1", rslts, probeRad, method, nsel, drawSpike, FALSE, numSkinDots, density); } } else { if (rawOutput) { writeRaw(outf, "2->1", rslts, probeRad, groupLabel?groupLabel:"", density,conFlag); } else if (OutputFmtType == 1) { writeAltFmtO(outf, !sayGroup, TRUE, "2->1", rslts, drawSpike); } else if (OutputFmtType == 2) { writeAltFmtXV(outf, !sayGroup, TRUE, "2->1", rslts, drawSpike); } else if (OutputFmtType == 3) { /*dcr041101 ONELINE :count:summaries:*/ countsummary(outf,"IntersectBothWays 2->1", 1, 2); /*dcr041101 Lines,Pass==2 output sum of 2 passes*/ } else { sprintf(extrastr,"%s",groupLabel?groupLabel:"dots"); /*060129*/ writeOutput(outf, "2->1", rslts, drawSpike, method, extrastr, probeRad); /*add probeRad 20111220dcr*/ if(ContactSUMMARY) /*dcr041101 Lines,Pass*/ countsummary(outf,"IntersectBothWays 2->1", 9, 2); /*dcr041101 Lines,Pass==2 output sum of 2 passes*/ } } }/*}}}(method == INTERSECTBOTHWAYS)____________________________*/ else if (method == EXTERNALSURFACE) {/*{{{(method == EXTERNALSURFACE)*******************************/ if (Verbose && !(countDots && rawOutput)) note("ExternalSurface"); genDotSurface(allMainAtoms, abins, allMovingAtoms, bbins, dots, probeRad, spikelen, SET1, rslts); if (countDots) { if (!rawOutput) descrCommand(outf, "program:", "command:", argc, argv); numSkinDots = enumDotSkin(allMainAtoms, abins, allMovingAtoms, bbins, dots, SET1); /*numSkinDots used to normalize output score*/ if (!rawOutput) { fprintf(outf, "selection: external\nname: %s\n",groupLabel?groupLabel:"dots"); fprintf(outf, "density: %.1f dots per A^2\nprobeRad: %.3f A\nVDWrad: (r * %.3f) + %.3f A\n", density, probeRad, RadScaleFactor, RadScaleOffset); } if (usesMovingAtoms) nsel = countSelected(allMainAtoms, SET1) + countSelected(allMovingAtoms, SET1); else nsel = countSelected(allMainAtoms, SET1); if (rawOutput) { rawEnumerate(outf, "", rslts, method, nsel, FALSE, TRUE, numSkinDots, density, groupLabel?groupLabel:"", rawname, scoreBias); } else { enumerate(outf, "extern dots", rslts, probeRad, method, nsel, FALSE, TRUE, numSkinDots, density); } } else { if (rawOutput) { writeRaw(outf, "1->none", rslts, probeRad, groupLabel?groupLabel:"", density,conFlag); } else if (OutputFmtType == 1) { writeAltFmtO(outf, TRUE, TRUE, "extern_dots", rslts, FALSE); } else if (OutputFmtType == 2) { descrCommand(outf, "# software:", "# command:", argc, argv); writeAltFmtXV(outf, TRUE, TRUE, "extern_dots", rslts, FALSE); } else { descrCommand(outf, "@caption", " command:", argc, argv); if (sayGroup) fprintf(outf, "@group dominant {%s}\n", groupLabel?groupLabel:"dots"); sprintf(extrastr,"%s",groupLabel?groupLabel:"dots"); /*060129*/ writeOutput(outf, "extern dots", rslts, FALSE, method, extrastr, probeRad); /*add probeRad 20111220dcr*/ } } }/*}}}(method == EXTERNALSURFACE)______________________________*/ else if (method == DUMPATOMCOUNT) {/*{{{(method == DUMPATOMCOUNT)*********************************/ if (Verbose && !rawOutput) { note("dumpAtomInfo"); descrCommand(outf, "program:", "command:", argc, argv); fprintf(outf, "selection: self\nname: %s\n", groupLabel?groupLabel:"dots"); fprintf(outf, "density: %.1f dots per A^2\nprobeRad: %.3f A\nVDWrad: (r * %.3f) + %.3f A\n", density, probeRad, RadScaleFactor, RadScaleOffset); fprintf(outf, "score weights: gapWt=%g, bumpWt=%g, HBWt=%g\n", GAPweight, BUMPweight, HBweight); } if (usesMovingAtoms) nsel = countSelected(allMainAtoms, SET1) + countSelected(allMovingAtoms, SET1); else nsel = countSelected(allMainAtoms, SET1); if (rawOutput) { if (groupLabel) { fprintf(outf, "%d %s %s%s\n", nsel, rawname, RAW_HEADER_COMMENT, groupLabel); } else { fprintf(outf, "%d %s\n", nsel, rawname); } } else { fprintf(outf, "atoms selected: %d\n", nsel); } }/*}}}(method == DUMPATOMCOUNT)_______________________________*/ freeResults(rslts); }/*doCommand()*/ /*}}}doCommand()_____________________________________________________________*/ /*{{{descrCommand()***********************************************************/ void descrCommand(FILE *outf, char *hdr1, char *hdr2, int argc, char **argv) { int i = 0; time_t t; char *ts = NULL; t = time(NULL); ts = asctime(localtime(&t)); /* system dependent time string */ fprintf(outf, "%s %s, run %s", hdr1, shortVersionStr, ts); fprintf(outf, "%s %s", hdr2, argv[0]); for (i = 1; i < argc; i++) { int dontQuote = (strpbrk(argv[i], " \t\n\b\r\v\f\a*\'") == NULL); fprintf(outf,(dontQuote?" %s":" \"%s\""), argv[i]); } fprintf(outf, "\n"); } /*}}}descrCommand()__________________________________________________________*/ /*{{{loadDotSpheres()*********************************************************/ void loadDotSpheres(pointSet dots[], float density) { int i = 0; for (i = 0; i < NUMATOMTYPES; i++) { dotSphere(&(dots[i]), getRadius(i, 0), density); } dotSphere(&COdots, getRadius(atomC, 1), density); } /*}}}loadDotSpheres()________________________________________________________*/ /*{{{unloadDotSpheres()*******************************************************/ void unloadDotSpheres(pointSet dots[]) { int i = 0; for (i = 0; i < NUMATOMTYPES; i++) { freeDotSphere(&(dots[i])); } } /*}}}unloadDotSpheres()______________________________________________________*/ /*{{{probehelp()**************************************************************/ void probehelp(int longlist) { fprintf(stderr, "\nSyntax: probe input.pdb >> out.kin\n"); fprintf(stderr, " or: probe [flags] \"src pattern\" [\"target pattern\"] pdbfiles... >> out.kin\n\n"); fprintf(stderr, "Flags:\n"); fprintf(stderr, " -SElf self intersection: src -> src (default)\n"); fprintf(stderr, " -Both intersect both ways: src <=> targ\n"); fprintf(stderr, " -ONce single intersection: src -> targ\n"); fprintf(stderr, " -OUt external van der Waals surface of src (solvent contact surface)\n"); fprintf(stderr, "\n"); fprintf(stderr, " -AUTObondrot filename read and process an autobondrot file\n"); if (! longlist) { fprintf(stderr, "\n shortcuts: <>, -SCSurface, -EXPOsed, -ASurface, -ACCESS, -SCAN0, -SCAN1\n"); fprintf(stderr, "\n most simple dotkin: probe -self all -kinemage input.pdb > out.kin\n"); } else { /*longlist option*/ fprintf(stderr, "\n shortcuts:\n"); fprintf(stderr, " <>same as: -4H -mc -het -self \"altA ogt33\"\n"); /*changed from 3 : dcr041017*/ fprintf(stderr, " -DEFAULTs same as: <>, but allows some other flags\n"); /*dcr041101*/ fprintf(stderr, " -SCSurface same as: -drop -rad1.4 -out \"not water\"\n"); fprintf(stderr, " -EXPOsed same as: -drop -rad1.4 -out (note: user supplies pattern)\n"); fprintf(stderr, " -ASurface same as: -drop -rad0.0 -add1.4 -out \"not water\"\n"); fprintf(stderr, " -ACCESS same as: -drop -rad0.0 -add1.4 -out (note: user supplies pattern)\n"); fprintf(stderr, " -SCAN0 same as: -4H -mc -self \"alta blt40 ogt33\"\n"); /*changed from 3 : dcr041017*/ fprintf(stderr, " -SCAN1 same as: -4H -once \"sc alta blt40 ogt33\" \"alta blt40 ogt65,(not water ogt33)\"\n"); /*changed from 3 : dcr041017*/ fprintf(stderr, "\n"); fprintf(stderr, " -DUMPAtominfo count the atoms in the selection: src\n\n"); fprintf(stderr, " (note that BOTH and ONCE require two patterns while\n"); fprintf(stderr, " OUT, SELF and DUMPATOMINFO require just one pattern)\n\n"); fprintf(stderr, " -Implicit implicit hydrogens\n"); fprintf(stderr, " -Explicit explicit hydrogens (default)\n"); fprintf(stderr, " -DEnsity# set dot density (default 16 dots/sq A)\n"); fprintf(stderr, " -Radius#.# set probe radius (default 0.25 A)\n"); fprintf(stderr, " -ADDvdw#.# offset added to Van der Waals radii (default 0.0)\n"); fprintf(stderr, " -SCALEvdw#.# scale factor for Van der Waals radii (default 1.0)\n"); fprintf(stderr, " -COSCale#.# scale C=O carbon Van der Waals radii (default 0.94)\n"); fprintf(stderr, " -SPike draw spike instead of dots (default)\n"); fprintf(stderr, " -SPike#.# set spike scale (default=0.5)\n"); fprintf(stderr, " -NOSpike draw only dots\n"); fprintf(stderr, " -HBRegular#.# max overlap for regular Hbonds(default=%.1f)\n", Min_regular_hb_cutoff); fprintf(stderr, " -HBCharged#.# max overlap for charged Hbonds(default=%.1f)\n", Min_charged_hb_cutoff); fprintf(stderr, " -Keep keep nonselected atoms (default)\n"); fprintf(stderr, " -DRop drop nonselected atoms\n"); fprintf(stderr, " -LIMit limit bump dots to max dist when kissing (default)\n"); fprintf(stderr, " -NOLIMit do not limit bump dots\n"); fprintf(stderr, " -LENs add lens keyword to kin file\n"); fprintf(stderr, " -NOLENs do not add lens keyword to kin file (default)\n"); fprintf(stderr, " -MC include mainchain->mainchain interactions\n"); fprintf(stderr, " -HETs include dots to non-water HET groups (default)\n"); fprintf(stderr, " -NOHETs exclude dots to non-water HET groups\n"); fprintf(stderr, " -WATers include dots to water (default)\n"); fprintf(stderr, " -NOWATers exclude dots to water\n"); fprintf(stderr, " -WAT2wat show dots between waters\n"); fprintf(stderr, " -DUMPH2O include water H? vectorlist in output\n"); fprintf(stderr, " -4H extend bond chain dot removal to 4 for H (default)\n"); fprintf(stderr, " -3 limit bond chain dot removal to 3\n"); fprintf(stderr, " -2 limit bond chain dot removal to 2\n"); fprintf(stderr, " -1 limit bond chain dot removal to 1\n"); fprintf(stderr, " -IGNORE \"pattern\" explicit drop: ignore atoms selected by pattern\n"); fprintf(stderr, " -DOCHO recognize CH..O Hbonds\n"); fprintf(stderr, " -CHO#.# scale factor for CH..O Hbond score (default=%.1f)\n", CHOHBfactor); fprintf(stderr, " -PolarH use short radii of polar hydrogens (default)\n"); fprintf(stderr, " -NOPolarH do not shorten radii of polar hydrogens\n"); fprintf(stderr, " -NUClear use nuclear position vdW radii\n"); // 130114 JJH fprintf(stderr, " (default is electron cloud positions)\n"); // 130114 JJH fprintf(stderr, " -NOFACEhbond do not identify HBonds to aromatic faces\n\n"); fprintf(stderr, " -Name \"name\" specify the group name (default \"dots\")\n"); fprintf(stderr, " -DOTMASTER group name used as extra master={name} on lists\n"); /*060129*/ fprintf(stderr, " -NOGroup do not generate @group statement in .kin format output\n"); fprintf(stderr, " -KINemage add @kinemage 1 statement to top of .kin format output\n\n"); fprintf(stderr, " -Countdots produce a count of dots-not a dotlist\n"); fprintf(stderr, " -Unformated output raw dot info\n"); fprintf(stderr, " name:pat:type:srcAtom:targAtom:mingap:gap:spX:spY:spZ:spikeLen:score:stype:ttype:x:y:z:sBval:tBval:\n"); fprintf(stderr, " -CONdense raw output in condensed format, i.e. one line per source->target atom. Also gives the dot count for that interaction. Works only with -Unformated flag\n"); fprintf(stderr, " name:pat:type:srcAtom:targAtom:dot-count:mingap:gap:spX:spY:spZ:spikeLen:score:stype:ttype:x:y:z:sBval:tBval:\n"); fprintf(stderr, " -OFORMAT output dot info formatted for display in O\n"); fprintf(stderr, " -XVFORMAT output dot info formatted for display in XtalView\n"); fprintf(stderr, " -ONELINE output one line :contacts:by:severity:type:\n"); /*dcr041101*/ fprintf(stderr, " -ONEDOTeach output one dot for each src-to-neighbor\n"); /*dcr111013*/ fprintf(stderr, " With -Unformated flag: for H gives angle parent-src-cause,\n"); fprintf(stderr, " name:pat:type:srcAtom:targAtom:XHTangle:mingap:gap:spX:spY:spZ:spikeLen:score:stype:ttype:x:y:z:sBval:tBval:\n"); fprintf(stderr, " -GAPBINs 2nd pointmaster character, -.5...0...+.5\n"); /*dcr111020*/ fprintf(stderr, " -NEAREST apply selection only on nearest neighbor\n"); /*dcr111020*/ fprintf(stderr, " -FILTERnearest apply 3rd selection on neighbor\n"); /*dcr111020*/ fprintf(stderr, " e.g. -nearest -filternearest patt1 all patt3\n"); /*dcr111020*/ fprintf(stderr, " -OCC1filter final filter exclude any w occ < 1\n"); /*dcr111020*/ fprintf(stderr, " -NOALTfilter final filter exclude any alts\n"); /*dcr111020*/ fprintf(stderr, " -GAPcolor color dots by gap amount (default)\n"); fprintf(stderr, " -ATOMcolor color dots by atom type\n"); fprintf(stderr, " -BASEcolor color dots by nucleic acid base type\n"); fprintf(stderr, " -COLORBase color dots by gap and nucleic acid base type\n"); fprintf(stderr, " -OUTCOLor \"name\" specify the point color for -OUT (default \"gray\")\n"); fprintf(stderr, " -GAPWeight# set weight for scoring gaps (default 0.25)\n"); fprintf(stderr, " -BUMPWeight# set relative scale for scoring bumps (default 10.0)\n"); fprintf(stderr, " -HBWeight# set relative scale for scoring Hbonds (default 4.0)\n"); fprintf(stderr, " -DIVLow#.# Division for Bump categories (default -0.4)\n"); fprintf(stderr, " -DIVHigh#.# Division for Contact categories (default 0.25)\n"); fprintf(stderr, " -DIVWorse#.# Division for regarding a clash as a worse clash (default -0.5)\n"); /* 04/21/2015 SJ*/ fprintf(stderr, " -SEPWORSE To separate bad overlaps and worse overlaps (default false, if true, deafult value to separate the worse clashes is -0.5)\n"); fprintf(stderr, " -MINOCCupancy#.# Occupancy below this is same as zero (default 0.02)\n"); fprintf(stderr, " -ELEMent add master buttons for different elements in kin output\n"); fprintf(stderr, " -NOHBOUT do not output contacts for HBonds\n"); fprintf(stderr, " -NOCLASHOUT do not output contacts for clashes\n"); fprintf(stderr, " -NOVDWOUT do not output contacts for van der Waals interactions\n"); fprintf(stderr, " -ONLYBADOUT onlybadout output bad clashes (severe overlap contacts)\n"); /*dcr041010*/ fprintf(stderr, " -SUMMARY output summary list of contacts and clashes\n"); /*dcr041101*/ fprintf(stderr, " -ONELINE output summary list on oneline\n"); /*dcr041101*/ fprintf(stderr, " -NOTICKs do not display the residue name ticker during processing\n"); fprintf(stderr, " -STDBONDs assume only standard bonding patterns in standard residues\n"); fprintf(stderr, " -NOPARENT do not bond hydrogens based on table of parent heavy atoms\n\n"); fprintf(stderr, " -SEGID use the PDB SegID field to descriminate between residues\n"); fprintf(stderr, " -OLDU generate old style -u output: kissEdge2BullsEye, etc\n"); fprintf(stderr, " -VErbose verbose mode (default)\n"); fprintf(stderr, " -DOTDUMP dump dot info while doing examineDots()\n"); fprintf(stderr, " -DUMPATOMS filename dump atom information to filename for regression testing\n"); fprintf(stderr, " -REFerence display reference string\n"); fprintf(stderr, " -CHANGEs display a list of program changes\n"); fprintf(stderr, " -Quiet quiet mode\n"); }/*longlist option*/ fprintf(stderr, "\n -Help show expanded help notice (includes other flags)\n"); fprintf(stderr, "\n -VERSION one line version to stdout\n\n");/*dcr041009*/ fprintf(stderr, "Pattern elements: (should be put in quotes on the command line)\n"); if (longlist) { /*longlist option continues*/ fprintf(stderr, " FILE# within file #\n"); fprintf(stderr, " MODEL# within model #\n"); fprintf(stderr, " CHAINaa within chain a\n"); fprintf(stderr, " SEGaaaa segment identifier aaaa (where _ represents blank)\n"); fprintf(stderr, " ALTa alternate conformation a\n"); fprintf(stderr, " ATOMaaaa atom name aaaa (where _ represents blank)\n"); fprintf(stderr, " (all 4 characters are used so H would be ATOM_H__)\n"); fprintf(stderr, " RESaaa residue aaa\n"); fprintf(stderr, " # residue #\n"); fprintf(stderr, " #a residue #, insert a\n"); fprintf(stderr, " #-# residue range # (insert codes ignored)\n"); fprintf(stderr, " a residue type by one letter codes (eg. y)\n"); fprintf(stderr, " aaa residue type by three letter codes (eg. tyr)\n"); fprintf(stderr, " ALL,PROTEIN,MC,SC,BASE,ALPHA,BETA,NITROGEN,CARBON,OXYGEN,\n"); fprintf(stderr, " SULFUR,PHOSPHORUS,HYDROGEN,METAL,POLAR,NONPOLAR,CHARGED,\n"); fprintf(stderr, " DONOR,ACCEPTOR,AROMATIC,METHYLENE,METHYL,HET,WATER,DNA,RNA\n"); fprintf(stderr, " all or a subset of the atoms\n"); fprintf(stderr, " OLT# Occupancy less than # (integer percent)\n"); fprintf(stderr, " OGT# Occupancy greater than # (integer percent)\n"); fprintf(stderr, " BLT# B-value less than # (integer)\n"); fprintf(stderr, " BGT# B-value greater than # (integer)\n"); fprintf(stderr, " INSa Insert code a (where _ represents blank)\n"); fprintf(stderr, " \n"); fprintf(stderr, " WITHIN #.# OF #.#, #.#, #.# atoms within distance from point\n"); }/*longlist option continues*/ else { /*NOT longlist*/ fprintf(stderr, " # residue number\n"); fprintf(stderr, " #a residue #, insert a\n"); fprintf(stderr, " #-# residue number range\n"); fprintf(stderr, " a OR aaa residue type by one or three letter codes\n"); fprintf(stderr, " FILE#,MODEL#,CHAINa,SEGaaaa,ALTa,ATOMaaaa,RESaaa,\n"); fprintf(stderr, " ALL,PROTEIN,MC,SC,BASE,ALPHA,BETA,NITROGEN,CARBON,OXYGEN,\n"); fprintf(stderr, " SULFUR,PHOSPHORUS,HYDROGEN,METAL,POLAR,NONPOLAR,CHARGED,\n"); fprintf(stderr, " DONOR,ACCEPTOR,AROMATIC,METHYLENE,METHYL,HET,WATER,DNA,RNA,\n"); fprintf(stderr, " OLT#, OGT#, BLT#, BGT#, INSa, WITHIN #.# OF #.#, #.#, #.#\n"); } fprintf(stderr, " \n"); fprintf(stderr, " Patterns can be combined into comma separated lists\n"); fprintf(stderr, " such as \"trp,phe,tyr\" meaning TRP or PHE or TYR.\n"); fprintf(stderr, " \n"); fprintf(stderr, " Patterns that are sepatated by blanks must all be true\n"); fprintf(stderr, " such as \"chainb 1-5\" meaning residues 1 to 5 in chain B.\n"); fprintf(stderr, " \n"); fprintf(stderr, " You can also group patterns with parenthesis, separate\n"); fprintf(stderr, " multiple patterns with | meaning 'or' and choose the\n"); fprintf(stderr, " complement with NOT as in \"not file1\" meaning not in file 1.\n"); fprintf(stderr, "\n"); fprintf(stderr, " An autobondrot file is similar to other PDB input files\n"); fprintf(stderr, " but it includes information identifying atoms subject to rotations\n"); fprintf(stderr, " and other transformations.\n"); if (longlist) { /*more longlist*/ fprintf(stderr, "\n Example autobondrot file fragment showing Calpha-Cbeta bond rotation\n"); fprintf(stderr, " and a periodic torsion penalty function for this rotation\n"); fprintf(stderr, " ATOM 1 CB TYR 61 34.219 17.937 4.659 1.00 0.00\n"); fprintf(stderr, " bondrot:chi1:78.7: 0:359:5:33.138:18.517: 5.531:34.219:17.937: 4.659\n"); fprintf(stderr, " cos:-3:60:3:\n"); fprintf(stderr, " ATOM 1 1HB TYR 61 34.766 18.777 4.206 1.00 0.00\n"); fprintf(stderr, " ATOM 1 2HB TYR 61 34.927 17.409 5.315 1.00 0.00\n"); fprintf(stderr, " ATOM 1 CG TYR 61 33.836 16.989 3.546 1.00 0.00\n"); fprintf(stderr, " ...\n"); fprintf(stderr, " Autobondrot commands use colons to separate values\n"); fprintf(stderr, " Transformations: BONDROT:id:currAng:start:end:stepSz:x1:y1:z1:x2:y2:z2\n"); fprintf(stderr, " TRANS: id:currpos:start:end:stepSz:x1:y1:z1:x2:y2:z2\n"); fprintf(stderr, " NULL # dummy\n"); fprintf(stderr, " Bias functions: COS:scale:phaseOffset:frequency\n"); fprintf(stderr, " POLY:scale:offset:polynomialDegree\n"); fprintf(stderr, " CONST:value\n"); fprintf(stderr, " Branching: SAVE and RESTORE or \"(\" and \")\"\n"); fprintf(stderr, " (e.g. to rotate each Chi and the methyls for isoleucine the\n"); fprintf(stderr, " sequence is: rotChi1/SAVE/rotChi2/rotCD1/RESTORE/rotCG2)\n"); fprintf(stderr, " Set orientation: GO:angle1:angle2:...\n"); fprintf(stderr, " Include files: @filename\n"); fprintf(stderr, " Comments: # comment text\n"); }/*more longlist*/ fprintf(stderr, "\n%s\n", versionString); exit(0); }/*probehelp()*/ /*}}}probehelp()_____________________________________________________________*/ /*{{{probeversion()***********************************************************/ void probeversion(void) /*VERSION dcr041009*/ {/*stdout one liner version for, e.g., MolProbity dcr041009*/ /*beware: orig defaults get superseded by shortcuts, like raw "probe"*/ printf("%s\n", shortVersionStr); exit(0); } /*}}}probeversion()__________________________________________________________*/ /*{{{processCommandline()**** called from mainProbeProc() ********************/ atom* processCommandline(int argc, char **argv, int *method, region *bboxA, float *density, float *probeRad, int *drawSpike, float *spikelen, int *countDots, int *keepUnselected, char **srcArg, char **targArg, char **extraArg, char **ignoreArg, char **groupLabel, int *rawOutput, int *conFlag, int *sayGroup, int *addKinToFile, movingAtomBuildInfo *mabip, residue **reslstptr) //conFlag added SJ - 10/07/2011, extraArg DCR {/*processCommandline()*/ /*gets running conditions from the commandline as well as */ /*loads atomlist from in file, returns atomlist which becomes allMainAtoms*/ /* atomlist = probe.c/loadAtoms() */ /*called with address of mabis, i.e. mabip moving atom build info ptr */ int file = 0, nargs = 0, argcnt = 0, i = 0, n = 0; char *p, message[200]; FILE *inf = stdin; atom *atomlist = NULL; *method = SELFINTERSECT; nargs = 1; *density = 16.0; *probeRad = 0.25; *drawSpike = TRUE; *spikelen = 0.50; *countDots = FALSE; *keepUnselected = TRUE; *rawOutput = FALSE; *conFlag = FALSE; // added by SJ - 10/07/2011 *sayGroup = TRUE; *addKinToFile = FALSE; *srcArg = NULL; *targArg = NULL; *ignoreArg=NULL; *groupLabel = NULL; /*&mabis of mainProbeProc()== *mabip here */ if (mabip) { mabip->filenum = 0; mabip->inf = NULL; mabip->close = FALSE; mabip->srcPat = NULL; mabip->targPat = NULL; mabip->reslstptr= reslstptr; mabip->scratchRes= NULL; } file = 0; argcnt = 0; for (i = 1; i < argc; i++) { p = argv[i]; if (p[0] == '-') { /* -flag item */ if ((argcnt >= nargs) && (p[1] == '\0')) { file++; atomlist = loadAtoms(stdin,atomlist,bboxA,file, reslstptr); } else if(compArgStr(p+1, "HELP", 1)){ probehelp(1); } else if(compArgStr(p+1, "VERSION", 7)){/*before VERBOSE, dcr041009*/ probeversion(); } else if(compArgStr(p+1, "VERBOSE", 2)){ Verbose = TRUE; } else if(compArgStr(p+1, "DOTDUMP", 7)){ Ldotdump = TRUE; /*111202dcr Ldotdump in examineOneDotEach() */ } else if(compArgStr(p+1, "DUMPATOMS", 9)){ if (++i < argc) { dumpFileName = argv[i]; } else { halt("no filename after -DUMPATOMS flag"); } } else if(compArgStr(p+1, "QUIET", 1)){ Verbose = FALSE; } else if(compArgStr(p+1, "EXPLICIT", 1)){ ImplicitH = FALSE; } else if(compArgStr(p+1, "IMPLICIT", 1)){ ImplicitH = TRUE; } else if(compArgStr(p+1, "DROP", 2)){ *keepUnselected = FALSE; } else if(compArgStr(p+1, "KEEP", 2)){ *keepUnselected = TRUE; } else if((n = compArgStr(p+1, "DENSITY", 2))){ *density = parseReal(p, n+1, 10); } else if((n = compArgStr(p+1, "RADIUS", 1))){ *probeRad = parseReal(p, n+1, 10); /* parameters correlated with radius */ if ((*probeRad) > GAPweight) { GAPweight = *probeRad; } if ((*probeRad) > HighGoodCut) { HighGoodCut = *probeRad; } } else if((n = compArgStr(p+1, "SCALEVDW", 5)) /* new name,*/ || (n = compArgStr(p+1, "VDWSCALE", 3))){ /* old name */ RadScaleFactor = parseReal(p, n+1, 10); if (RadScaleFactor < 0.0001 || RadScaleFactor > 1000.0) { sprintf(message, "invalid scale factor: -scalevdw%g", RadScaleFactor); halt(message); } } else if((n = compArgStr(p+1, "ADDVDW", 3))){ RadScaleOffset = parseReal(p, n+1, 10); if (RadScaleOffset < -10.0 || RadScaleOffset > 1000.0) { sprintf(message, "invalid vdw offset: -addvdw%g", RadScaleOffset); halt(message); } } else if((n = compArgStr(p+1, "SPIKE", 2))){ *drawSpike = TRUE; if (p[n+1]) { *spikelen = parseReal(p, n+1, 10); } } else if(compArgStr(p+1, "NOSPIKE", 3)){ *drawSpike = FALSE; } else if(compArgStr(p+1, "KINEMAGE", 3)){ *addKinToFile = TRUE; *countDots = FALSE; /* also forces kin format */ *rawOutput = FALSE; *conFlag = FALSE; // by SJ 10/07/2011 OutputFmtType = 0; } else if(compArgStr(p+1, "NOGROUP", 3)){ *sayGroup = FALSE; } else if(compArgStr(p+1, "ELEMENT", 4)){ AtomMasters = TRUE; } else if(compArgStr(p+1, "COUNTDOTS", 1)){ *countDots = TRUE; } else if(compArgStr(p+1, "UNFORMATED", 1)){ *rawOutput = TRUE; } else if(compArgStr(p+1, "CONDENSE", 3)){ *conFlag = TRUE;// added by SJ 10/07/2011 - to parse the flag for condened raw output } else if(compArgStr(p+1, "WEAKHbonds", 5)){ LweakHbonds = TRUE; /*111215dcr*/ } else if(compArgStr(p+1,"SEPWORSE",8)){ LworseOverlap = TRUE; /*04/09/2015 SJ - for separating the worse overlaps(>=0.5A) into their own list*/ } else if(compArgStr(p+1, "XHvector", 2)){ LXHvector = TRUE; /*20120120dcr*/ } else if(compArgStr(p+1, "ONEDOTeach", 6)){ LOneDotEach = TRUE; /*111013dcr*/ } else if(compArgStr(p+1, "GAPBINs", 6)){ Lgapbins = TRUE; /*111020dcr*/ } else if(compArgStr(p+1, "NEAREST", 7)){ Lnearest = TRUE; /*111021dcr*/ } else if(compArgStr(p+1, "FILTERnearest", 6)){ Lfilternearest = TRUE; /*111022dcr*/ Lnearest = TRUE; /*to make the logic explicit*/ } else if(compArgStr(p+1, "honlyocclude", 5)){ Lhonlyocclude = TRUE; /*130427dcr re ala-dipept gamma turn vs C occlusion*/ } else if(compArgStr(p+1, "OCC1filter", 4)){ Lfilterocc1 = TRUE; /*20120120dcr*/ } else if(compArgStr(p+1, "NOALTfilter", 5)){ Lfilternoalt = TRUE; /*20120120dcr*/ } else if(compArgStr(p+1, "OUTSIDE", 2)){ *method = EXTERNALSURFACE; nargs = 1; } else if(compArgStr(p+1, "ONCE", 2)){ *method = INTERSECTONCE; nargs = 2; } else if(compArgStr(p+1, "BOTH", 1)){ *method = INTERSECTBOTHWAYS; nargs = 2; } else if(compArgStr(p+1, "SELF", 2)){ *method = SELFINTERSECT; nargs = 1; } else if(compArgStr(p+1, "DUMPATOMINFO", 5)){ *method = DUMPATOMCOUNT; nargs = 1; } else if(compArgStr(p+1, "SCSURFACE", 3)){ *method = EXTERNALSURFACE; nargs = 1; argcnt = 1; * srcArg = "not water"; *keepUnselected = FALSE; *probeRad = 1.4f; *groupLabel = "SCS"; } else if(compArgStr(p+1, "EXPOSED", 4)){ *method = EXTERNALSURFACE; nargs = 1; *keepUnselected = FALSE; *probeRad = 1.4f; *groupLabel = "SCS"; } else if(compArgStr(p+1, "ASURFACE", 2)){ *method = EXTERNALSURFACE; nargs = 1; argcnt = 1; * srcArg = "not water"; *keepUnselected = FALSE; RadScaleOffset = 1.4f; *probeRad = 0.0; *groupLabel = "AS"; } else if(compArgStr(p+1, "ACCESS", 6)){ *method = EXTERNALSURFACE; nargs = 1; *keepUnselected = FALSE; RadScaleOffset = 1.4f; *probeRad = 0.0; *groupLabel = "AS"; } else if((n = compArgStr(p+1, "SCAN", 4))){ int scanType = parseInteger(p, n+1, 10); switch(scanType) { case 0: *method = SELFINTERSECT; nargs = 1; argcnt = 1; * srcArg = "alta blt40 ogt33"; Maxbonded = 4; /*changed from 3 : dcr041017*/ DoMcMc = TRUE; break; case 1: *method = INTERSECTONCE; nargs = 2; argcnt = 2; * srcArg = "sc alta blt40 ogt33"; *targArg = "alta blt40 ogt65,(not water ogt33)"; Maxbonded = 4; /*changed from 3 : dcr041017*/ break; default: sprintf(message, "invalid scan type %d, for more info use -help", scanType); halt(message); break; } } else if(compArgStr(p+1, "SEGID", 5)){ SplitBySegID = TRUE; } else if(compArgStr(p+1, "LIMIT", 3)){ LimitDots = TRUE; } else if(compArgStr(p+1, "NOLIMIT", 5)){ LimitDots = FALSE; } else if(compArgStr(p+1, "LENS", 3)){ LensDots = TRUE; } else if(compArgStr(p+1, "NOLENS", 5)){ LensDots = FALSE; } else if(compArgStr(p+1, "MC", 2)){ DoMcMc = TRUE; } else if(compArgStr(p+1, "DUMPH2O", 5)){ DumpNewHO = TRUE; } else if(compArgStr(p+1, "HETS", 3)){ DoHet = TRUE; } else if(compArgStr(p+1, "NOHETS", 5)){ DoHet = FALSE; } else if(compArgStr(p+1, "WAT2WAT", 4)){ DoWatWat = TRUE; /* should be before WATERS */ DoH2O = TRUE; } else if(compArgStr(p+1, "WATERS", 3)){ DoH2O = TRUE; } else if(compArgStr(p+1, "NOWATERS", 5)){ DoH2O = FALSE; } else if(compArgStr(p+1, "4H", 1)){ Maxbonded = 4; } else if(compArgStr(p+1, "3", 1)){ Maxbonded = 3; } else if(compArgStr(p+1, "2", 1)){ Maxbonded = 2; } else if(compArgStr(p+1, "1", 1)){ Maxbonded = 1; } else if(compArgStr(p+1, "POLARH", 1)){ UsePolarH = TRUE; } else if(compArgStr(p+1, "NOPOLARH", 3)){ UsePolarH = FALSE; } else if(compArgStr(p+1, "NUCLEAR", 3)){ NuclearRadii = TRUE; } else if(compArgStr(p+1, "NOFACEHBOND", 6)){ HB2aromFace = FALSE; } else if(compArgStr(p+1, "NAME", 1)){ if (++i < argc) { *groupLabel = argv[i]; } else { halt("no group name after -Name flag"); } } else if(compArgStr(p+1, "DOTMASTER", 9)){ LMasterName = TRUE; /*extra master={name} on lists 060129*/ } else if(compArgStr(p+1, "IGNORE", 6)){ if (++i < argc) { *ignoreArg = argv[i]; } else { halt("no pattern after -IGNORE flag"); } } else if((n = compArgStr(p+1, "COSCALE", 4))){ CORadScale = parseReal(p, n+1, 10); } else if((n = compArgStr(p+1, "HBREGULAR", 3))){ Min_regular_hb_cutoff = parseReal(p, n+1, 10); } else if((n = compArgStr(p+1, "HBCHARGED", 3))){ Min_charged_hb_cutoff = parseReal(p, n+1, 10); } else if(compArgStr(p+1, "DOCHO", 5)){ PermitCHXHB = TRUE; } else if((n = compArgStr(p+1, "CHO", 3))){ CHOHBfactor = parseReal(p, n+1, 10); } else if(compArgStr(p+1, "ATOMCOLOR", 4)){ ColorGap = FALSE; } else if(compArgStr(p+1, "GAPCOLOR", 3)){ ColorGap = TRUE; } else if(compArgStr(p+1, "BASECOLOR", 4)){ ByNABaseColor = TRUE; ColorGap = FALSE; /* forces base color to be only type of display */ } else if(compArgStr(p+1, "COLORBASE", 6)){ ByNABaseColor = TRUE; } else if((n = compArgStr(p+1, "GAPWEIGHT", 4))){ GAPweight = parseReal(p, n+1, 10); } else if((n = compArgStr(p+1, "BUMPWEIGHT", 5))){ BUMPweight = parseReal(p, n+1, 10); } else if((n = compArgStr(p+1, "HBWEIGHT", 3))){ HBweight = parseReal(p, n+1, 10); } else if((n = compArgStr(p+1, "DIVLow", 4))){ LowGoodCut = parseReal(p, n+1, 10); } else if((n = compArgStr(p+1, "DIVHigh", 4))){ HighGoodCut = parseReal(p, n+1, 10); } else if((n = compArgStr(p+1, "DIVWorse", 4))){// 04/21/2015 SJ to change the default of the worse overlap WorseBumpCut = parseReal(p, n+1, 10); } else if((n = compArgStr(p+1, "MINOCCupancy", 6))){ OccupancyCutoff = parseReal(p, n+1, 10); } else if((n = compArgStr(p+1, "NOOCCupancy", 5))){ OccupancyCriteria = FALSE; /*060831*/ } else if(compArgStr(p+1, "OFORMAT", 7)){ OutputFmtType = 1; } else if(compArgStr(p+1, "XVFORMAT", 8)){ OutputFmtType = 2; } else if(compArgStr(p+1, "ONELINE", 7)){ /*oneline dcr041101*/ ContactSUMMARY = TRUE; OutputFmtType = 3; argcnt = 1; /*recognize next naked item as an input file name*/ } else if(compArgStr(p+1, "SUMMARY", 7)){ /*summary dcr041101*/ ContactSUMMARY = TRUE; argcnt = 1; /*recognize next naked item as an input file name*/ } else if(compArgStr(p+1, "DEFAULTS", 7)){ /*defaults dcr041101*/ *srcArg = "altA ogt33"; *method = SELFINTERSECT; Maxbonded = 4; /*changed from 3 : dcr041017*/ DoMcMc = TRUE; DoHet = TRUE; argcnt = 1; /*recognize next naked item as an input file name*/ /*like naked defaults but allow flags like -summary or -oneline*/ } else if(compArgStr(p+1, "NOHBOUT", 7)){ OutputHBs = FALSE; } else if(compArgStr(p+1, "NOCLASHOUT", 10)){ OutputClashes = FALSE; } else if(compArgStr(p+1, "NOVDWOUT", 8)){ OutputVDWs = FALSE; } else if(compArgStr(p+1, "ONLYBADOUT", 10)){ /*dcr041010*/ OnlyBadOut = TRUE; } else if(compArgStr(p+1, "NOTICKS", 6)){ ShowTicks = FALSE; } else if(compArgStr(p+1, "STDBONDS", 7)){ UseStdBond = TRUE; } else if(compArgStr(p+1, "OLDU", 4)){ OldStyleU = TRUE; *rawOutput = TRUE; } #ifdef JACKnANDREA else if(compArgStr(p+1, "WRITEHBONDS", 6)){ writeHbonds = TRUE; } #endif else if(compArgStr(p+1, "NOPARENT", 8)){ UseHParent = FALSE; } else if(compArgStr(p+1, "OUTCOLOR", 6)){ if (++i < argc) { OutPointColor = argv[i]; } else { halt("no color name after -OUTCOLolor flag"); } } else if(compArgStr(p+1, "AUTObondrot", 4)){ if (++i < argc) { p = argv[i]; if (p[0] == '-') { if (p[1] == '\0') { file++; if (mabip) { mabip->filenum = file; mabip->inf = stdin; /*input of autobondrot atoms*/ mabip->close = FALSE; } } else { sprintf(message, "%s flag instead of filename after -AUTObondrot", p); halt(message); } } else { file++; inf = fopen(p, "rb"); /*p holds autobondrot input file name*/ if (inf) { if (mabip) { mabip->filenum = file; mabip->inf = inf; mabip->close = TRUE; } } else { sprintf(message, "could not open -AUTOBONDROT file: %s", p); halt(message); } } /*Lautobondrot = TRUE;*/ /*050111 probe.c global logical */ *countDots = TRUE; *rawOutput = TRUE; /*autobondrot mode does NOT read from input file at this time*/ } else { halt("no filename after -AUTObondrot flag"); } } else if(compArgStr(p+1, "REFerence", 3)){ sprintf(message, "Please cite: %s", referenceString); note(message); sprintf(message, "For more information see %s", electronicReference); note(message); exit(0); } else if((n = compArgStr(p+1, "CHANGEs", 6))){ dump_changes(stderr); } else if((n = compArgStr(p+1, "DEBUG", 5))){ DebugLevel = parseInteger(p, n+1, 10); } else { /*either naked - or unrecognized character string*/ sprintf(message, "unrecognized flag, %s", p); halt(message); } }/* -flag item */ else if (argcnt <= 0) { /*lonely naked item: store as a source argument*/ *srcArg = p; /*if nothing else comes in, use as input file name*/ argcnt = 1; } else if (argcnt == 1 && nargs == 2) { if((p[0] == '=') && (p[1] == '\0')) { *targArg = *srcArg; } else { *targArg = p; } argcnt = 2; } else if(Lfilternearest && argcnt == 2) { *extraArg = p; /*111021dcr*/ /* -onedotonly -nearest "srcPat" "targPat" "extraPat" output filter*/ argcnt = 3; } else { file++; inf = fopen(p, "rb"); if (inf) { strcpy(inputfilename,p); /*dcr041023*/ /*read atoms from input file now*/ atomlist = loadAtoms(inf, atomlist, bboxA, file, reslstptr); fclose(inf); } else { sprintf(message, "could not open file %d: %s", file, p); halt(message); } } } if (file < 1) { if (argc == 2 && argcnt == 1 && *srcArg) { /*naked item stored as srcArg*/ p = *srcArg; /*interpret as an input file name, then presume defaults*/ note("*SIMPLE SELF DOTS (for more info use -help flag and look for <>)"); *srcArg = "altA ogt33"; *method = SELFINTERSECT; Maxbonded = 4; /*changed from 3 : dcr041017*/ DoMcMc = TRUE; DoHet = TRUE; file++; inf = fopen(p, "rb"); if (inf) { strcpy(inputfilename,p); /*dcr041023*/ /*read atoms from input file now*/ atomlist = loadAtoms(inf, atomlist, bboxA, file, reslstptr); fclose(inf); } else { sprintf(message, "could not open file %d: %s", file, p); halt(message); } } else { errmsg("no input files"); probehelp(0); } } return atomlist; /*loaded under some conditions, but not if autobondrot!*/ }/*processCommandline()*/ /*}}}processCommandline()____________________________________________________*/ /*{{{initEndData()************************************************************/ void initEndData(chainEndData_t *ed) { int i = 0; for (i = 0; i < 8; i++) { ed->ambigN[i] = NULL; ed->ambigO[i] = NULL; } ed->res_mc_oxy_cnt = 0; ed->first = 0; ed->last = 0; ed->Ntmarkers = 0; ed->Ctmarkers = 0; for (i = 0; i < 4; i++) { ed->hisNHcount[i] = 0; } } /*}}}initEndData()___________________________________________________________*/ /*{{{hisNHcheck()*************************************************************/ /* His ring atoms start out with a positive charge. */ /* Here we remove the charge for His residues without two ring NH protons.*/ void hisNHcheck(chainEndData_t *ed, atom *atomlist, int rescnt) { int i = 0; atom *rat = NULL; if (ed->hisNHcount[0]) { /* the zeroth counter flags: Is this a HIS? */ /*is a histidine*/ for(rat = atomlist; rat; rat = rat->next) { if ((rat->r->rescnt == rescnt) && (rat->props & POSITIVE_PROP)) { char altc = toupper(rat->altConf); int ialt = 0; if ((altc==' ')||(altc=='A')||(altc=='1')) { ialt = 1; } else if ((altc=='B')||(altc=='2')) { ialt = 2; } else if ((altc=='C')||(altc=='3')) { ialt = 3; } if (ialt > 0) { if (ed->hisNHcount[ialt] < 2) { /* his not positive */ rat->props &= ~POSITIVE_PROP; } else { /* his positive but not an acceptor */ rat->props &= ~ACCEPTOR_PROP; } } } } for (i = 0; i < 4; i++) { ed->hisNHcount[i] = 0; } }/*is a histidine*/ } /*}}}hisNHcheck()____________________________________________________________*/ /*{{{resCheck()***************************************************************/ /* called after each residue is read */ void resCheck(chainEndData_t *ed, atom *atomlist, int rescnt) { hisNHcheck(ed, atomlist, rescnt); } /*}}}resCheck()______________________________________________________________*/ /*{{{CtermCheck()*************************************************************/ /* called when a new residue begins and at the end of all processing */ /* to process the previous residue */ void CtermCheck(chainEndData_t *ed, int rescnt, int isChainEnd) { int i = 0; /* avoid processing the first time 'cause there ain't no chain pending */ if (isChainEnd && rescnt) { ed->last = rescnt; /* last residue */ /* see if we can put the pieces together to determine end charge */ for (i = 0; i < 4; i++) { /* only array[0-3] contains first Ns */ if (ed->ambigN[i]) { if (ed->ambigN[i]->r->rescnt == ed->first && ed->Ntmarkers == ed->first) { ed->ambigN[i]->props |= POSITIVE_PROP; } } else { break; } } for (i = 0; i < 8; i++) { /* array[0-7] contains last Os */ if (ed->ambigO[i]) { if (ed->ambigO[i]->r->rescnt == ed->last && ed->Ctmarkers == ed->last) { ed->ambigO[i]->props |= NEGATIVE_PROP; } } else { break; } } initEndData(ed); /* reset the end data record */ } } /*}}}CtermCheck()____________________________________________________________*/ /*{{{noticedNt()**************************************************************/ /* called when we have an indicator that residue is at N terminus */ /* to look for ambigN and mark charged */ void noticedNt(chainEndData_t *ed, int rescnt) { int i = 0; for (i = 4; i < 8; i++) { /* only array[4-7] contains last Ns */ if (ed->ambigN[i]) { if (ed->ambigN[i]->r->rescnt == rescnt) { ed->ambigN[i]->props |= POSITIVE_PROP; } } else { break; } } } /*}}}noticedNt()_____________________________________________________________*/ /*{{{noticedCt()**************************************************************/ /* called when we have an indicator that residue is at C terminus */ /* to look for ambigO and mark charged */ void noticedCt(chainEndData_t *ed, int rescnt) { int i = 0; for (i = 0; i < 8; i++) { /* array[0-7] contains last Os */ if (ed->ambigO[i]) { if (ed->ambigO[i]->r->rescnt == rescnt) { ed->ambigO[i]->props |= NEGATIVE_PROP; } } else { break; } } } /*}}}noticedCt()_____________________________________________________________*/ /*{{{NtermCheck()*************************************************************/ /* called whenever a new residue begins */ void NtermCheck(chainEndData_t *ed, int rescnt, int isChainEnd) { int i = 0; ed->res_mc_oxy_cnt = 0; /* reset Oxygen data for each residue */ for (i = 0; i < 8; i++) { ed->ambigO[i] = NULL; } /* we have reserved the top 4 nitrogen slots for most recent res */ /* so we have to clear this here */ for (i = 4; i < 8; i++) { ed->ambigN[i] = NULL; } if (isChainEnd) { ed->first = rescnt; /* first residue */ } } /*}}}NtermCheck()____________________________________________________________*/ /*{{{ProcessResInfo()*** called from loadAtoms() & movingAtomListProcessing()*/ /* called for each atom read in */ void ProcessResInfo(chainEndData_t *ed, atom *a) { int i = 0; if (strstr("HIS:his", a->r->resname)) {/*is a histidine*/ ed->hisNHcount[0] = 1; /* indicates this is a HIS */ if ( strstr(" HD1: hd1: HE2: he2", a->atomname)) { char altc = toupper(a->altConf); if ((altc==' ')||(altc=='A')||(altc=='1')) { ed->hisNHcount[1]++; } else if ((altc=='B')||(altc=='2')) { ed->hisNHcount[2]++; } else if ((altc=='C')||(altc=='3')) { ed->hisNHcount[3]++; } } }/*is a histidine*/ /* look for indicators that we are at the end of a chain */ if (a->props & CHECK_ENDS_PROP) { if (isHatom(a->elem)) { if (! ed->Ntmarkers) { ed->Ntmarkers = a->r->rescnt; } noticedNt(ed, a->r->rescnt); } else if ( (a->elem == atomO) && ((a->altConf==' ')||(a->altConf=='A')||(a->altConf=='1')) ){ ed->res_mc_oxy_cnt++; if (ed->res_mc_oxy_cnt == 2) { ed->Ctmarkers = a->r->rescnt; noticedCt(ed, a->r->rescnt); } } else if (strstr(a->atomname, "OXT")) { /* check SUBSET of name */ ed->Ctmarkers = a->r->rescnt; noticedCt(ed, a->r->rescnt); } } /* save pointers to the first ambigNs and the last ambigOs and Ns */ /* we save multiple atoms to both find both possible mc oxygens */ /* and to handle multiple conformations (we can do 4 max) */ /* All atoms must have the same residue counter */ if (a->props & MAYBECHG_PROP) { if (!strcmp(a->atomname, " N ")) { if (ed->first == a->r->rescnt) { /* mark the Nterm -- jmw 20011001 */ a->props |= CHECK_ENDS_PROP; } for (i = 0; i < 4; i++) {/* first Ns[0-3] (cleared at chain ends)*/ if (ed->ambigN[i] == NULL) { if ((i == 0) || (ed->ambigN[i-1]->r == a->r)) { ed->ambigN[i] = a; } break; } } for (i = 4; i < 8; i++) { /* last Ns[4-7] (cleared each res)*/ if (ed->ambigN[i] == NULL) { if ((i == 4) || (ed->ambigN[i-1]->r == a->r)) { ed->ambigN[i] = a; } break; } } } if (!strcmp(a->atomname, " O ")) { /* last Os[0-7] (cleared each res)*/ for (i = 0; i < 8; i++) { if (ed->ambigO[i] == NULL) { if ((i == 0)|| (ed->ambigO[i-1]->r == a->r)) { ed->ambigO[i] = a; } break; } } } } } /*}}}ProcessResInfo()________________________________________________________*/ /*{{{movingAtomListProcessing()***********************************************/ void movingAtomListProcessing(atom *initAtomLst, void *userdata) { int previd = 0, rescnt = 0; char prevInsCode = ' '; char prevChain[5]; atom *a = NULL, *prevAtomLst = NULL, *nexta = NULL; chainEndData_t endData; initEndData(&endData); prevChain[0] = prevChain[1] = '?'; prevChain[2] = '\0'; for(a=initAtomLst; a; a = nexta) { nexta = a->next; a->next = prevAtomLst; /* remove a from initial atom list */ prevAtomLst = a; /* and build up a reversed list of */ /* only atoms previously seen */ if (a->r->resid != previd || a->r->resInsCode != prevInsCode || strcmp(a->r->chain, prevChain) != 0) { if (rescnt){ resCheck(&endData, prevAtomLst, rescnt); } CtermCheck(&endData, rescnt, (strcmp(a->r->chain, prevChain) != 0)); ++rescnt; NtermCheck(&endData, rescnt, (strcmp(a->r->chain, prevChain) != 0)); previd = a->r->resid; prevInsCode = a->r->resInsCode; strcpy(prevChain, a->r->chain); //prevChain[2] = '\0'; } a->r->rescnt = rescnt; ProcessResInfo(&endData, a); } if (rescnt) { resCheck(&endData, prevAtomLst, rescnt); } CtermCheck(&endData, rescnt, TRUE); } /*}}}movingAtomListProcessing()______________________________________________*/ /*{{{loadAtoms()************* called from processCommandline() ***************/ atom* loadAtoms(FILE *fp, atom *atomlist, region *boundingBox, int file, residue **reslstptr) { int previd = 0, rescnt = 0, model = 0; char *rec, prevInsCode = ' '; char prevChain[5]; atom *a = NULL; residue *scratchRes; chainEndData_t endData; prevChain[0] = prevChain[1] = '?'; prevChain[2] = '\0'; scratchRes = newResidueData(); initEndData(&endData); while((rec = getPDBrecord(fp))) { /*loop over all records in pdb file*/ if (isTer(rec)) { prevChain[0] = prevChain[1] = '?'; prevChain[2] = '\0'; } if ((isAtom(rec) || isHet(rec)) && ! isPseudoAtom(rec)) {/*atom or htatm*/ a = newAtom(rec, file, model, scratchRes); if (!a) { break;} if (a->r->resid != previd || a->r->resInsCode != prevInsCode || strcmp(a->r->chain, prevChain) != 0) { if (rescnt){ resCheck(&endData, atomlist, rescnt); } CtermCheck(&endData, rescnt, (strcmp(a->r->chain, prevChain) != 0)); ++rescnt; NtermCheck(&endData, rescnt, (strcmp(a->r->chain, prevChain) != 0)); previd = a->r->resid; prevInsCode = a->r->resInsCode; strcpy(prevChain, a->r->chain); //prevChain[2] = '\0'; } a->r->rescnt = rescnt; ProcessResInfo(&endData, a); if (ImplicitH && isHatom(a->elem)) { deleteAtom(a); /* filter out implicit hydrogens */ continue; } if (isDummyAtom(*a)) { deleteAtom(a); /* filter out dummy Xplor atoms */ continue; } if (!atomlist) { boundingBox->min = a->loc; boundingBox->max = a->loc; } a->next = atomlist; atomlist = a; updateBoundingBox(&(a->loc), boundingBox); if (*reslstptr == NULL) {/* first residue goes on residue list */ *reslstptr = scratchRes; scratchRes = newResidueData(); } else { /* otherwise we have to compare residue data blocks */ if (resDiffersFromPrev(*reslstptr, scratchRes)) { scratchRes->nextRes = *reslstptr; *reslstptr = scratchRes; scratchRes = newResidueData(); } else { a->r = *reslstptr; /* same, so point to prior block */ a->r->a = a; /* makes this atom the head of atom list for this res */ } } }/*atom or htatm*/ if (isModel(rec)) { model = parseModel(rec); modelNumber[++modelCount] = model; /*041114*/ if(Verbose) { fprintf(stderr,"modelNumber[%d]==%d\n",modelCount,modelNumber[modelCount]); } } }/*loop over all records in pdb file*/ if (rescnt) { resCheck(&endData, atomlist, rescnt); } CtermCheck(&endData, rescnt, TRUE); deleteResidueData(scratchRes); /* clean up any extra residue */ return atomlist; } /*}}}loadAtoms()_____________________________________________________________*/ /*{{{binAtoms()***************************************************************/ atomBins* binAtoms(atom *theAtoms, region *boundingBox, char serialNum, float probeRad, int keepUnselected, int selflags) { atom *a = NULL; atomBins *bins = NULL; bins = initBins(serialNum, boundingBox, 2.0f*(getMaxRadius(ImplicitH)+probeRad)+0.2f); if (bins) { for(a = theAtoms; a; a = a->next) { if (keepUnselected || (a->flags & selflags)) { if (! (a->flags & IGNORE_FLAG)) { addNeighbor(a, bins); } /*050118 this seems the only place that uses IGNORE_FLAG*/ } } } /* note: addNeighbor() also called in updateHydrogenInfo() */ return bins; } /*}}}binAtoms()______________________________________________________________*/ /*{{{getRadius()**************************************************************/ float getRadius(int at, int useCOScale) { float rad = 0.0, sf = RadScaleFactor; if (useCOScale){ sf *= CORadScale; } if (ImplicitH) { rad = getImplRad(at); } else { rad = getExplRad(at); } return (rad*sf) + RadScaleOffset; } /*}}}getRadius()_____________________________________________________________*/ /*{{{newRingInfo()************************************************************/ ringInfo * newRingInfo(ringInfo **head, point3d* ctr, point3d* norm) { ringInfo *ri = NULL; ri = (ringInfo *)malloc(sizeof(ringInfo)); if (ri) { ri->nextRing = NULL; ri->ringCenter = *ctr; ri->ringNormal = *norm; if (head) { /* if list passed, link new ring as new head */ ri->nextRing = *head; *head = ri; } } return ri; } /*}}}newRingInfo()___________________________________________________________*/ /*{{{deleteRingInfoList()*****************************************************/ void deleteRingInfoList(ringInfo *ri) { /* kill entire list */ ringInfo *p = NULL, *nxt = NULL; p = ri; while(p) { nxt = p->nextRing; free(p); p = nxt; } } /*}}}deleteRingInfoList()____________________________________________________*/ /*{{{newResidueData()*********************************************************/ residue * newResidueData() { /* new blank residue */ residue *r = NULL; r = (residue *)malloc(sizeof(residue)); if (r) { r->nextRes = NULL; r->a = NULL; r->file = 0; r->model = 0; //r->chain[0] = ' '; /*RMI I don't know what to do here*/ r->chain[0] = ' '; r->chain[1] = ' '; r->chain[2] = '\0'; r->resid = ' '; r->resInsCode = ' '; r->rescnt = 0; r->segid[0] = '\0'; r->ring = NULL; } return r; } /*}}}newResidueData()________________________________________________________*/ /*{{{resDiffersFromPrev()*****************************************************/ int resDiffersFromPrev(residue *r1, residue *r2) { if (r1 == NULL || r2 == NULL) { return 1; } return (! IS_THE_SAME_RES(r1, r2)); } /*}}}resDiffersFromPrev()____________________________________________________*/ /*{{{deleteResidueData()******************************************************/ void deleteResidueData(residue *r) { if (r) { deleteRingInfoList(r->ring); free(r); } } /*}}}deleteResidueData()_____________________________________________________*/ /*{{{disposeListOfResidues()**************************************************/ void disposeListOfResidues(residue *theRes) { residue *r = NULL, *nextr = NULL; for(r=theRes; r; r = nextr) { nextr = r->nextRes; deleteResidueData(r); } } /*}}}disposeListOfResidues()_________________________________________________*/ /*{{{dumpRes()****************************************************************/ /* dumpRes() used for debugging */ void dumpRes(residue *theRes) { residue *r = NULL; atom *a = NULL; for(r = theRes; r; r = r->nextRes) { fprintf(stderr, "RES{%s%d%c}[%d] ring(%p)\n", r->resname, r->resid, r->resInsCode, r->rescnt, r->ring); fprintf(stderr, " "); for(a = r->a; a && (a->r == r); a = a->next) { fprintf(stderr, "%s ", a->atomname); } fprintf(stderr, "\n"); } } /*}}}dumpRes()_______________________________________________________________*/ /*{{{deleteAtom()*************************************************************/ void deleteAtom(atom *a) { if (a) { free(a); } } /*}}}deleteAtom()____________________________________________________________*/ /*{{{disposeListOfAtoms()*****************************************************/ void disposeListOfAtoms(atom *theAtoms) { atom *a = NULL, *nexta = NULL; for(a=theAtoms; a; a = nexta) { nexta = a->next; deleteAtom(a); } } /*}}}disposeListOfAtoms()____________________________________________________*/ /*{{{newAtom() <--loadAtoms(),newMovingAtom(); -->various property tests *****/ atom * newAtom(char *rec, int file, int model, residue * resDataBuf) {/*newAtom()*/ atom *a = NULL; char msg[100]; a = (atom *)malloc(sizeof(atom)); if (a) { a->r = resDataBuf; a->next = NULL; a->nextInBin= NULL; a->scratch = NULL; a->mark = 0; /* used to mark bonds */ a->flags = 0; a->props = 0; parseResidueName(rec, a->r->resname); parseAtomName(rec, a->atomname); a->r->a = a; /* residue points back to this atom (can change for res) */ a->r->file = file; a->r->model = model; parseChain(rec, a->r->chain); /* a->r->chain = parseChain(rec);*/ a->r->resid = parseResidueNumber(rec); parseResidueHy36Num(rec, a->r->Hy36resno); a->r->resInsCode = parseResidueInsertionCode(rec); a->r->rescnt = 0; a->altConf = parseAltLocCode(rec); a->binSerialNum = '?'; /* set when we bin */ if (strstr(":ASX:GLX:ASN:GLN:", a->r->resname) /* special case treats undecided */ && (a->atomname[1] == 'A')) { /* as an oxygen */ a->elem = identifyAtom(" O ", a->r->resname, Verbose); /*dcr041007 allow warning add resname to call */ sprintf(msg, "atom %s will be treated as oxygen", a->atomname); warn(msg); } else { /* normal case */ a->elem = identifyAtom(a->atomname, a->r->resname, Verbose);/*dcr041007 allow warning add resname to call */ } /*next section seems to be the only place where atom->bondedto is set.*/ /*select.c/setHydrogenParentName() and select.c/setMainchainBonding()*/ /* both merely call stdconntable.c/searchForStdBondingPartner() */ /* which constructs a search string from (a->r->resname, a->atomname) */ /* and calls stdconntable.c/SearchStdResConnTable() */ /* which returns the string found in the StdResTblBucket[] hash table.*/ /*This hash table is a marvelous tour de force listing of names of*/ /*residue/atoms and what atoms they can be bonded to */ a->bondedto = NULL; if (isHatom(a->elem)) { a->atomclass = -1; /* to be specified later */ if (UseHParent) { a->bondedto = setHydrogenParentName(a->r->resname, a->atomname); } } else { a->atomclass = a->elem; if (UseStdBond) { a->bondedto = setMainchainBonding(a->r->resname, a->atomname); /*setMainchainBonding() does not seem to have a distance limit*/ } } setProperties(a, isHet(rec), HB2aromFace, PermitCHXHB); /*select.c*/ /*This is where (e.g.) both HET_PROP and DNA_PROP are set for the atom*/ if (ByNABaseColor) { /* forces coloring by base rather than atom type */ a->atomclass = naBaseCategory(a); } parseXYZ(rec, &(a->loc)); a->ix = a->iy = a->iz = 0; /* note: H DONOR_PROP assignment done later in updateHydrogenInfo() */ a->occ = parseOccupancy(rec); a->bval = parseTempFactor(rec); parseSegID(rec, a->r->segid); //use SEGID instead of chainID if it exists and chainID is blank if ( (strncmp(a->r->segid, " ",4) != 0) && (strncmp(a->r->chain, " ",2) == 0)){ //parseSegID(rec, a->r->chain); a->r->chain[0] = a->r->segid[0]; a->r->chain[1] = a->r->segid[1]; a->r->chain[2] = a->r->segid[2]; a->r->chain[3] = a->r->segid[3]; a->r->chain[4] = '\0'; } a->radius = getRadius(a->elem, isCarbonylAtom(a)); a->covRad = getCovRad(a->elem); } else { warn("could not allocate space for new atom"); } return a; }/*newAtom()*/ /*}}}newAtom()_______________________________________________________________*/ /*{{{selectSource() <--mainProbeProc(),newMovingAtom(); -->select/matchPat() */ void selectSource(atom *theAtoms, pattern *srcPat, int srcFlag, pattern *targPat, int targFlg, pattern *ignorePat) { atom *src = NULL; for(src = theAtoms; src; src = src->next) { if( (DoHet || ! (src->props & HET_PROP) ) && (DoH2O || ! (src->props & WATER_PROP) ) ) { /*060212 hets marked as prot,dna,rna could be excluded by these tests*/ /* if src||target was specified as not-prot,dna,rna,... */ /* seemingly because atoms of hets are also assigned those types */ if (srcPat && matchPat(src, srcPat)) { src->flags |= srcFlag; } if(targPat && matchPat(src, targPat)) { src->flags |= targFlg; } /*if ((FABS(src->occ) <= OccupancyCutoff) */ /* || (ignorePat && matchPat(src, ignorePat))) */ /*050119 not use occ cutoff here, so 0-occ atom will be put into bins*/ if (ignorePat && matchPat(src, ignorePat)) { src->flags = IGNORE_FLAG; /* overwrite any other settings */ /*050118 this seems the only place that sets IGNORE_FLAG*/ } } } } /*}}}selectSource()__________________________________________________________*/ /*{{{atomsClose()*************************************************************/ /*atomsClose() only called from findTouchingAtoms() */ int atomsClose(atom *a, atom *b, float probeRad) { int nearpt = FALSE; float lim = 0.0, dsq = 0.0; lim = a->radius + b->radius + probeRad + probeRad; dsq = (float)v3distanceSq(&(a->loc), &(b->loc)); /* if too close they must be the same atom...dummy atom?? */ /*removing the (dsq > 0.001) test actually removes one side of a clash!!*/ /*20111207 dcr changed test to .003*.003 = .00001 i.e. .001 in each coord */ /* to minimize false same-atom-ness but still allow filtering of dummy atoms*/ if ((dsq > 0.00001) && (dsq <= (lim*lim))) { nearpt = TRUE; } return nearpt; } /*}}}atomsClose()____________________________________________________________*/ /*{{{inRange()****************************************************************/ int inRange(point3d *p, point3d *q, float lim) { return v3distanceSq(p, q) <= (lim*lim); } /*}}}inRange()_______________________________________________________________*/ /*{{{perpendicular()**********************************************************/ void perpendicular(point3d *p, point3d *q, point3d *r, point3d *s) { /*cribbed from prekin perpendicular, 111018 */ /*construct perpendicular from point to line*/ /*perpendicular from 3rd point to line of points 1 and 2 */ /* fourth point is on line at perpendicular to third point*/ float A1X,A1Y,A1Z,A2X,A2Y,A2Z,B1X,B1Y,B1Z,ADX,ADY,ADZ,B2X,B2Y,B2Z,ART; A1X = p->x; A1Y = p->y; A1Z = p->z; A2X = q->x; A2Y = q->y; A2Z = q->z; B1X = r->x; B1Y = r->y; B1Z = r->z; ADX = A2X - A1X;ADY = A2Y - A1Y;ADZ = A2Z - A1Z; /*arbitrary point on line A B2X = A1X + ADX*ART;B2Y = A1Y + ADY*ART;B2Z = A1Z + ADZ*ART; direction components of line B1 --- B2 BDX = A1X-B1X+ADX*ART; BDY = A1Y-B1Y+ADY*ART; BDZ = A1Z-B1Z+ADZ*ART; for this to be perpendicular to line A, sum of product of components ==0 ADX*(A1X-B1X+ADX*ART)+ADY*(A1Y-B1Y+ADY*ART)+ADZ*(A1Z-B1Z+ADZ*ART)=0; rearrange ADX*A1X-ADX*B1X+ADX*ADX*ART +ADY*A1Y-ADY*B1Y+ADY*ADY*ART +ADZ*A1Z-ADZ*B1Z+ADZ*ADZ*ART=0; regroup ART*(ADX*ADX++ADY*ADY+ADZ*ADZ) =ADX*B1X-ADX*A1X+ADY*B1Y-ADY*A1Y+ADZ*B1Z-ADZ*A1Z; */ ART = (ADX*B1X-ADX*A1X+ADY*B1Y-ADY*A1Y+ADZ*B1Z-ADZ*A1Z) /(ADX*ADX+ADY*ADY+ADZ*ADZ); B2X = A1X + ADX*ART; B2Y = A1Y + ADY*ART; B2Z = A1Z + ADZ*ART; s->x = B2X; s->y = B2Y; s->z = B2Z; } /*}}}perpendicular()_________________________________________________________*/ /*{{{gapSize()****************************************************************/ float gapSize(point3d *p, point3d *q, float qrad) { return v3distance(p, q) - qrad; } /*}}}gapSize()_______________________________________________________________*/ /*{{{findTouchingAtoms()******************************************************/ /* findTouchingAtoms() - Note: resets bond marks! Returns list of close atoms */ /*that are in the set of bins here called abins */ /*in the form of the first atom that is close to the src-atom where */ /*a linked list of all close atoms is made by the scratch member of each atom*/ /*which indicates another neighbor of the src-atom */ atom* findTouchingAtoms(atom *src, atom *head, atomBins *abins, float probeRad, int targFlg, int *ok) {/*findTouchingAtoms()*/ atom *a = NULL; int i = 0, j = 0, k = 0, nearpt = 0, targcount = 0; int jx = (src->ix), jy = (src->iy), jz = (src->iz); int imin = 0, jmin = 0, kmin = 0; int imax = 0, jmax = 0, kmax = 0; if (src->binSerialNum != abins->binSerialNum) {/*in foreign bins */ /*dcr?: why do they have to be in foreign==different bins??*/ /*this requirement seems not to affect loss of close clashes*/ /* must look up the positions */ if( (src->loc.x < (abins->min.x - abins->delta)) || (src->loc.y < (abins->min.y - abins->delta)) || (src->loc.z < (abins->min.z - abins->delta)) || (src->loc.x > (abins->max.x + abins->delta)) || (src->loc.y > (abins->max.y + abins->delta)) || (src->loc.z > (abins->max.z + abins->delta)) ) { return NULL; /* nothing anywhere nearby */ } binLoc(&(src->loc), &jx, &jy, &jz, abins); #ifdef DEBUG_A2B fprintf(stderr, "f(%c != %c) %s%d %s [%d, %d, %d] <1..%d, 1..%d, 1..%d>\n", src->binSerialNum, abins->binSerialNum, src->r->resname, src->r->resid, src->atomname, jx, jy, jz, abins->nx-2, abins->ny-2, abins->nz-2); #endif }/*in foreign bins */ imin = jx-1; imax = jx+1; /* bin ranges */ jmin = jy-1; jmax = jy+1; kmin = jz-1; kmax = jz+1; /* trim any excess edges */ if (imin < 1) { imin = 1; } if (jmin < 1) { jmin = 1; } if (kmin < 1) { kmin = 1; } if (imax > (abins->nx - 2)) { imax = (abins->nx - 2); } if (jmax > (abins->ny - 2)) { jmax = (abins->ny - 2); } if (kmax > (abins->nz - 2)) { kmax = (abins->nz - 2); } for(i = imin; i <= imax; i++) { for(j = jmin; j <= jmax; j++) { for(k = kmin; k <= kmax; k++) { for(a = abins->list[i][j][k]; a; a = a->nextInBin) { /*Lmodeltest for interaction of diff. models 041112*/ /*050121 Lmodeltest superseded, mage sends model# */ /*...within subroutine findTouchingAtoms()... */ if( (src->r->model != a->r->model) ) { continue; } if( modelLimit > 1 ) /*041114*/ {/*looping through multiple models*/ if(src->r->model != modelToProcess) { continue; } } /*if (( src->r == a->r) && -- alternate */ /* alternate */ if (( src->altConf != a->altConf) && /* conformations */ ( src->altConf != ' ') && /* don't interact */ ( a->altConf != ' ') ) { continue; } /*nearpt seems to be essential to get any dots or spikes */ /*...seems not to be where too-close atoms fail to get any dots or spikes*/ nearpt = atomsClose(src, a, probeRad);/*only call to atomsClose*/ if (nearpt) {/*add this atom to the linked-list of those touching src atom*/ a->mark = 0; /* clear the bonding marker. */ a->scratch = head; /* link using scratch pointers */ head = a; /* with at head of list */ if (a->flags & targFlg) { targcount++; } } } } } } src->mark = 0; /* clear the src atom bonding marker. */ *ok = !(targFlg && (targcount < 1)); return head; }/*findTouchingAtoms*/ /*}}}findTouchingAtoms()_____________________________________________________*/ /*{{{dotClassIndex()**********************************************************/ /* dotClassIndex() - maps dot type to class index */ /* OLD t == 0 --> 0,1 : contact dot */ /* t < 0 --> 2,3 : bump */ /* t > 0 --> 4 : hbond */ /* NEW - SJ t == 0 --> 0,1 : contact dot */ /* 04/10/2015 t < 0 --> 3,4,5 : small, bad, worse bump*/ /* t > 0 --> 2,6 : weak Hbonds, Hbonds*/ int dotClassIndex(int t, float mingap) { int idx = 0; // Old code for old definition //if (t == 0) { idx = (mingap > HighGoodCut) ? 0 : 1; } /* contact */ //else if (t < 0) { idx = (mingap > LowGoodCut) ? 2 : 3; } /* clash */ //else if (LweakHbonds) //{ // if (t==1) {idx = 4; } /*good hbond*/ // else {idx = 5; } /*weak hbond 20111215dcr*/ //} //else { idx = 4; } /* hbonds */ /* new code SJ 04/10/2105 to implement the new definiion above*/ if (t == 0) { /* contact */ if (mingap > HighGoodCut) idx = 0; /* wide contact*/ else idx = 1; /* close contact*/ } else if(t < 0){ /* clash*/ if (mingap > LowGoodCut) idx = 3; /* small overlap*/ else if(LworseOverlap) /* divide bad and worse overlap only if this is true*/ { if(mingap > WorseBumpCut) idx = 4; /* bad overlap */ else idx = 5; /* worse overlap */ } else idx = 4; /* bad overlap*/ } else{ /* H bond */ if(LweakHbonds){/* divide H bonds into weak only if this is true.*/ if(mingap > 0) idx = 2; /* weak H bond*/ else idx = 6; /* H bonds*/ } else idx = 6; /* H bonds*/ } return idx; } /*}}}dotClassIndex()_________________________________________________________*/ /*{{{saveDot() called from examineDots() and from SurfDots() *****************/ void saveDot(atom *src, atom *targ, int type, point3d *loc, point3d *spike, dotNode *results[][NODEWIDTH], int ovrlaptype, float mingap, char ptmaster, float angle) {/*saveDot()*/ /*ptmaster dcr041009*/ dotNode* dot = NULL; int which = 0, idx = 0; /*overlaptype: -1 bump, 0 touch, +1 H bond */ /*entering from SurfDots needs no further Logical filtering*/ /*move this decision making back up into examineDots()*/ /*where we now will accummulate count numbers also*/ #ifdef OLDCODE /*041020 saveDot() decisions moved up into examineDots*/ if ((OutputHBs && ovrlaptype > 0) || (OutputClashes && ovrlaptype < 0) || (OutputVDWs && ovrlaptype == 0)) { #endif idx = dotClassIndex(ovrlaptype, mingap); /*OLD: 0 wide contact,1 close contact,2 small overlap,3 bad overlap,4 H-bonds*/ /*OLD: LweakHbonds: 4 good H-bonds, 5 weak (distant) H-bonds 20111215dcr*/ /* NEW 04/10/2015 SJ: 0 wide contact, 1 close contact, 2 weak H-bonds, 3 small overlap, 4 bad overlap, 5 worse overlap, 6 H-bonds*/ /*idx=1 for EXTERNALSURFACE ovrlaptype==0, mingap==0.0 */ dot = newDot(src, targ, loc, spike, ovrlaptype, mingap, ptmaster, angle); /*ptmaster dcr041009*/ /*XHTangle dcr20120120*/ if (dot) {/*add new dot to head of lists arrayed by type and severity */ which = type; dot->next = results[which][idx]; results[which][idx] = dot; } #ifdef OLDCODE } /* otherwise we forget about the dot */ #endif }/*saveDot()*/ /*}}}saveDot_________________________________________________________________*/ /*{{{newDot()*****************************************************************/ dotNode* newDot(atom *src, atom* targ, point3d *loc, point3d *spike, int ovrlaptype, float gap, char masterchr, float angle) { /*ptmaster dcr041009, XHTangle dcr20120120*/ dotNode* d = NULL; d = (dotNode *)malloc(sizeof(dotNode)); if (d) { d->next = NULL; d->a = src; d->t = targ; d->loc = *loc; d->spike = *spike; d->type = ovrlaptype; /* -1 bump, 0 touch, +1 H bond */ d->gap = gap; d->ptmaster = masterchr; /*dcr041009*/ d->angle = angle; /* XHTangle dcr20120120*/ } else warn("could not allocate space for new dot"); return d; } /*}}}newDot()________________________________________________________________*/ /*{{{INRANGEINLINE()**********************************************************/ #define INRANGEINLINE(a, b, lim) \ (((((a).x-(b).x)*((a).x-(b).x))+\ (((a).y-(b).y)*((a).y-(b).y))+\ (((a).z-(b).z)*((a).z-(b).z))) <= ((lim)*(lim))) /*returns true when distance-sq between a and b <= limit-sq */ /*}}}INRANGEINLINE()_________________________________________________________*/ /*{{{examineOneDotEach()****** called from genDotIntersect() LOneDotEach *****/ void examineOneDotEach(atom *src, int type, atom *scratch, pointSet dots[], float probeRad, float spikelen, int targFlg, dotNode *results[][NODEWIDTH], atom *allMainAtoms) /*allMainAtoms20120120*/ { /*called from genDotIntersect() with same names except scratch==atomList2 */ /*cloned and modified from examineDots to do only the OneDotEach option 111205*/ /*111013 switch to dump details... 111202 now input param: -dotdump */ atom *holdsrc; atom *holdcause; int holdtype; point3d holddotloc; point3d holdspikeloc; int holdovrlaptype; float holdmingap = 999.9f; char holdptmaster; int Lholdingdot = FALSE; /*111021dcr src and targ can have diff patterns*/ int Lbetweenatoms = FALSE; /*111202dcr*/ atom *targ = NULL, *cause = NULL; atom *extratarg = NULL; /*extra outer loop over neighbors, 111013*/ int i = 0, nearpt = 0, within = 0, ok = 0, maskHB = 0; int ndots = 1; /*for only One Dot for Each neighbor LOneDotEach 111008*/ int skipthisMCMC = 0, isaHB = 0, tooCloseHB = 0, ovrlaptype = 0; float gap = 0.0, mingap = 0.0, sl = 0.0, hbondbumpgap = 0.0; float srcCauseDist = 0.0; /*check for disaster overlap condition*/ float hbcutoff = 0.0, targVDW = 0.0; point3d targetloc, dotloc, dotvect, exploc, spikeloc; point3d gapvect, gaploc; point3d perploc; /*LOneDotEach occluder perpendicular to src-targ 111018*/ pointSet *srcDots = NULL; char ptmaster = ' '; /*pointmaster character dcr041009*/ int idx = 0; /*dcr041020*/ atom *parent = NULL; /*dcr20120120*/ /*int LXHvector = TRUE;*/ /*dcr20120120 eventually commandline param*/ float XHdistance = 99.9f; /*dcr20120120*/ float XHTangle = 0.0; /*dcr20120120*/ if (src->elem == ignoreAtom) { return; } /*050119 but there seems no way to ever set an atom->elem = ignoreAtom !?*/ if (OccupancyCriteria && src->occ < OccupancyCutoff) /*Occ.Crit. 060831*/ {/*050119*/ return; } cause = NULL; /* who is responsible for the dot */ srcDots = &(dots[src->elem]); if (src->elem == atomC && isCarbonylAtom(src)) {srcDots = &COdots;} maskHB = 0; if(src->props & DONOR_PROP) { maskHB = ACCEPTOR_PROP; if(src->props & ACCEPTOR_PROP) { /* ambig */ maskHB |= DONOR_PROP; } } else if (src->props & ACCEPTOR_PROP) { maskHB = DONOR_PROP; } /*111013 original mode: for each src atom (each call to this subroutine) */ /* place dot ball at src, */ /* loop over all dots, for each dot loop over neighbors*/ /*OneDotEach mode: for each src atom (each call to this subroutine) */ /* Outer loop over all neighbors, each in turn taken to cause a dot*/ /* for each neighbor, calculate dot at src surface in line with neighbor */ /* pass dot to Inner loop for this single dot,*/ /* then loop over neighbors to see if any occlude the line-of-sight src--cause */ /* I.E. wrap orig "Inner" dot loop with extra "Outer" neighbor loop, */ /* alternate modes do only one pass at one or the other loop. */ if(Ldotdump) fprintf(stderr,"\n\n******* SOURCE Atom at src:%7.3f,%7.3f,%7.3f, name %s %s %s SOURCE\n",src->loc.x,src->loc.y,src->loc.z,src->atomname,src->r->resname,src->r->Hy36resno); /* */ /*{{{OUTER LOOP over target atoms*/ for(extratarg = scratch; extratarg; extratarg = extratarg->scratch) /*OUTER LOOP*/ {/* OUTER LOOP over target atoms */ /* scan over neighbors and calc in-sight-line dot position for each */ if(Ldotdump) fprintf(stderr,"[src %s %s %s] TARGET:%7.3f,%7.3f,%7.3f, name %s %s %s TARGET\n",src->atomname,src->r->resname,src->r->Hy36resno,extratarg->loc.x,extratarg->loc.y,extratarg->loc.z,extratarg->atomname,extratarg->r->resname,extratarg->r->Hy36resno); /* */ /*operate on a single dot from outer loop */ /*first get src to target vector, call this working dot vector*/ v3sub(&(extratarg->loc), &(src->loc), &dotvect); /*then scale to src vdW to get src ctr to close point vector*/ v3scale(&dotvect, src->radius); /*the actual dot vector*/ /*add this to src position to get close point position in 3 space*/ v3add(&(src->loc), &dotvect, &dotloc);/*dotlocation on line to target*/ /*save actual dotvector, scaled later to get dot on other sde of gap*/ gapvect.x = dotvect.x; gapvect.y = dotvect.y; gapvect.z = dotvect.z; /*scale this out by probeRadius to get exploc */ v3scale(&dotvect, src->radius + probeRad); v3add(&(src->loc), &dotvect, &exploc); /*probeRad out from dotloc*/ if(Ldotdump) fprintf(stderr,"[src %s %s %s] DOT [%s %s %s] dotloc:%7.3f,%7.3f,%7.3f DOT\n",src->atomname,src->r->resname,src->r->Hy36resno,extratarg->atomname,extratarg->r->resname,extratarg->r->Hy36resno,dotloc.x,dotloc.y,dotloc.z); /* */ mingap = 999.9f; isaHB = tooCloseHB = FALSE; hbondbumpgap = 0.0; ok = FALSE; /*expect dot to be rejected*/ /*{{{attributes of dot wrt causal target: sets ok */ /*LOneDotEach only considers target atom that caused this dot*/ /* to assign attributes that might filter out this dot */ /*so at the end of this section, we will break out and then */ /*if the dot is of the right sort so far, we will scan over */ /*neighbors to see if this dot is occluded. */ targ = extratarg; /*override: evaluate causual target for this dot*/ hbcutoff = Min_regular_hb_cutoff; /* redefined later on */ targetloc = targ->loc; targVDW = targ->radius; #ifdef INLINE_FOR_SPEED nearpt = INRANGEINLINE(exploc, targetloc, probeRad+targVDW); #else nearpt = inRange(&exploc, &targetloc, probeRad+targVDW); #endif /*returns true when distance-sq between a and b <= limit-sq */ /*exploc is location of probe ctr with probe surface at dotloc*/ /*Fails erratically when atoms superimpose so distance back from a */ /*probe distance beyond surface of src*/ /*is a full probe radius + vdW to the center of target atom */ /*So full overlap disaster condition needs to be tested*/ /*before filter on nearpt value. */ if(Ldotdump) { fprintf(stderr,"nearpt=1 if within reach of probe\n"); fprintf(stderr," targ->mark>=1 if target within BONDED set of src\n"); fprintf(stderr,"nearpt %d, targ->mark %d [src: %s %s %s] %8.3f,%8.3f,%8.3f ; [targ %s %s %s] %8.3f,%8.3f,%8.3f ; exp:%8.3f,%8.3f,%8.3f\n",nearpt,targ->mark,src->atomname,src->r->resname,src->r->Hy36resno,src->loc.x,src->loc.y,src->loc.z,targ->atomname,targ->r->resname,targ->r->Hy36resno,targetloc.x,targetloc.y,targetloc.z,exploc.x,exploc.y,exploc.z); } skipthisMCMC = 0; /*just checking if this target atom acceptable for current src dot*/ if( DoMcMc == FALSE && ( src->props & MC_PROP) && (targ->props & MC_PROP) ) {/*potential MC-MC for us to skip this target atom*/ if (! ( ( src->props & HET_PROP) || (targ->props & HET_PROP) ) ) {/* neither can be a HET */ if (strcmp(src->r->chain, targ->r->chain) == 0) {/* must be the same chain*/ /*problem for NMR models ????*/ skipthisMCMC = 1; } } } /*now enter if...if else... chain of decisions...*/ /*drop through defaults to ok=0, dot is REJECTED*/ if (targ->elem == ignoreAtom) {/*NOP: condition not contribute to decision*/} else if (OccupancyCriteria && targ->occ < OccupancyCutoff) {/*050119 will not contribute */} /*Occ.Crit. 060831*/ else if (skipthisMCMC) {/*NOP: condition not contribute to decision*/} else if( DoWatWat == FALSE && ( src->props & WATER_PROP) && (targ->props & WATER_PROP)) { ; /* water will not contribute to decision when not being considered*/ } if(targ->mark == 0) {/*not flagged as a bonded atom*/ /*20111207dcr: unfortunately, can't distinguish between atoms*/ /* that by element type might be bonding and those in separate*/ /* groups that really are too close*/ /*overlap disaster: when src and target(==cause) atoms on top of each other*/ /*tests done from the edge of src to ctr of target are ill behaved here*/ /*double v3distance(point3d* a, point3d* b) */ /*{return(sqrt(v3distanceSq(a, b)));}*/ srcCauseDist = v3distance(&(extratarg->loc), &(src->loc)); if( (srcCauseDist <= (src->radius))||(srcCauseDist <= (extratarg->radius))) { nearpt = 2; /*disasterously near*/ ok = TRUE; /*really want to see this clash*/ cause = targ; /*targ gets reused, use cause as identifier of extratarg atom*/ if(Ldotdump) fprintf(stderr,"nearpt==2, atoms overlap!\n"); } } /*do next as plain if, rather than else if 20111207 dcr*/ if((nearpt==1) && (targ->mark == 0)) /*filter on nearpt value*/ {/* close but nonbonded, so might be an OK dot...*/ if(Ldotdump) fprintf(stderr,"close but nonbonded: calc gap...\n"); gap = gapSize(&dotloc, &targetloc, targVDW); /*111013*/ if(Ldotdump) fprintf(stderr,"[src %s %s %s] gap== %8.3f,mingap== %8.3f [etarg %s %s %s]\n",src->atomname,src->r->resname,src->r->Hy36resno,gap,mingap,targ->atomname,targ->r->resname,targ->r->Hy36resno); /*identify target with minimum gap*/ /*set properties*/ if (gap < mingap) /*otherwise ok wll remain FALSE*/ {/*if: gap < mingap*/ int bothCharged = FALSE; int chargeComplement = FALSE; if( ( src->props & (NEGATIVE_PROP|POSITIVE_PROP)) &&(targ->props & (NEGATIVE_PROP|POSITIVE_PROP)) ) { bothCharged = TRUE; chargeComplement = ( ( ( src->props & POSITIVE_PROP) && (targ->props & NEGATIVE_PROP)) || ( ( src->props & NEGATIVE_PROP) && (targ->props & POSITIVE_PROP)) ); } if( (targ->props & maskHB) /*test for H-BOND*/ && ( (! bothCharged) || chargeComplement) ) { hbcutoff = bothCharged ? Min_charged_hb_cutoff : Min_regular_hb_cutoff; if (gap < -hbcutoff) { tooCloseHB = TRUE; /* treat like a bump */ hbondbumpgap = gap; isaHB = TRUE; } else { /* hbond or contact */ isaHB = TRUE; tooCloseHB = FALSE; } } else {/* clash or contact, not an H-BOND */ /*NOTE: atomHOd : hb-only-dummy, phantom H atom*/ if (src->elem == atomHOd) { continue; } /*with for: loop*/ if (targ->elem == atomHOd) { continue; } /*with for: loop*/ isaHB = FALSE; } mingap = gap; /* update minimum gap */ ok = TRUE; cause = targ; /*targ is loop counter, about to go out of scope*/ }/*if: gap < mingap*/ }/*else if: close but nonbonded*/ if(nearpt == 2) /*disasterously near, i.e. on top of each other*/ { ok = TRUE; if(Ldotdump) fprintf(stderr,"OK by courtesy, nearpt==2, keep horrible clash\n"); } if(ok) { if(Ldotdump) fprintf(stderr,"OK so far ... (targ->mark== %d)\n",targ->mark); } else { if(Ldotdump) fprintf(stderr,"NOT OK: REJECT (targ->mark== %d)\n",targ->mark); } /*}}}attributes of dot wrt causal target: sets ok */ if(Ldotdump) fprintf(stderr,"Mark1 ok== %d, nearpt== %d\n",ok,nearpt); /*cause is atom identified as extratarg which invokes the dot on src atom*/ /*targ is current atom being considered, overloaded name: reused below */ /*as putative occluding 3rd atom from the surrounding neighbor atoms*/ /*Does OneDotEach always identify cause atom as a targFlg atom, need to check??*/ /*BUT probe has several ways to define sets of atoms, some as src, some as targ*/ if (ok && !(cause->flags & targFlg)) { /* drop non-target atom dot */ ok = FALSE; } if (ok && (nearpt==1)) {/*ok so far!: check for occluding atoms */ if(Ldotdump) fprintf(stderr,"DOT OK so far, check for OCCLUDING neighbors\n"); /*{{{scan over putative occluding neighbors*/ /*targ is an atom*, scratch is an atom*, atom* defined in abin.h */ /* atom* scratch is a member of atom*: which has member scratch...*/ /* so for-loop moves through all these atoms while an atom has a scratch*/ for(targ = scratch; targ; targ = targ->scratch) {/*for: scan over target atoms*/ if(LOneDotEach && (targ == extratarg))continue;/*to interate loop*/ /*Looking from src atom toward extratarg, check neighbors for occlusion*/ /*111013*/ if(Ldotdump) fprintf(stderr,"----------OCCLUSION test: [src %s %s %s]dot vs [neighbor %s %s %s elem==%d]\n",src->atomname,src->r->resname,src->r->Hy36resno,targ->atomname,targ->r->resname, targ->r->Hy36resno, targ->elem); /*elem 130427dcr*/ /* "targ" atoms here are neighbors potentially occluding*/ if(Lhonlyocclude && !isHatom(targ->elem) )continue;/*to interate loop 130427dcr*/ if (targ->mark && (targ->elem != atomHOd) &&!LOneDotEach)/*111122 ODE test*/ /*if (targ->mark && (targ->elem != atomHOd))*/ /*111122 ODE test*/ {/*bonded atom*/ /*NOTE: atomHOd : hb-only-dummy, phantom H atom*/ #ifdef INLINE_FOR_SPEED within = INRANGEINLINE(dotloc, (targ->loc), targ->radius); #else /*returns true when distance-sq between a and b <= limit-sq */ within = inRange(&dotloc, &(targ->loc), targ->radius); #endif if (within) { if(Ldotdump) fprintf(stderr," NEIGHBOR [%s %s %s] OCCLUDS src---targ interaction\n",targ->atomname, targ->r->resname,targ->r->Hy36resno); ok = FALSE; break; } else { if(Ldotdump) fprintf(stderr," NEIGHBOR [%s %s %s] BENIGN for src---targ interaction\n",targ->atomname, targ->r->resname,targ->r->Hy36resno); } }/*bonded atom*/ else if ((src->elem == atomHOd) && !(targ->props & ACCEPTOR_PROP)) {/*dummy HYDROGEN*/ /* atomHOd : hb-only-dummy, phantom H atom*/ /*eliminate if H? within non-Acceptor atom! WIERD*/ /*111122dcr maybe this is for dummy H on water: only make H-bonds????*/ #ifdef INLINE_FOR_SPEED within = INRANGEINLINE(dotloc, (targ->loc), targ->radius); #else /*returns true when distance-sq between a and b <= limit-sq */ within = inRange(&dotloc, &(targ->loc), targ->radius); #endif if (within) { if(Ldotdump) fprintf(stderr,"WIERD:HYDROGEN?[%s %s %s] OCCLUDS src---target\n",targ->atomname,targ->r->resname,targ->r->Hy36resno); ok = FALSE; break; } else { if(Ldotdump) fprintf(stderr,"WIERD:HYDROGEN?[%s %s %s] BENIGN src---target\n",targ->atomname,targ->r->resname,targ->r->Hy36resno); } }/*dummy HYDROGEN*/ /* atomHOd : hb-only-dummy, phantom H atom*/ else if(LOneDotEach) {/*LOneDotEach*/ {/*project putative occluding atom onto src--target line*/ /*drop perpendicular from putative occluder to src--target line*/ /*NOTE: one atom is src, the other is extratarg == target causes dot*/ /*putative occluding atom == "targ", i.e. scan all targ for each src*/ if(Ldotdump) fprintf(stderr,"check PERPENDICULAR distance of occluding atom to src---targ line\n"); /*get position on line defined by atoms where 3rd atom is perpendicular*/ perpendicular(&(src->loc),&(extratarg->loc),&(targ->loc),&perploc); /*now check distance of putative occluder (targ) to perploc*/ #ifdef INLINE_FOR_SPEED within = INRANGEINLINE(perploc, (targ->loc), targ->radius); #else within = inRange(&perploc, &(targ->loc), targ->radius); #endif if(within) {/*occluder within vdW of line defined by atoms (src and extratarg)*/ if(Ldotdump) fprintf(stderr,"occluding atom within it's own vdW of src---targ line\n"); /*TWO conditions can cause occlusion: */ /* 1:block gap */ /*close enough to occlude interaction if perploc within the gap*/ /*src loc -- perp loc -- target-causing-dot loc */ /* 2:eat dots */ /*3rd atom overlaps either one of the gap defining points*/ /* which can happen even if perploc not between atoms*/ /*TWO possible occluding situations: block gap or eat dots*/ /*need to know location of gap dot on target atom side*/ v3scale(&gapvect, src->radius + gap); v3add(&(src->loc), &gapvect, &gaploc); if( ( v3distance(&gaploc,&perploc) < fabs(gap) ) &&( v3distance(&dotloc,&perploc) < fabs(gap) ) ) {/*perp hits within gap so does occlude line of sight*/ if(Ldotdump) fprintf(stderr,"atom [%s %s %s] OCCLUDES src---targ line at %7.3f %7.3f %7.3f, DOT REJECTED\n" ,targ->atomname,targ->r->resname,targ->r->Hy36resno,perploc.x,perploc.y,perploc.z); ok = FALSE; break; }/*perp hits within gap so does occlude line of sight*/ /*check if neighbor close enough to occlude gapdots*/ /*need to know if perploc actually between atom centers*/ /*check that perpendicular drop point is between the two atoms*/ if( (v3distance(&perploc,&(src->loc)) < v3distance(&(src->loc),&(extratarg->loc))) &&(v3distance(&perploc,&(extratarg->loc)) < v3distance(&(src->loc),&(extratarg->loc)))) {/*perp hits line inbetween atoms*/ Lbetweenatoms = TRUE; }/*perp hits line inbetween atoms*/ else {Lbetweenatoms = FALSE;} /*now see if 3rd atom eats dot*/ if( ( (v3distance(&dotloc,&(targ->loc)) < targ->radius) ||(v3distance(&gaploc,&(targ->loc)) < targ->radius) ) &&Lbetweenatoms) {/*neighbor vdW engulfs gap dot(s),i.e. OCCLUDES*/ if(Ldotdump) fprintf(stderr,"atom [%s %s %s] OCCLUDES gap dots,distance: %7.3f & %7.3f,vs vdW %7.3f: DOT REJECTED\n" ,targ->atomname,targ->r->resname,targ->r->Hy36resno, v3distance(&dotloc,&(targ->loc)),v3distance(&gaploc,&(targ->loc)),targ->radius ); ok = FALSE; break; }/*neighbor vdW engulfs gap dot(s),i.e. OCCLUDES*/ else { if(Ldotdump) fprintf(stderr,"atom [%s %s %s] NOT occlude dots,distance: %7.3f & %7.3f,vs vdW %7.3f: DOT OK\n" ,targ->atomname,targ->r->resname,targ->r->Hy36resno, v3distance(&dotloc,&(targ->loc)),v3distance(&gaploc,&(targ->loc)),targ->radius ); } /*check if neighbor close enough to occlude gapdots*/ }/*occluder within vdW of line defined by atoms (src and extratarg)*/ else { if(Ldotdump) fprintf(stderr,"this neighbor atom too far away to occlude src---targ line\n"); } }/*project putative occluding atom onto src--target line*/ }/*LOneDotEach*/ }/*for: scan over target atoms*/ /*}}}scan over putative occluding neighbors*/ }/*ok so far?: check for occluding atoms */ if(Ldotdump) fprintf(stderr,"OK2 == %d, DOT to be outputed if ok\n",ok); if (ok) {/*dot is ok, atom that this dot hits is called cause */ /* Contact */ /* ovrlaptype : -1 bump, 0 touch, +1 H bond */ /* mingap originally minimum distance of src dot to cause atom*/ if(nearpt == 2) /*disasterously near, i.e. on top of each other*/ { sl = -((src->radius) + (extratarg->radius))/2.0f; mingap = 2.0f*LowGoodCut; /*force arbitrary bad clash*/ ovrlaptype =-1; /*disaster overlap == clash!! */ /*spikeloc calc from src->radius + sl, sl=spikelen*mingap,*/ /*spikelen = 0.50 OR input param*/ /*usually: sl = spikelen*mingap*/ } else if (mingap > 0.0) { sl = 0.0; if(LweakHbonds && isaHB) {ovrlaptype = 2;} /*20111215dcr*/ /* TODO SJ: Why is this 2*/ else {ovrlaptype = 0;} } else if (isaHB && tooCloseHB) { mingap = hbondbumpgap + hbcutoff; sl = spikelen*mingap; ovrlaptype =-1; /* Hbond which is too close is a clash */ } else if (isaHB) /* Hbond */ { sl = spikelen*mingap; ovrlaptype =+1; /* must test angle for Hbond in some cases */ if ( src->props & TEST_ACCEPT_ANGLE_PROP) { /* cause->loc vs src->aromaticRing; ovrlaptype =-1;?? */ } else if ( cause->props & TEST_ACCEPT_ANGLE_PROP) { /* src->loc vs cause->aromaticRing; ovrlaptype =-1;?? */ } if( src->props & CH_DONOR_PROP || cause->props & CH_DONOR_PROP) {/* CH..O type Hbond */ sl *= CHOHBfactor; /*down scale score for these types of Hbonds*/ } #ifdef JACKnANDREA if (writeHbonds) printf("{%s %s %s %d}P %f %f %f {%s %s %s %d} %f %f %f\n", src->atomname, src->r->resname, src->r->chain, src->r->resid, src->loc.x, src->loc.y, src->loc.z, cause->atomname, cause->r->resname, cause->r->chain, cause->r->resid, cause->loc.x, cause->loc.y, cause->loc.z); /*this makes a monster list of all possible H-bond dots as heavy-atom--H vecs*/ /*which must then be pared down to the unique atom-pair vectors*/ /*also includes all the atomOHd dummy phantoms which junks up the display*/ #endif } else /* Clash */ { sl = spikelen*mingap; ovrlaptype =-1; } if(Ldotdump) fprintf(stderr,"[src %s %s %s] OK2 passed, ovrlaptype == %d********************\n",src->atomname,src->r->resname,src->r->Hy36resno,ovrlaptype); /* */ /*NOTE: atomHOd : hb-only-dummy, phantom H atom*/ /* ovrlaptype : -1 bump, 0 touch, +1 H bond */ if((ovrlaptype == 1)||(src->elem != atomHOd && cause->elem != atomHOd)) { /* TODO SJ: Figure out what will happen is weak H bonds is true*/ v3scale(&dotvect, src->radius + sl); v3add(&(src->loc), &dotvect, &spikeloc); /* possibly limit contact dots... */ /* ovrlaptype : -1 bump, 0 touch, +1 H bond */ /*kissEdge2bullsEye seems not responsible for loss of very close clashes*/ if ( ovrlaptype != 0 || LimitDots == FALSE || dot2srcCenter(&dotloc, src, cause) <= kissEdge2bullsEye(src->radius, cause->radius, probeRad)) {/*could use this dot/spike, */ /*saveDot decisions moved here to examineDots 041020 */ /*overlaptype: -1 bump, 0 touch, +1 H bond */ idx = dotClassIndex(ovrlaptype, mingap); /* OLD: idx: 0 wide contact,1 close contact,2 small overlap,3 bad overlap,4 H-bonds*/ /* New SJ 04/21/2015: idx: 0 wide contact,1 close contact,2 weak H-bonds, 3 small overlap,4 bad overlap, 5 worse overlaps, 6 H-bonds*/ /*OnlyBadOut option dcr041010*/ /* ovrlaptype : -1 bump, 0 touch, +1 H bond */ if( (!OnlyBadOut && OutputHBs && ovrlaptype > 0) ||(!OnlyBadOut && OutputClashes && ovrlaptype < 0) ||( OnlyBadOut && OutputClashes && /*idx==3*/ (idx ==4 || idx == 5))/*Bad Clash dcr041020*/ /* SJ 04/21/2015 changed idx according to new values*/ ||(!OnlyBadOut && OutputVDWs && ovrlaptype == 0)) { /*logicals allow the overlaptype of this dot*/ /*now derive ptmaster*/ if( (src->props & MC_PROP) && (cause->props & MC_PROP) ) {/*potential MC-MC to flag*/ if( !(( src->props & HET_PROP) || (cause->props & HET_PROP)) ) {/*neither can be a HET, what wierd condition is avoided?*/ ptmaster = 'M'; /*McMc interaction*/ mcmccont[0][idx]++; //mcmccont[0][5]++; mcmccont[0][7]++; /* SJ 04/21/2015 changing the second dimension to 7 to match new initialization*/ } } else if( (src->props & SC_PROP) && (cause->props & SC_PROP) ) {/*potential SC-SC to flag*/ if( !(( src->props & HET_PROP) || (cause->props & HET_PROP)) ) {/*neither can be a HET, what wierd condition is avoided?*/ ptmaster = 'S'; /*ScSc interaction*/ scsccont[0][idx]++; //scsccont[0][5]++; scsccont[0][7]++;/* SJ 04/21/2015 changing the second dimension to 7 to match new initialization*/ } } else if( ((src->props & SC_PROP) && (cause->props & MC_PROP)) ||((src->props & MC_PROP) && (cause->props & SC_PROP))) {/*potential MC-SC or SC-MC to flag*/ if( !(( src->props & HET_PROP) || (cause->props & HET_PROP)) ) {/*neither can be a HET, what wierd condition is avoided?*/ ptmaster = 'P'; /*McSc or ScMc interaction*/ mcsccont[0][idx]++; //mcsccont[0][5]++; mcsccont[0][7]++;/* SJ 04/15/2015 changing the second dimension to 7 to match new initialization*/ } } else { ptmaster = 'O'; /*Oh for Other, general default*/ othrcont[0][idx]++; //othrcont[0][5]++; othrcont[0][7]++;/* SJ 04/15/2015 changing the second dimension to 7 to match new initialization*/ } /*111013*/ if(Ldotdump) fprintf(stderr,"[src %s %s %s] :idx==%d: [neigh %s %s %s]*** saveDot:: gap %8.3f, mingap %8.3f, holdmingap %8.3f, dotloc: %7.3f,%7.3f,%7.3f; spikeloc: %7.3f,%7.3f,%7.3f \n",src->atomname,src->r->resname,src->r->Hy36resno,idx,cause->atomname, cause->r->resname,cause->r->Hy36resno,gap,mingap,holdmingap, dotloc.x,dotloc.y,dotloc.z,spikeloc.x,spikeloc.y,spikeloc.z); if(Lnearest) { if(mingap < holdmingap) { holdsrc = src; holdcause = cause; holdtype = type; holddotloc.x = dotloc.x; holddotloc.y = dotloc.y; holddotloc.z = dotloc.z; holdspikeloc.x = spikeloc.x; holdspikeloc.y = spikeloc.y; holdspikeloc.z = spikeloc.z; holdovrlaptype = ovrlaptype; holdmingap = mingap; holdptmaster = ptmaster; Lholdingdot = TRUE; if(Ldotdump) fprintf(stderr,"reset holdings: holdmingap %8.3f\n",holdmingap); } } else { if(Ldotdump) fprintf(stderr,"saveDot CALL+++++++++++++++++++++++\n"); saveDot(src, cause, type, &dotloc, &spikeloc, results, ovrlaptype, mingap, ptmaster,XHTangle); /*dcr041009*/ } }/*logicals allow the overlaptype of this dot*/ }/*could use this dot/spike, */ } }/*dot is ok, atom that this dot hits is called cause */ /*operate on a single dot from outer loop */ /*111013*/ if(Ldotdump) fprintf(stderr,"+++++++++++END of pass: loop on extratarg, finished with %s %s %s\n\n",extratarg->atomname,extratarg->r->resname,extratarg->r->Hy36resno); /* */ }/* OUTER LOOP over target atoms */ /*}}}OUTER LOOP over target atoms*/ if(Lnearest && Lholdingdot) {/*Lnearest mode and has an OK dot which is nearest that is allowed by targPat*/ ok = TRUE; if(Lfilternearest) /*behavior if extraPat==NULL ???? */ { if(!matchPat(holdcause, extraPat)) /*use 3rd Pattern 111022dcr*/ {/*src deemed OK before call examineDots, unless cause OK, ignore atom*/ ok = FALSE; } } if(ok && Lfilterocc1) /*final filter exclude any with occ < 1 20120120dcr*/ { if(holdsrc->occ < 1 || holdcause->occ < 1) { ok = FALSE; fprintf(stderr,"EXCLUDE[src %s OCC:%4.2f, %s %s] :: [cause %s OCC:%4.2f, %s %s]\n",holdsrc->atomname,holdsrc->occ,holdsrc->r->resname,holdsrc->r->Hy36resno ,holdcause->atomname,holdcause->occ,holdcause->r->resname,holdcause->r->Hy36resno); } } if(ok && Lfilternoalt) /*final filter exclude any alts 20120120dcr*/ { if(holdsrc->altConf != ' ' || holdcause->altConf != ' ') { ok = FALSE; fprintf(stderr,"EXCLUDE[src %s ALT:%c, %s %s] :: [cause %s ALT:%c, %s %s]\n" ,holdsrc->atomname,holdsrc->altConf,holdsrc->r->resname,holdsrc->r->Hy36resno ,holdcause->atomname,holdcause->altConf,holdcause->r->resname,holdcause->r->Hy36resno); } } if(ok) { if(Ldotdump) fprintf(stderr,"[src %s %s %s] :: [neigh %s %s %s]*YES**about to saveDot:: mingap %8.3f HOLD\n",holdsrc->atomname,holdsrc->r->resname,holdsrc->r->Hy36resno,holdcause->atomname,holdcause->r->resname,holdcause->r->Hy36resno,holdmingap); if(LXHvector && holdsrc->elem == atomH) { for(parent = allMainAtoms; parent; parent = parent->next) /*dcr20120120*/ {/*for: loop over all main atoms taking each in turn as the parent-atom*/ XHdistance = v3distance(&(parent->loc), &(src->loc)); if(XHdistance < 1.2 && XHdistance > 0.1 ) {/*parent atom for this Hydrogen*/ XHTangle = v3angle(&(parent->loc), &(src->loc), &(holdcause->loc)); XHTangle = XHTangle*360/(2*3.14159265f); if(Ldotdump) fprintf(stderr,"[src %s %s %s] :: [parent %s %s %s]*YES**about to saveDot:: mingap %8.3f, XHdistance: %8.3f, XHTangle: %8.3f HOLD\n",holdsrc->atomname,holdsrc->r->resname,holdsrc->r->Hy36resno,parent->atomname,parent->r->resname,parent->r->Hy36resno,holdmingap,XHdistance,XHTangle); /*save the XHvector and XHTangle construction*/ /* saveDot(holdsrc, parent, holdtype, &parent->loc, &holdsrc->loc, results, -1, holdmingap, holdptmaster, XHTangle);*/ } } } else {XHTangle = 0.0;} /*now save the real dot*/ saveDot(holdsrc, holdcause, holdtype, &holddotloc, &holdspikeloc, results, holdovrlaptype, holdmingap, holdptmaster,XHTangle); } }/*Lnearest mode and has an OK dot*/ }/*examineOneDotEach()*/ /*}}}examineOneDotEach()_____________________________________________________*/ /*{{{examineDots()*********** called from genDotIntersect() ******************/ void examineDots(atom *src, int type, atom *scratch, pointSet dots[], float probeRad, float spikelen, int targFlg, dotNode *results[][NODEWIDTH]) { /*called from genDotIntersect() with same names except scratch==atomList2 */ atom *targ = NULL, *cause = NULL; int i = 0, nearpt = 0, within = 0, ok = 0, maskHB = 0; int skipthisMCMC = 0, isaHB = 0, tooCloseHB = 0, ovrlaptype = 0; float gap = 0.0, mingap = 0.0, sl = 0.0, hbondbumpgap = 0.0; float hbcutoff = 0.0, targVDW = 0.0; point3d targetloc, dotloc, dotvect, exploc, spikeloc; pointSet *srcDots = NULL; char ptmaster = ' '; /*pointmaster character dcr041009*/ int idx = 0; /*dcr041020*/ if (src->elem == ignoreAtom) { return; } /*050119 but there seems no way to ever set an atom->elem = ignoreAtom !?*/ if (OccupancyCriteria && src->occ < OccupancyCutoff) /*Occ.Crit. 060831*/ {/*050119*/ return; } cause = NULL; /* who is responsible for the dot */ srcDots = &(dots[src->elem]); if (src->elem == atomC && isCarbonylAtom(src)) { srcDots = &COdots; } maskHB = 0; if(src->props & DONOR_PROP) { maskHB = ACCEPTOR_PROP; if(src->props & ACCEPTOR_PROP) { /* ambig */ maskHB |= DONOR_PROP; } } else if (src->props & ACCEPTOR_PROP) { maskHB = DONOR_PROP; } for(i=0; i < srcDots->n; i++) {/*loop over dots on src dot-ball */ dotvect = srcDots->p[i]; v3add(&(src->loc), &dotvect, &dotloc); v3scale(&dotvect, src->radius + probeRad); v3add(&(src->loc), &dotvect, &exploc); ok = FALSE; mingap = 999.9f; isaHB = tooCloseHB = FALSE; hbondbumpgap = 0.0; /*targ is an atom*, scratch is an atom*, atom* defined in abin.h */ /* atom* scratch is a member of atom*: which has member scratch...*/ /* so for-loop moves through all these atoms while an atom has a scratch*/ for(targ = scratch; targ; targ = targ->scratch) {/*for: loop over target atoms: set ok=1 for some subset*/ hbcutoff = Min_regular_hb_cutoff; /* redefined later on */ targetloc = targ->loc; targVDW = targ->radius; #ifdef INLINE_FOR_SPEED nearpt = INRANGEINLINE(exploc, targetloc, probeRad+targVDW); #else nearpt = inRange(&exploc, &targetloc, probeRad+targVDW); #endif /*returns true when distance-sq between a and b <= limit-sq */ /*...in subroutine examineDots()...*/ skipthisMCMC = 0; /*just checking if this target atom acceptable for current src dot*/ /*(premature to assign pointmasters when target not yet determined*/ if( DoMcMc == FALSE && ( src->props & MC_PROP) && (targ->props & MC_PROP) ) { /*potential MC-MC for us to skip this target atom*/ if (! ( ( src->props & HET_PROP) || (targ->props & HET_PROP) ) ) {/* neither can be a HET */ if (strcmp(src->r->chain, targ->r->chain) == 0) {/* must be the same chain*/ /*problem for NMR models ????*/ skipthisMCMC = 1; } } } /*now enter if...if else... chain of decisions...*/ if (targ->elem == ignoreAtom) {/* will not contribute */} else if (OccupancyCriteria && targ->occ < OccupancyCutoff) {/*050119 will not contribute */} /*Occ.Crit. 060831*/ else if (skipthisMCMC) {/* will not contribute */} else if( DoWatWat == FALSE && ( src->props & WATER_PROP) && (targ->props & WATER_PROP)) { ; /* water will not contribute when not being considered*/ } /*Lautobondrot 050111 investigation point */ else if (nearpt && (targ->mark == 0)) {/*else if: close but nonbonded*/ gap = gapSize(&dotloc, &targetloc, targVDW); /*identify target with minimum gap*/ if (gap < mingap) {/*if: gap < mingap*/ int bothCharged = FALSE; int chargeComplement = FALSE; if( ( src->props & (NEGATIVE_PROP|POSITIVE_PROP)) &&(targ->props & (NEGATIVE_PROP|POSITIVE_PROP)) ) { bothCharged = TRUE; chargeComplement = ( ( ( src->props & POSITIVE_PROP) && (targ->props & NEGATIVE_PROP)) || ( ( src->props & NEGATIVE_PROP) && (targ->props & POSITIVE_PROP)) ); } if( (targ->props & maskHB) && ( (! bothCharged) || chargeComplement) ) { hbcutoff = bothCharged ? Min_charged_hb_cutoff : Min_regular_hb_cutoff; if (gap < -hbcutoff) { tooCloseHB = TRUE; /* treat like a bump */ hbondbumpgap = gap; isaHB = TRUE; } else { /* hbond or contact */ isaHB = TRUE; tooCloseHB = FALSE; } } else { /* clash or contact */ /*NOTE: atomHOd : hb-only-dummy, phantom H atom*/ if (src->elem == atomHOd) { continue; } /*with for: loop*/ if (targ->elem == atomHOd) { continue; } /*with for: loop*/ isaHB = FALSE; } mingap = gap; /* update minimum gap */ ok = TRUE; cause = targ; }/*if: gap < mingap*/ }/*else if: close but nonbonded*/ }/*for: loop over target atoms: set ok=1 for some subset*/ if (ok && !(cause->flags & targFlg)) { /* drop non-target atom dot */ ok = FALSE; } if (ok) {/*ok: apply bonded atom dot filters */ /*targ is an atom*, scratch is an atom*, atom* defined in abin.h */ /* atom* scratch is a member of atom*: which has member scratch...*/ /* so for-loop moves through all these atoms while an atom has a scratch*/ for(targ = scratch; targ; targ = targ->scratch) {/*for: scan over target atoms*/ /* eliminate dot if within bonded atom! */ if (targ->mark && (targ->elem != atomHOd)) { /*NOTE: atomHOd : hb-only-dummy, phantom H atom*/ #ifdef INLINE_FOR_SPEED within = INRANGEINLINE(dotloc, (targ->loc), targ->radius); #else within = inRange(&dotloc, &(targ->loc), targ->radius); #endif /*returns true when distance-sq between a and b <= limit-sq */ /*...in subroutine examineDots()...*/ if (within) { ok = FALSE; break; } } else if ((src->elem == atomHOd) && !(targ->props & ACCEPTOR_PROP)) { /*NOTE: atomHOd : hb-only-dummy, phantom H atom*/ /* eliminate if H? within non-Acceptor atom! */ #ifdef INLINE_FOR_SPEED within = INRANGEINLINE(dotloc, (targ->loc), targ->radius); #else within = inRange(&dotloc, &(targ->loc), targ->radius); #endif /*returns true when distance-sq between a and b <= limit-sq */ /*...in subroutine examineDots()...*/ if (within) { ok = FALSE; break; } } }/*for: scan over target atoms*/ }/*ok: apply bonded atom dot filters */ if (ok) {/*dot is ok, atom that this dot hits is called cause */ /* Contact */ /* ovrlaptype : -1 bump, 0 touch, +1 H bond */ /*dcr?: mingap seems to do what ???? 050111*/ /* 04/10/2015 - SJ, including weak H bonds in overlaptype +1, if LweakHbonds is true*/ if (mingap > 0.0) { sl = 0.0; /* TODO: SJ not sure waht this sl is. Find out and change accordingly*/ if(LweakHbonds && isaHB) /* 04/10/2015 ovrlaptype is H bond if detecting weak H bonds*/ ovrlaptype = 1; else ovrlaptype = 0; } else if (isaHB && tooCloseHB) { mingap = hbondbumpgap + hbcutoff; sl = spikelen*mingap; ovrlaptype =-1; /* Hbond which is too close is a clash */ } else if (isaHB) /* Hbond */ { sl = spikelen*mingap; ovrlaptype =+1; /* must test angle for Hbond in some cases */ if ( src->props & TEST_ACCEPT_ANGLE_PROP) { /* cause->loc vs src->aromaticRing; ovrlaptype =-1;?? */ } else if ( cause->props & TEST_ACCEPT_ANGLE_PROP) { /* src->loc vs cause->aromaticRing; ovrlaptype =-1;?? */ } if( src->props & CH_DONOR_PROP || cause->props & CH_DONOR_PROP) {/* CH..O type Hbond */ sl *= CHOHBfactor; /*down scale score for these types of Hbonds*/ } #ifdef JACKnANDREA if (writeHbonds) printf("{%s %s %s %d}P %f %f %f {%s %s %s %d} %f %f %f\n", src->atomname, src->r->resname, src->r->chain, src->r->resid, src->loc.x, src->loc.y, src->loc.z, cause->atomname, cause->r->resname, cause->r->chain, cause->r->resid, cause->loc.x, cause->loc.y, cause->loc.z); /*this makes a monster list of all possible H-bond dots as heavy-atom--H vecs*/ /*which must then be pared down to the unique atom-pair vectors*/ /*also includes all the atomOHd dummy phantoms which junks up the display*/ #endif } else /* Clash */ { sl = spikelen*mingap; ovrlaptype =-1; } /*NOTE: atomHOd : hb-only-dummy, phantom H atom*/ /* ovrlaptype : -1 bump, 0 touch, +1 H bond */ if( (ovrlaptype == 1)/* TODO: SJ Have to figue out what will happen here is LweakHbonds is true*/ ||(src->elem != atomHOd && cause->elem != atomHOd) ) { v3scale(&dotvect, src->radius + sl); v3add(&(src->loc), &dotvect, &spikeloc); /* possibly limit contact dots... */ /* ovrlaptype : -1 bump, 0 touch, +1 H bond */ /*kissEdge2bullsEye seems not responsible for loss of very close clashes*/ if ( ovrlaptype != 0 || LimitDots == FALSE || dot2srcCenter(&dotloc, src, cause) <= kissEdge2bullsEye(src->radius, cause->radius, probeRad)) {/*could use this dot/spike, */ /*saveDot decisions moved here to examineDots 041020 */ /*overlaptype: -1 bump, 0 touch, +1 H bond */ idx = dotClassIndex(ovrlaptype, mingap); /*OLD idx: 0 wide contact,1 close contact,2 small overlap,3 bad overlap,4 H-bonds*/ /* new 04/10/2015 SJ: 0 wide contact, 1 close contact, 2 weak H bonds, 3 small overlap, 4 bad overlap, 5 worse overlap, 6 H-bonds*/ /*OnlyBadOut option dcr041010*/ /* ovrlaptype : -1 bump, 0 touch, +1 H bond */ if( (!OnlyBadOut && OutputHBs && ovrlaptype > 0) ||(!OnlyBadOut && OutputClashes && ovrlaptype < 0) ||( OnlyBadOut && OutputClashes && /*idx==3 SJ*/ (idx == 4 || idx == 5))/*Bad Clash dcr041020*/ /* idx value changes according to new above, SJ*/ ||(!OnlyBadOut && OutputVDWs && ovrlaptype == 0)) { /*logicals allow the overlaptype of this dot*/ /*now derive ptmaster*/ if( (src->props & MC_PROP) && (cause->props & MC_PROP) ) {/*potential MC-MC to flag*/ if( !(( src->props & HET_PROP) || (cause->props & HET_PROP)) ) {/*neither can be a HET, what wierd condition is avoided?*/ ptmaster = 'M'; /*McMc interaction*/ mcmccont[0][idx]++; //mcmccont[0][5]++; mcmccont[0][7]++; /* SJ 04/15/2015 changing the second dimension to 7 to match new initialization*/ } } else if( (src->props & SC_PROP) && (cause->props & SC_PROP) ) {/*potential SC-SC to flag*/ if( !(( src->props & HET_PROP) || (cause->props & HET_PROP)) ) {/*neither can be a HET, what wierd condition is avoided?*/ ptmaster = 'S'; /*ScSc interaction*/ scsccont[0][idx]++; //scsccont[0][5]++; scsccont[0][7]++;/* SJ 04/15/2015 changing the second dimension to 7 to match new initialization*/ } } else if( ((src->props & SC_PROP) && (cause->props & MC_PROP)) ||((src->props & MC_PROP) && (cause->props & SC_PROP))) {/*potential MC-SC or SC-MC to flag*/ if( !(( src->props & HET_PROP) || (cause->props & HET_PROP)) ) {/*neither can be a HET, what wierd condition is avoided?*/ ptmaster = 'P'; /*McSc or ScMc interaction*/ mcsccont[0][idx]++; //mcsccont[0][5]++; mcsccont[0][7]++;/* SJ 04/15/2015 changing the second dimension to 7 to match new initialization*/ } } else { ptmaster = 'O'; /*Oh for Other, general default*/ othrcont[0][idx]++; //othrcont[0][5]++; othrcont[0][7]++;/* SJ 04/15/2015 changing the second dimension to 7 to match new initialization*/ } saveDot(src, cause, type, &dotloc, &spikeloc, results, ovrlaptype, mingap, ptmaster,0); /* dcr041009 */ /*XHTangle dummy dcr20120120*/ }/*logicals allow the overlaptype of this dot*/ }/*could use this dot/spike, */ } }/*dot is ok, atom that this dot hits is called cause */ }/*loop over dots on src dot-ball */ }/*examineDots()*/ /*}}}examineDots()___________________________________________________________*/ /*{{{markBonds()**************************************************************/ /*markBonds() decides which atoms are bonded to each other and thus */ /*decides under what conditions close atoms will NOT be considered to clash*/ /*this seems to be the only subroutine that uses atom->bondedto information*/ /*for a given src atom: mark all atoms that are close enough to bond and are */ /*allowed to bond and extend this to thier neighbors through Maxbonded bonds*/ /*including marking the src atom itself if such a Maxbonded-loop returns to it*/ void markBonds(atom *src, atom *neighbors, int distcount, int max) {/*markBonds reformated 041112*/ atom *targ = NULL, *theHatom = NULL, *theOtherAtom = NULL; int nearpt, tooclose, isAanH, isBanH, letItBond = FALSE; /*when making surface dots, max==1 */ /*when looking for contacts and clashes, max==Maxbonded */ for(targ = neighbors; targ; targ = targ->scratch) {/*for(loop over neighbors)*/ /*first pass look at neighbors of external calling atom: targ->mark == 0*/ /*later look at neighbors of neighbors: targ->mark > max */ /* until Maxbonded neighbor when targ->mark == max */ if (targ->mark == 0 || targ->mark > distcount) {/*if(targ->mark == 0 || targ->mark > distcount) */ /* note: when we are relying on the StdBond table */ /* to figure out which bonds are allowed */ /* we can be more generous with the distance */ /* cuttofs to permit more oddball distortions */ /*dcr: so it seems that with UseStdBond 2*COVRADFUDGE is added */ /*atomprops.h/defines COVRADFUDGE 0.2 050111 */ #ifdef INLINE_FOR_SPEED nearpt = INRANGEINLINE((src->loc), (targ->loc), src->covRad + targ->covRad + COVRADFUDGE + (UseStdBond ? COVRADFUDGE : 0.0) ); tooclose = INRANGEINLINE((src->loc), (targ->loc), src->covRad + targ->covRad - 0.6 - (UseStdBond ? COVRADFUDGE : 0.0) ); #else nearpt = inRange(&(src->loc), &(targ->loc), src->covRad + targ->covRad + COVRADFUDGE + (UseStdBond ? COVRADFUDGE : 0.0) ); tooclose = inRange(&(src->loc), &(targ->loc), src->covRad + targ->covRad - 0.6 - (UseStdBond ? COVRADFUDGE : 0.0) ); #endif /*returns true when distance-sq between a and b <= limit-sq */ /*050111 tooclose= does NOT seem to be the call where */ /*atoms that move too close together suddenly stop showing bad clashes ! */ /*i.e. force tooclose = 0 here does not affect very close clash cut-off*/ /* and indeed atoms not bonded because of tooclose should show clashes !? */ /*but actual close distance of cut-off seems to be dependent on UseStdBond */ /*so since this is the only place where COVRADFUDGE is used, one suspects */ /*that this markBonds() is responsible for showing as bonded those atoms */ /*that fail to show clashes as they get near to each other*/ /*However, although this markBonds seems to check each atom pair many times, */ /*it does look as if it declares the letItBond flag correctly even when close!*/ /*050111 BUT this seems where all bonding info is calc, so maybe letItBond is */ /*not interpreted correctly by this analysis...*/ /*050117 specific case of ala_dipeptide at phi==0,psi==0, tau<= 102.25 seems */ /*the N of res 2 is just within bonding criteria of C of res 0 and with*/ /*UseStdBond N -- C bonds are allowed between atoms of different residues*/ /*thus the covalent neighbors of that N and that C, e.g. NH of 2 and O of 0 */ /*are also considered to NOT clash! */ /*For autobondrot and mage/probe strict succession would cure this...*/ /*BUT one does need to allow for deletions in residue order!*/ /*SO this is a FEATURE, NOT a BUG */ /*and the work-around is to restrict the tau deviation to what is more */ /*reasonable for tau anyway of -8, thus -5 to +10 around ideal 111.1 */ /*might be the best range for getting a feel for allowable phi,psi regions*/ /* conditions for allowing a bond to be formed... */ if (nearpt && (! tooclose) ) {/*close enough for a bond*/ isAanH = isHatom(src->elem); isBanH = isHatom(targ->elem); #ifdef OLD_BONDING_CALC if( ( ! isAanH && ! isBanH ) /* neither atom is a hydrogen */ ||( (! isAanH || ! isBanH) /*or one is an H */ && ATOMS_IN_THE_SAME_RES(src, targ) /* and both are in same*/ /* res with parent name*/ &&( (theHatom->bondedto == NULL) /*and either unknown */ or as expected*/ || strstr(theHatom->bondedto, theOtherAtom->atomname) ) ) ) { targ->mark = distcount; if (distcount < max) { markBonds(targ, neighbors, distcount+1, max); } } #else /*NOT OLD_BONDING_CALC*/ letItBond = FALSE; /* we start out skeptical */ if ((! isAanH) && (! isBanH)) {/* neither atom is a hydrogen - both are heavy atoms */ if (UseStdBond) {/* if we are allowing only "standard" bonding */ if (! ATOMS_IN_THE_SAME_RES(src, targ)) { /* and are in diff res */ if ((src->props & MC_PROP) && (targ->props & MC_PROP)) {/* both are mc in diff res - */ /*only specific bonding patterns allowed */ /*strcmp(a,b) returns 0 when strings == (a==b) */ /*strstr returns 1 when any str in colon delineated*/ /* list matches reference str */ letItBond = ( !strcmp(" N ", src->atomname) && !strcmp(" C ",targ->atomname) ) || ( !strcmp(" N ",targ->atomname) && !strcmp(" C ", src->atomname) ) || ( !strcmp(" P ", src->atomname) && strstr(" O3*: O3'", targ->atomname) ) || ( !strcmp(" P ",targ->atomname) && strstr(" O3*: O3'", src->atomname) ); } else if( !(src->props & MC_PROP) && !(targ->props & MC_PROP) ) { /* both are sc in diff res they must be CYS-CYS SS */ letItBond = ( !strcmp(src->r->resname, "CYS") && !strcmp(src->atomname, " SG ") && !strcmp(targ->r->resname, "CYS") && !strcmp(targ->atomname, " SG ") ); } /*050121 remove dcr041112 code insertion */ /*inserted test to allow mage/probe fitting of a*/ /*particular sidechain in context of an NMR ensemble*/ /*Lmodeltest not needed, mage sends model # when needed*/ }/* and are in diff res */ else {/* both heavy atoms in same res */ /*strstr returns 1 when any str in colon delineated*/ /* list matches reference str */ #ifdef STRICT_CONN_TABLE /*050111 seems to be not defined*/ /* strict - heavy atom must be as exactly expected */ letItBond = ( (src->bondedto != NULL) && strstr(src->bondedto, targ->atomname) ); #else /*NOT STRICT_CONN_TABLE*/ /* heavy atom either unknown or as expected */ letItBond = ( (src->bondedto == NULL) || strstr(src->bondedto, targ->atomname) ); #endif /*def? STRICT_CONN_TABLE*/ } }/* if we are allowing only "standard" bonding */ else {/* not using std bond table - */ /*so we let heavy atoms bond when they are close enough */ letItBond = TRUE; } }/* neither atom is a hydrogen - both are heavy atoms */ else if( (! isAanH || ! isBanH) /* only one atom is a hydrogen */ && ATOMS_IN_THE_SAME_RES(src, targ) ) /* and both atoms*/ /* are in the same res*/ { if (isAanH) { theHatom = src; theOtherAtom = targ; } else { theHatom = targ; theOtherAtom = src; } if( (theHatom->bondedto == NULL) /* heavy atom either unknown*/ /* or as expected */ || strstr(theHatom->bondedto, theOtherAtom->atomname) || (UseHParent == FALSE) ) { letItBond = TRUE; } } else {/*either both are hydrogens in the same residue- no bond possible*/ /* or one is a H and they are in diff res - no bond possible */ letItBond = FALSE; } if (letItBond) /*check connectivity through Maxbonded neighbors*/ { targ->mark = distcount; /*initially == 1 as called externally*/ if (distcount < max) /*max == Maxbonded as input or defaulted*/ {/*reentrant: no clashes between atoms within Maxbonded bonds*/ markBonds(targ, neighbors, distcount+1, max); } /*else targ->mark == max > discount and loop will end*/ #ifdef DUMP_DEBUG_EXTRA if (DebugLevel > 7) { fprintf(stdout, "{%4.4s%c%3.3s%s%4.4s%c}P %8.3f%8.3f%8.3f\n", src->atomname, src->altConf, src->r->resname, src->r->chain, src->r->Hy36resno, src->r->resInsCode, src->loc.x, src->loc.y, src->loc.z); fprintf(stdout, "{%4.4s%c%3.3s%s%4.4s%c}L %8.3f%8.3f%8.3f\n", targ->atomname, targ->altConf, targ->r->resname, targ->r->chain, targ->r->Hy36resno, targ->r->resInsCode, targ->loc.x, targ->loc.y, targ->loc.z); } #endif /*DUMP_DEBUG_EXTRA*/ } #endif /*def? OLD_BONDING_CALC not defined as of 050111*/ }/*close enough for a bond*/ }/*if(targ->mark == 0 || targ->mark > distcount) */ }/*for(loop over neighbors)*/ }/*markBonds*/ /*}}}markBonds()_____________________________________________________________*/ /*{{{fixupLongBondChains()****************************************************/ /* fixupLongBondChains() - unmark and remove from the bonded set */ /* atoms which are more than cutoff bonds */ /* from the source when neither src or targ */ /* is a hydrogen */ void fixupLongBondChains(atom *src, atom *neighbors, int cutoff) { atom *targ = NULL; if (! isHatom(src->elem)) { for(targ = neighbors; targ; targ = targ->scratch) { if (targ->mark > cutoff && ! isHatom(targ->elem)) { targ->mark = 0; } } } } /*}}}fixupLongBondChains()___________________________________________________*/ /*{{{dotType()****************************************************************/ int dotType(atom *src, atom *atomList, int recalcOnly) { atom *a = NULL; int rslt = src->atomclass; /* When ByNABaseColor is TRUE, the atomclass is not the real element type */ /* but instead is the nucleic acid base type. In this case, the times we */ /* need the actual element (or parent element for hydrogens) recalcOnly is*/ /* TRUE, we determine the element & do not molest the atomclass assignment*/ if (src->atomclass < 0 || recalcOnly) { rslt = src->elem; if (isHatom(src->elem)) { for(a = atomList; a; a = a->scratch) { if (a->mark == 1 && ! isHatom(a->elem) && src->r == a->r) { rslt = a->elem; break; } } } if (!recalcOnly) { /* usually reassign atomclass, but not for recalc */ src->atomclass = rslt; } } return rslt; } /*}}}dotType()_______________________________________________________________*/ /*{{{debugBondingLists()******************************************************/ void debugBondingLists(atom *src, atom *neighbors) { atom *targ = NULL; int i = 0; fprintf(stdout, "%4.4s%c%3.3s%s%4.4s%c(%c%d)", src->atomname, src->altConf, src->r->resname, src->r->chain, src->r->Hy36resno, src->r->resInsCode, src->binSerialNum, src->mark); for(targ = neighbors; targ; targ = targ->scratch) { fprintf(stdout, ", %d: %4.4s%c%3.3s%s%4.4s%c(%c%d)", ++i, targ->atomname, targ->altConf, targ->r->resname, targ->r->chain, targ->r->Hy36resno, targ->r->resInsCode, targ->binSerialNum, targ->mark); } fprintf(stdout, "\n"); } /*}}}debugBondingLists()_____________________________________________________*/ /*{{{genDotIntersect()*** called from doCommand() *** calls examineDots() ****/ /*genDotIntersect() called from doCommand during processing of modes: */ /*once with srcFlag==1, targFlag==2 in modes SELFINTERSECT,INTERSECTONCE, */ /*and twice from INTERSECTBOTHWAYS first with srcFlag,targFlag == 1,2 then 2,1*/ /*autobondrot static atoms are allMainAtoms with abin bins*/ /* mobile atoms are allMovingAtoms with bbins bins*/ /* NOTE: allMainAtoms/abins & allMovingAtoms/bbins must be disjoint */ /* sets of atoms (none in common) or allMovingAtoms/bbins can */ /* be NULL. */ /* allMovingAtoms refers to autobondrot set of atoms */ void genDotIntersect(atom *allMainAtoms, atomBins *abins, atom *allMovingAtoms, atomBins *bbins, pointSet dots[], float probeRad, float spikelen, int srcFlag, int targFlg, dotNode *results[][NODEWIDTH]) {/*genDotIntersect()*/ atom *src = NULL, *atomList = NULL, *atomList2 = NULL; int type = 0, usesMovingAtoms = FALSE; int oktargsA = TRUE, oktargsB = TRUE; int i=0,j=0;/* SJ - 05/06/2015 moved the variables here from below for compatibility with windows C compiler*/ atom *source = NULL, *target=NULL; dotNode *head = NULL, *curNode = NULL, *prevNode = NULL; if(Ldotdump) { fprintf(stderr,"entered genDotIntersect()\n"); } usesMovingAtoms = ((allMovingAtoms != NULL) && (bbins != NULL)); for(src = allMainAtoms; src; src = src->next) /*main==autobondrotstatic*/ {/*for: loop over all main atoms taking each in turn as the src-atom*/ if (src->flags & srcFlag) /*dcr?: seems srcFlag always either 1 or 2 ??*/ {/*for each src atom*/ oktargsA = TRUE; oktargsB = TRUE; atomList = findTouchingAtoms(src, NULL,abins, probeRad, targFlg, &oktargsA); /*20111206 findTouchingAtoms() sets ~exactly overlapping atoms to be NOT oktargsA*/ /*this is where too close things fail to show dots !! ?? */ if (usesMovingAtoms) /*autobondrot*/ { atomList2 = findTouchingAtoms(src, atomList, bbins, probeRad, targFlg, &oktargsB); } else { atomList2 = atomList; } if(Ldotdump) fprintf(stderr,"findTouchingAtoms() sets oktargsA== %d \n",oktargsA); if (atomList2 && (oktargsA || (usesMovingAtoms && oktargsB) )) { markBonds(src, atomList2, 1, Maxbonded); /*in genDotIntersect()*/ /*markBonds identifies bonds between atoms - */ /*seems not to have lower distance limit*/ /*but this is where autobondrot fails to spike very close atoms*/ if (Maxbonded > 3) fixupLongBondChains(src, atomList2, 3); type = dotType(src, atomList2, FALSE); if(LOneDotEach) { if(Ldotdump) fprintf(stderr,"genDotIntersect() calls examineOneDotEach()\n"); examineOneDotEach(src,type,atomList2,dots,probeRad,spikelen,targFlg,results,allMainAtoms); /*allMainAtoms20120120*/ } else { examineDots(src, type, atomList2, dots, probeRad, spikelen, targFlg, results); } if(Verbose && ShowTicks) { fprintf(stderr, "%s%d \r", src->r->resname, src->r->resid); } } }/*for each src atom*/ }/*for: loop over all static atoms taking each in turn as the src-atom*/ if(usesMovingAtoms) {/*if: usesMovingAtoms==autobondrot*/ for(src = allMovingAtoms; src; src = src->next) {/*for: loop over allMovingAtoms taking each in turn as the src-atom*/ if (src->flags & srcFlag) /*dcr?: what is being flagged here??*/ { oktargsA = TRUE; oktargsB = TRUE; atomList = findTouchingAtoms(src, NULL, abins, probeRad, targFlg, &oktargsA); atomList2 = findTouchingAtoms(src, atomList, bbins, probeRad, targFlg, &oktargsB); if (atomList2 && (oktargsA || oktargsB)) { markBonds(src, atomList2, 1, Maxbonded); /*markBonds identifies bonds between atoms - */ /*seems not to have lower distance limit*/ /*but is where autobondrot fails to spike very close atoms*/ if (Maxbonded > 3) fixupLongBondChains(src, atomList2, 3); type = dotType(src, atomList2, FALSE); /*usesMovingAtoms==autobondrot does not use LOneDotEach option 111205*/ examineDots(src,type,atomList2,dots,probeRad,spikelen,targFlg,results); if(Verbose && ShowTicks) { fprintf(stderr, "%s%d \r", src->r->resname, src->r->resid); } } } }/*for: loop over allMovingAtoms taking each in turn as the src-atom*/ }/*if: usesMovingAtoms==autobondrot*/ // SJ - 10/03/2014 - Code to delete contact dots (except H bonds) between atom pairs when atleast one of them has an altConf " " and the sum of their occpancies is <= 1 // UNKNOWN BUGS: Not sure what this code will do if the results is genearted using autobondrot for(i=0;ia; target = curNode->t; if((source->altConf == ' ' || target->altConf == ' ') && source->occ + target->occ <= 1) { //this node has to be deleted if(curNode==head) { head=curNode->next; free(curNode); prevNode=NULL; curNode=head; } else { prevNode->next=curNode->next; free(curNode); curNode=prevNode->next; } } else {//move forward when no deletion prevNode=curNode; curNode=curNode->next; } }/*while loop*/ results[i][j]=head; }/*loop over wc, cc, weak H bonds (wh, if LweakHbonds is true), so, bo, and worse overlap (wo, if LworseOverlap is true)*/ }/*for loop over atomtypes*/ //results now is filtered }/*genDotIntersect()*/ /*}}}genDotIntersect()_______________________________________________________*/ /*{{{genDotSurface() only called when method == EXTERNALSURFACE **************/ /* NOTE: allMainAtoms/abins & allMovingAtoms/bbins must be disjoint */ /* sets of atoms (none in common) or allMovingAtoms/bbins can */ /* be NULL. */ /* allMovingAtoms refers to autobondrot set of atoms */ void genDotSurface(atom *allMainAtoms, atomBins *abins, atom *allMovingAtoms, atomBins *bbins, pointSet dots[], float probeRad, float spikelen, int srcFlag, dotNode *results[][NODEWIDTH]) { atom *src = NULL, *atomList = NULL, *atomList2 = NULL; int type = 0, usesMovingAtoms = FALSE; int dummy = TRUE; usesMovingAtoms = ((allMovingAtoms != NULL) && (bbins != NULL)); for(src = allMainAtoms; src; src = src->next) { if (src->flags & srcFlag) { atomList = findTouchingAtoms(src, NULL, abins, probeRad, 0, &dummy); if (usesMovingAtoms) { atomList2 = findTouchingAtoms(src, atomList, bbins, probeRad, 0, &dummy); } else { atomList2 = atomList; } if (atomList2) { markBonds(src, atomList2, 1, 1); /*in genDotSurface()*/ } type = dotType(src, atomList2, FALSE); surfDots(src, type, atomList2, dots, probeRad, spikelen, results); if(Verbose && ShowTicks) { fprintf(stderr, "%s%d \r", src->r->resname, src->r->resid); } } } if (usesMovingAtoms) { for(src = allMovingAtoms; src; src = src->next) { if (src->flags & srcFlag) { atomList = findTouchingAtoms(src, NULL, abins, probeRad, 0, &dummy); atomList2 = findTouchingAtoms(src, atomList, bbins, probeRad, 0, &dummy); if (atomList2) { markBonds(src, atomList2, 1, 1); /*in genDotSurface()*/ } type = dotType(src, atomList2, FALSE); surfDots(src, type, atomList2, dots, probeRad, spikelen, results); if(Verbose && ShowTicks) { fprintf(stderr, "%s%d \r", src->r->resname, src->r->resid); } } } } } /*}}}genDotIntersect()_______________________________________________________*/ /*{{{surfDots() only called from genDotSurface(), method == EXTERNALSURFACE **/ void surfDots(atom *src, int type, atom *scratch, pointSet dots[], float probeRad, float spikelen, dotNode *results[][NODEWIDTH]) { atom *targ; int i, nearpt, ok; point3d dotloc, dotvect, exploc, spikeloc; pointSet *srcDots; char ptmaster = ' '; /*pointmaster character dcr041009*/ if (src->elem == ignoreAtom) { return; } srcDots = &(dots[src->elem]); if (src->elem == atomC && isCarbonylAtom(src)) { srcDots = &COdots; #ifdef DUMP_DEBUG_EXTRA if (DebugLevel>8) { fprintf(stdout, "DEBUG surfDots(%s%d %s radius = %.3f)\n", src->r->resname,src->r->resid,src->atomname,src->radius); } #endif } for(i=0; i < srcDots->n; i++) { dotvect = srcDots->p[i]; v3add(&(src->loc), &dotvect, &dotloc); v3scale(&dotvect, src->radius + probeRad); v3add(&(src->loc), &dotvect, &exploc); v3scale(&dotvect, src->radius + 0.0); v3add(&(src->loc), &dotvect, &spikeloc); ok = TRUE; /*targ is an atom*, scratch is an atom*, atom* defined in abin.h */ /* atom* scratch is a member of atom*: which has member scratch...*/ /* so for-loop moves through all these atoms while an atom has a scratch*/ for(targ = scratch; targ; targ = targ->scratch) { if ( (targ->elem == ignoreAtom) || ((!DoWatWat) && (src->props & WATER_PROP) && (targ->props & WATER_PROP) && (src->r != targ->r)) /* we do see HOH hydrogens */ || ((!DoHet) && (targ->props & HET_PROP)) || ((!DoH2O) && (targ->props & WATER_PROP)) ) { /* ignore */ } else { #ifdef INLINE_FOR_SPEED nearpt = INRANGEINLINE(exploc, (targ->loc), probeRad+targ->radius); #else nearpt = inRange(&exploc, &(targ->loc), probeRad+targ->radius); #endif if (nearpt) { ok = FALSE; break; } } } if (ok) { saveDot(src, NULL, type, &dotloc, &spikeloc, results, 0, 0.0,ptmaster,0); /*XHTangle dummy dcr20120120*/ /*ptmaster not used to distinguish anykind of surface dots dcr041009*/ /*Note: here only when method == EXTERNALSURFACE */ /*and ovrlaptype==0, mingap==0.0, ptmaster==' ' */ } } }/*surfDots()*/ /*}}}genDotIntersect()_______________________________________________________*/ /*{{{initResults()************************************************************/ void initResults(dotNode *results[][NODEWIDTH]) { int i, j; for (i = 0; i < NUMATOMTYPES; i++) { for (j = 0; j < NODEWIDTH; j++) { results[i][j] = NULL; } } } /*}}}initResults()___________________________________________________________*/ /*{{{freeResults()************************************************************/ void freeResults(dotNode *results[][NODEWIDTH]) { int i, j; dotNode *node, *next; for (i = 0; i < NUMATOMTYPES; i++) { for (j = 0; j < NODEWIDTH; j++) { for (node = results[i][j]; node; node = next) { next = node->next; free(node); } results[i][j] = NULL; } } } /*}}}freeResults()___________________________________________________________*/ /*{{{assignGapColorForKin()***************************************************/ char* assignGapColorForKin(float gap, int class) { char *colorValue = ""; //if (class == 4) { colorValue = "greentint "; } /* hbond */ /* 04/16/2015 SJ changed 4 to 6, as H-bond are now 6 in NODEWIDTH, see probe.h*/ if (class == 6) { colorValue = "greentint "; } else if (gap > 0.35){ colorValue = "blue "; } else if (gap > 0.25){ colorValue = "sky "; } else if (gap > 0.15){ colorValue = "sea "; } else if (gap > 0.0) { colorValue = "green "; } else if (gap >-0.1) { colorValue = "yellowtint ";} else if (gap >-0.2) { colorValue = "yellow "; } else if (gap >-0.3) { colorValue = "orange "; } else if (gap >-0.4) { colorValue = "red "; } else { colorValue = "hotpink "; } return colorValue; } /*}}}assignGapColorForKin()__________________________________________________*/ /*{{{assignGapColorForO()*****************************************************/ char* assignGapColorForO(float gap, int class) { char *colorValue = ""; //if (class == 4) { colorValue = "pale_green "; } /* hbond *//* 04/16/2015 SJ changed 4 to 6, as H-bond are now 6 in NODEWIDTH, see probe.h*/ if (class == 6) { colorValue = "pale_green "; } else if (gap > 0.35){ colorValue = "cornflower_blue "; } else if (gap > 0.25){ colorValue = "sky_blue "; } else if (gap > 0.15){ colorValue = "aquamarine "; } else if (gap > 0.0) { colorValue = "green "; } else if (gap >-0.1) { colorValue = "yellow_green "; } else if (gap >-0.2) { colorValue = "yellow "; } else if (gap >-0.3) { colorValue = "orange "; } else if (gap >-0.4) { colorValue = "red "; } else { colorValue = "orange_red "; } return colorValue; } /*}}}assignGapColorForO()____________________________________________________*/ /*{{{assignGapColorForXV()****************************************************/ char* assignGapColorForXV(float gap, int class) { return assignGapColorForKin(gap, class); } /*}}}assignGapColorForXV()___________________________________________________*/ /*{{{convertKinColorToO()*****************************************************/ char* convertKinColorToO(char* incolor) { char *outcolor = "light_gray"; if (strcmp(incolor, "red" ) == 0) { outcolor = "red"; } else if (strcmp(incolor, "green" ) == 0) { outcolor = "green"; } else if (strcmp(incolor, "blue" ) == 0) { outcolor = "cornflower_blue"; } else if (strcmp(incolor, "cyan" ) == 0) { outcolor = "cyan"; } else if (strcmp(incolor, "yellow" ) == 0) { outcolor = "yellow"; } else if (strcmp(incolor, "magenta" ) == 0) { outcolor = "magenta"; } else if (strcmp(incolor, "white" ) == 0) { outcolor = "white"; } else if (strcmp(incolor, "pink" ) == 0) { outcolor = "salmon"; } else if (strcmp(incolor, "orange" ) == 0) { outcolor = "orange"; } else if (strcmp(incolor, "purple" ) == 0) { outcolor = "purple"; } else if (strcmp(incolor, "sky" ) == 0) { outcolor = "sky_blue"; } else if (strcmp(incolor, "brown" ) == 0) { outcolor = "brown"; } else if (strcmp(incolor, "gray" ) == 0) { outcolor = "light_gray"; } else if (strcmp(incolor, "black" ) == 0) { outcolor = "black"; } else if (strcmp(incolor, "gold" ) == 0) { outcolor = "gold"; } else if (strcmp(incolor, "yellowtint") == 0) { outcolor = "yellow_green"; } else if (strcmp(incolor, "sea" ) == 0) { outcolor = "aquamarine"; } else if (strcmp(incolor, "pinktint" ) == 0) { outcolor = "pink"; } else if (strcmp(incolor, "bluetint" ) == 0) { outcolor = "light_blue"; } else if (strcmp(incolor, "greentint" ) == 0) { outcolor = "pale_green"; } else if (strcmp(incolor, "hotpink" ) == 0) { outcolor = "orange_red"; } else if (strcmp(incolor, "invisible" ) == 0) { outcolor = "black"; } else { outcolor = "light_gray"; } return outcolor; } /*}}}convertKinColorToO()____________________________________________________*/ /*{{{convertKinColorToXV()****************************************************/ char* convertKinColorToXV(char* incolor) { return incolor; } /*}}}convertKinColorToXV()___________________________________________________*/ /*{{{writeOutput()************************************************************/ void writeOutput(FILE *outf, char* groupname, dotNode *results[][NODEWIDTH], int spike, int method, char* extrastr, float probeRad) { /*writeOutput for kinemage*/ /*add probeRad 20111220dcr*/ /*dcr041020 need method for better kinemage keywords*/ /* 060129 extrastr for extra master to control orig vs fitted dots */ int i, j; int k; /*counter for gapbins 111020dcr*/ char ptm[2] = {' ','\0'}; /*20120120dcr ptmaster for bval*/ dotNode *node; atom *a; /*owner of dot*/ atom *t; /*cause of dot, output for OneDotEach mode 110113*/ char *color = ""; // OLD: char *mast[NODEWIDTH] = {"wide contact", "close contact", "small overlap", // "bad overlap", "H-bonds", "H-weaks"}; /*111215dcr weakHbonds option*/ char *mast[NODEWIDTH] = {"wide contact", "close contact", "weak H-bonds", "small overlap", "bad overlap", "worse overlap", "H-bonds"}; /* 04/13/2015 SJ changed to match new NODEWIDTH, see probe.h*/ char *surfacestr = "surface"; /*041020*/ char *contactstr = "vdw contact"; /*060129*/ char extraMstr[32]; /*fill below with extra master name 060129*/ char pointid[100], lastpointid[100]; char ptmast[7]={'\0','\0','\0','\0','\0','\0','\0'}; /*111020dcr 6 to 7 */ /*dcr041009 string to hold 'M' (McMc), 'S' (ScSc), 'P' (McSc), 'O' (other)*/ /*dcr041009 ptmaster part of each dotNode member of results[][]*/ char masterchr=' '; /*dcr041017*/ /*111020dcr gapbin ptmasters*/ int ngapbins = 22; /*(std probe diameter .5 +clash overlap .5=1.0/0.05interval) +2*/ /*+1 for 0th fence post, another +1 for X excess*/ char gapnamestr[22] = {'z','y','x','w','v','u','t','g','r','q','f','F','Q','R','G','T','U','V','W','X','Y','Z'}; /*0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 ,10 ,11 ,12 ,13 ,14 ,15 ,16 ,17 ,18 ,19 ,20 ,21*/ /*std gapbins scope at least -.5 to +.5, wider if probeRad > 0.25 standard*/ int gaplimit = ((2*(probeRad>0.25f?probeRad:0.25f) +0.5f)/0.05f)+2; /*20111220dcr*/ long *gapcounts = (long *) malloc(gaplimit * sizeof(long)); long maxgapcounts = 0; int Lgotgapbin = FALSE; char strcName[6]; /*first nendstr char of inputfilename 20120120dcr*/ int nendstr = 6; if(Lgapbins) /*111020dcr gapbins pointmasters, 111215dcr count bin entries*/ { for(k=0; k 0) if(mcmccont[0][7] > 0) { masterchr = MCMCCHR; /*single char*/ fprintf(outf, "@pointmaster '%c' {McMc contacts}\n",masterchr); } //if(scsccont[0][5] > 0) if(scsccont[0][7] > 0) { masterchr = SCSCCHR; /*single char*/ fprintf(outf, "@pointmaster '%c' {ScSc contacts}\n",masterchr); } //if(mcsccont[0][5] > 0) if(mcsccont[0][7] > 0) { masterchr = MCSCCHR; /*single char*/ fprintf(outf, "@pointmaster '%c' {McSc contacts}\n",masterchr); } //if(othrcont[0][5] > 0) if(othrcont[0][7] > 0) { masterchr = OTHERCHR; /*single char*/ fprintf(outf, "@pointmaster '%c' {Hets contacts}\n",masterchr); } if(Lgapbins) /*111020dcr*/ { for(i = 0; i < ngapbins; i++) { fprintf(outf, "@pointmaster '%c' {gap %3.2f}\n" ,gapnamestr[i],(((float)i-11.0)/20.0)+0.05); /*20111220dcr +0.05*/ } } for (i = 0; i < NUMATOMTYPES; i++) {/*i: Hatom varients + rest of periodic table + base types, see atomprops.h*/ for (j = 0; j < NODEWIDTH; j++) {/* OLD j: wide contact, close contact, small overlap, bad overlap, H-bonds*/ /*OLD optionally: H-weaks if LweakHbonds 111215dcr*/ /* NEW 04/13/2015 - SJ j: wide contact, close contact, weak H bonds, small overlap, bad overlap, worse overlap, H-bonds*/ node = results[i][j]; /*head of list of dot nodes, i.e. first node*/ if (node) {/*node exists, write list headers*/ if (j == 0 || j == 1) /*wide and close contacts*/ {/* contact */ if (AtomMasters) { if(LOneDotEach) { fprintf(outf, "@balllist {x} radius=0.2 color=%s master={%s dots} master={%s}%s%s\n", getColor(i), getAtomName(i), mast[j],extraMstr, LensDots ? " lens" : ""); } else { fprintf(outf, "@dotlist {x} color=%s master={%s dots} master={%s}%s%s\n", getColor(i), getAtomName(i), mast[j],extraMstr, LensDots ? " lens" : ""); } } else { if(LOneDotEach) { fprintf(outf, "@balllist {x} radius=0.2 color=%s master={%s}%s%s\n", getColor(i), mast[j],extraMstr, LensDots ? " lens" : ""); } else { fprintf(outf, "@dotlist {x} color=%s master={%s}%s%s\n", getColor(i), mast[j],extraMstr, LensDots ? " lens" : ""); } } }/* contact */ // else if ((j == 2 || j == 3) && spike) else if ((j == 3 || j == 4 || j == 5) && spike) /* 04/13/2015 SJ changed to reflect new numbers for different categories*/ {/* bump w/ spike */ if (AtomMasters) { if(LOneDotEach) { fprintf(outf, "@vectorlist {x} width=7 color=%s master={%s dots} master={%s}%s\n", getColor(i), getAtomName(i), mast[j],extraMstr); } else { fprintf(outf, "@vectorlist {x} color=%s master={%s dots} master={%s}%s\n", getColor(i), getAtomName(i), mast[j],extraMstr); } } else { if(LOneDotEach) { fprintf(outf, "@vectorlist {x} width=7 color=%s master={%s}%s\n", getColor(i), mast[j],extraMstr); } else { fprintf(outf, "@vectorlist {x} color=%s master={%s}%s\n", getColor(i), mast[j],extraMstr); } } }/* bump w/ spike */ //else if (j == 4) else if (j == 6)/* 04/13/2015 SJ changed to reflect new numbers for different categories*/ {/* bump w/o spike or H-bond */ if (AtomMasters) { if(LOneDotEach) { fprintf(outf, "@ringlist {x} radius=0.2 color=%s master={%s dots} master={%s}%s\n", getColor(i), getAtomName(i), mast[j],extraMstr); } else { fprintf(outf, "@dotlist {x} color=%s master={%s dots} master={%s}%s\n", getColor(i), getAtomName(i), mast[j],extraMstr); } } else { if(LOneDotEach) { fprintf(outf, "@ringlist {x} radius=0.2 color=%s master={%s}%s\n", getColor(i), mast[j],extraMstr); } else { fprintf(outf, "@dotlist {x} color=%s master={%s}%s\n", getColor(i), mast[j],extraMstr); } } }/* bump w/o spike or H-bond */ else if (LweakHbonds) {/* H-weak */ if (AtomMasters) { if(LOneDotEach) { fprintf(outf, "@ringlist {x} radius=0.2 color=%s master={%s dots} master={%s}%s\n", getColor(i), getAtomName(i), mast[j],extraMstr); } else { fprintf(outf, "@dotlist {x} color=%s master={%s dots} master={%s}%s\n", getColor(i), getAtomName(i), mast[j],extraMstr); } } else { if(LOneDotEach) { fprintf(outf, "@ringlist {x} radius=0.2 color=%s master={%s}%s\n", getColor(i), mast[j],extraMstr); } else { fprintf(outf, "@dotlist {x} color=%s master={%s}%s\n", getColor(i), mast[j],extraMstr); } } }/* H-weak */ lastpointid[0] = '\0'; /* reset */ }/*node exists, write list headers*/ while(node) {/*kinemage point or point-line for each node in (ij)th dot-node list*/ a = node->a; /*owner of dot*/ t = node->t; /*cause of dot, use for OneDotEach 111013dcr*/ if(Lgapbins) /*111020dcr*/ {/*two characters for two pointmasters*/ ptmast[0] = ' '; /*insure leading space*/ ptmast[1] = '\''; /*surround char pointmasters with*/ ptmast[2] = node->ptmaster; ptmast[3] = ' '; /*space for later gapbin character 111020dcr*/ ptmast[4] = '\''; /*single quote marks*/ ptmast[5] = ' '; /*need trailing space*/ ptmast[6] = '\0'; /*end string*/ } else if(node->ptmaster == ' ') /*dcr041009*/ {/*blank means NO pointmaster*/ ptmast[0] = '\0'; /*which will print as a zero-length str*/ } else { ptmast[0] = ' '; /*insure leading space*/ ptmast[1] = '\''; /*surround single char pointmaster with*/ ptmast[2] = node->ptmaster; ptmast[3] = '\''; /*single quote marks*/ ptmast[4] = ' '; /*need trailing space*/ ptmast[5] = '\0'; /*end string*/ } /* sprintf(pointid, "%s%c %s%d%c", */ /* a->atomname, a->altConf, */ /* a->r->resname, a->r->resid, */ /* a->r->resInsCode); */ if(LOneDotEach) /*report cause of dot after owner of dot 111013dcr*/ { sprintf(pointid, "%s%c%s%s%c%s %s%c%s%s%c%s GAP %6.3f", a->atomname, a->altConf, a->r->resname, a->r->Hy36resno, a->r->resInsCode, a->r->chain, t->atomname, t->altConf, t->r->resname, t->r->Hy36resno, t->r->resInsCode, t->r->chain, node->gap); } else {/*usual just report owner of dot*/ sprintf(pointid, "%s%c%s%s%c%s", a->atomname, a->altConf, a->r->resname, a->r->Hy36resno, a->r->resInsCode, a->r->chain); } if (strcmp(pointid, lastpointid)) { strcpy(lastpointid, pointid); fprintf(outf, "{%s}", pointid); } else { fprintf(outf, "{\"}"); } if (ColorGap) { if (node->t) { color = assignGapColorForKin(node->gap, j); fprintf(outf, "%s", color); } else { fprintf(outf, "%s ", OutPointColor); } /* added "%s" string format to color and */ /* edited "%s " to "%s" in OutPointColor, wba 110909 */ /* edited "%s" back to "%s " in OutPointColor, gjk 1200907 */ } if(Lgapbins) /*111020dcr gapbins pointmasters*/ { Lgotgapbin = FALSE; /*until identify which gapbin*/ for(k=0; kgap < (((float)k-11.0)/20.0)+0.05 ) {/*pt master intervals of 0.05 from -0.5 to +0.5*/ ptmast[3] = gapnamestr[k]; gapcounts[k] = gapcounts[k] +1; if(gapcounts[k] > maxgapcounts){maxgapcounts = gapcounts[k];} if(k < ngapbins) { Lgotgapbin = TRUE; /*got the gapbin idenity for this dot*/ } break; /*out of for loop*/ } } if(!Lgotgapbin) {/*assign this node, aka dot, to overflow gapbin*/ ptmast[3] = gapnamestr[ngapbins-1]; /*the overflow bin-master*/ gapcounts[gaplimit-1] = gapcounts[gaplimit-1] +1; /*last gapbin*/ /*if(gapcounts[gaplimit-1] > maxgapcounts) {maxgapcounts = gapcounts[gaplimit-1];}*/ } } //if ((j == 2 || j == 3) && spike) if ((j == 3 || j == 4 || j == 5) && spike)/* 04/13/2015 SJ changed to reflect new numbers for different categories*/ {/* bump */ fprintf(outf, "P %s%.3f,%.3f,%.3f {\"}%s %s%.3f,%.3f,%.3f\n", /*dcr041009*/ ptmast,node->loc.x, node->loc.y, node->loc.z, /*dcr041009*/ color, /*** note: second color reference */ ptmast,node->spike.x, node->spike.y, node->spike.z); /*dcr041009*/ } else {/* contact or H-bond */ fprintf(outf, "%s%.3f,%.3f,%.3f\n", /*dcr041009*/ ptmast,node->loc.x, node->loc.y, node->loc.z); /*dcr041009*/ } node = node->next; /*in this (ij)th dot-node list*/ }/*kinemage point or point-line for each node in (ij)th dot-node list*/ }/*OLD j: wide contact, close contact, small overlap, bad overlap, H-bonds*/ /* see new def in NODEWIDTH probe.h*/ }/*i: Hatom varients + rest of periodic table + base types, see atomprops.h*/ if(Lgapbins) /*111020dcr gapbins pointmasters, 111215dcr count bin entries*/ {/*Lgapbins TRUE*/ fprintf(outf, "@text\n"); for(k=0; k50}\n"); fprintf(outf, "@group {AXES} dominant\n"); fprintf(outf, "@labellist {gapbins}\n"); for(k=0; k<(gaplimit-1);k++) { fprintf(outf, "{%5.2f} %5.2f, 0.00, 0.00\n", ((((float)k-11.0)/20.0)+0.05), 100*(((((float)k-11.0)/20.0)+0.05)) ); } for(k=90; k<=180; k=k+10) { fprintf(outf, "{%5.2f} -70.0, %5.2f, 0.00\n",(float)k,(float)(k-90)); } if(sizeof(inputfilename) < nendstr) {nendstr = sizeof(inputfilename);} for(k=0; ka; /*owner of dot*/ t = node->t; /*cause of dot, use for OneDotEach 111013dcr*/ /* Slow way to determine integer truncation of k/5 with a maximum of 9. */ for(k=0; k<10; k++){if(a->bval < k*5) {break;}} /*2:10A,3:15B,4:20C,5:25D,6:30E,8:40F,10:50G*/ if(k<=2){ptm[0] = 'A';} else if(k<=3){ptm[0] = 'B';} else if(k<=4){ptm[0] = 'C';} else if(k<=5){ptm[0] = 'D';} else if(k<=6){ptm[0] = 'E';} else if(k<=7){ptm[0] = 'F';} else if(k<=8){ptm[0] = 'G';} else{ptm[0] = 'H';} /*Note first nemdstr char of inputfilename 20120120dcr*/ fprintf(outf,"{%6s %s%c%s%s%c%s:O=%4.2f, B=%6.2f; %s%c%s%s%c%s:O=%4.2f, B=%6.2f; GAP %6.3f}'%s' %8.3f, %8.3f , %8.3f\n", strcName, a->atomname, a->altConf, a->r->resname, a->r->Hy36resno, a->r->resInsCode, a->r->chain, a->occ, a->bval, t->atomname, t->altConf, t->r->resname, t->r->Hy36resno, t->r->resInsCode, t->r->chain, t->occ, t->bval, node->gap, ptm, 100*(node->gap), (node->angle)-90.0, a->bval); node = node->next; /*in this (ij)th dot-node list*/ }/*kinemage point or point-line for each node in (ij)th dot-node list*/ fprintf(outf, "\n"); }/*node exists*/ }/*NODEWIDTH*/ }/*NUMATOMTYPES*/ }/*LXHvector*/ }/*Lgapbins TRUE*/ }/*writeOutput for kinemage*/ /*}}}writeOutput()___________________________________________________________*/ /*{{{writeAltFmtO()***********************************************************/ void writeAltFmtO(FILE *outf, int showBegin, int showEnd, char* groupname, dotNode *results[][NODEWIDTH], int spike) { int i, j, numGroups, gn; dotNode *node; char *color = "", *prevColor = ""; //Old: char *mast[NODEWIDTH] = {"WC", "CC", "SO", "BO", "HB"}; char *mast[NODEWIDTH] = {"WC", "CC", "WH", "SO", "BO", "WO", "HB"}; /* 04/15/2015 SJ changed to match new NODEWIDTH, see probe.h*/ numGroups = (showBegin && showEnd) ? 1 : 2; gn = (!showBegin && showEnd) ? 2 : 1; for (j = 0; j < NODEWIDTH; j++) { fprintf(outf, "begin_object %s%d\n", mast[j], gn); fprintf(outf, "mode solid\n"); for (i = 0; i < NUMATOMTYPES; i++) { prevColor = "--none--"; node = results[i][j]; if (node) { color = convertKinColorToXV(getColor(i)); } while(node) { if (ColorGap) { if (node->t) { color = assignGapColorForO(node->gap, j); } else { color = "124";} } if (color != prevColor) { fprintf(outf, "colour %s\n", color); prevColor = color; } if ((j == 3 || j == 4 || j ==5) && spike) {/* bump */ /* 04/15/2015 SJ changed j numbers to match new ones for so, bo, and wo*/ fprintf(outf, "move %.3f %.3f %.3f\nline %.3f %.3f %.3f\n", node->loc.x, node->loc.y, node->loc.z, node->spike.x, node->spike.y, node->spike.z); } else {/* contact or H-bond */ fprintf(outf, "sphere_xyz %.3f %.3f %.3f 0.03\n", node->loc.x, node->loc.y, node->loc.z); } node = node->next; } } fprintf(outf, "end_object\n"); } if (showEnd) { if (showBegin) { fprintf(outf, "begin_object %s\n", groupname); } else { fprintf(outf, "begin_object contsurf\n"); } for (j = 0; j < NODEWIDTH; j++) { for (i = 0; i < numGroups; i++) { fprintf(outf, "instance %s%d\n", mast[j], i+1); } } fprintf(outf, "end_object\n"); } } /*}}}writeAltFmtO()__________________________________________________________*/ /*{{{writeAltFmtXV()**********************************************************/ void writeAltFmtXV(FILE *outf, int showBegin, int showEnd, char* groupname, dotNode *results[][NODEWIDTH], int spike) { int i, j; dotNode *node; char *color = ""; //Old: char *mast[NODEWIDTH] = {"wide contact", "close contact", "small overlap", "bad overlap", "H-bonds"}; /* 04/15/2015 SJ changed to reflect new NODEWIDTH, see probe.h*/ char *mast[NODEWIDTH] = {"wide contact", "close contact", "weak H-bonds", "small overlap", "bad overlap", "worse overlap", "H-bonds"}; if (showBegin) { if (showEnd) { fprintf(outf, "# begin %s\n", groupname); } else { fprintf(outf, "# begin object\n# group: %s\n", groupname); } } else { fprintf(outf, "# group: %s\n", groupname); } for (j = 0; j < NODEWIDTH; j++) { for (i = 0; i < NUMATOMTYPES; i++) { node = results[i][j]; if (node) { fprintf(outf, "# (%s %s)\n", getAtomName(i), mast[j]); color = convertKinColorToXV(getColor(i)); } while(node) { if (ColorGap) { if (node->t) { color = assignGapColorForXV(node->gap, j); } else { color = "124";} } if ((j == 3 || j == 4 || j == 5) && spike) {/* bump */ /* 04/15/2015 SJ changed j numbers to match new ones for so, bo, and wo*/ fprintf(outf, "%.3f %.3f %.3f %.3f %.3f %.3f %s\n", node->loc.x, node->loc.y, node->loc.z, node->spike.x, node->spike.y, node->spike.z, color); } else {/* contact or H-bond */ fprintf(outf, "%.3f %.3f %.3f %.3f %.3f %.3f %s\n", node->loc.x, node->loc.y, node->loc.z, node->loc.x, node->loc.y, node->loc.z, color); } node = node->next; } } } if (showEnd) { if (showBegin) { fprintf(outf, "# end %s\n", groupname); } else { fprintf(outf, "# endgroup: %s\n# end object\n", groupname); } } else { fprintf(outf, "# endgroup: %s\n", groupname); } } /*}}}writeAltFmtXV()_________________________________________________________*/ /*{{{dot2bullsEye()***********************************************************/ float dot2bullsEye(point3d *dot, atom *src, atom *targ) { point3d targ2srcVec, targSurfacePoint; v3sub(&(src->loc), &(targ->loc), &targ2srcVec); v3scale(&targ2srcVec, targ->radius); v3add(&targ2srcVec, &(targ->loc), &targSurfacePoint); return v3distance(dot, &targSurfacePoint); } /*}}}dot2bullsEye()__________________________________________________________*/ /*{{{dot2srcCenter()**********************************************************/ float dot2srcCenter(point3d *dot, atom *src, atom *targ) { point3d src2targVec, srcSurfacePoint; v3sub(&(targ->loc), &(src->loc), &src2targVec); v3scale(&src2targVec, src->radius); v3add(&src2targVec, &(src->loc), &srcSurfacePoint); return v3distance(dot, &srcSurfacePoint); } /*}}}dot2srcCenter()_________________________________________________________*/ /*{{{kissEdge2bullsEye()******************************************************/ float kissEdge2bullsEye(float ra, float rb, float rp) { return 2.0*ra*sqrt(rb*rp/((ra+rb)*(ra+rp))); } /*}}}kissEdge2bullsEye()_____________________________________________________*/ /*{{{writeRaw()***************************************************************/ void writeRaw(FILE *outf, char* groupname, dotNode *results[][NODEWIDTH], float probeRad, char* label, float density, int conFlag) //conFlag added SJ 10/07/2011 { int i, j; dotNode *node; atom *a, *t; // char *mast[NODEWIDTH] = {"wc", "cc", "so", "bo", "hb"}; char *mast[NODEWIDTH] = {"wc","cc","wh","so","bo","wo","hb"}; /* 04/10/2015 SJ to changed to match new categories of weak H bonds and worse overlap (these lists are made only if LweakHbonds and LworseOverlap are true respectively). see NODEWIDTH probe.h*/ float gap, sl, dtgp, ke2be, d2be, d2sc, score; double scaledGap; for (i = 0; i < NUMATOMTYPES; i++) {/*loop over NUMATOMTYPES*/ for (j = 0; j < NODEWIDTH; j++) {/*loop over NODEWIDTH*/ node = results[i][j]; if(node) // added 10/05/11 - SJ for condensing rawOutput { node = Condense(node,conFlag); results[i][j]=node; } while(node) {/*while node...*/ fprintf(outf, "%s:%s:%s:", label, groupname, mast[j]); a = node->a; t = node->t; fprintf(outf, "%s%4.4s%c%s %s%c:", a->r->chain, a->r->Hy36resno, a->r->resInsCode, a->r->resname,a->atomname, a->altConf); if(t) {/*t node exists*/ fprintf(outf, "%s%4.4s%c%s %s%c:", t->r->chain, t->r->Hy36resno, t->r->resInsCode, t->r->resname,t->atomname, t->altConf); gap = gapSize(&(a->loc), &(t->loc),(a->radius + t->radius)); dtgp = node->gap; if (OldStyleU) { ke2be = kissEdge2bullsEye(a->radius, t->radius, probeRad); d2be = dot2bullsEye(&(node->loc), a, t); d2sc = dot2srcCenter(&(node->loc), a, t); } sl = gapSize(&(node->loc), &(node->spike), 0.0); score = 0.0; switch(j) /* SJ 04/10/2015 changed to include new categories*/ { case 0: case 1: scaledGap = dtgp/GAPweight; score = exp(-scaledGap*scaledGap); break; case 2: /* TODO: SJ weak H bonds, don't know what to do here, because they can be both wc and cc, so will have to check*/ case 3: /* small overlap, doing nothing, as before*/ case 4: score = -BUMPweight * sl; break; case 5: score = -BUMPweight * sl; break; /* worse overlap, same as bad overlap*/ case 6: score = HBweight * sl; break; } if(conFlag) { fprintf(outf,"%i:",node->dotCount); } else if(LOneDotEach) /*dcr20120120*/ { fprintf(outf,"%.2f:",node->angle); } if (OldStyleU) { fprintf(outf, "%.3f:%.3f:%.3f:%.3f:%.3f:%.3f:%.4f", gap, dtgp, ke2be, d2be, d2sc, sl, score/density); } else { fprintf(outf, "%.3f:%.3f:%.3f:%.3f:%.3f:%.3f:%.4f", gap, dtgp, node->spike.x, node->spike.y, node->spike.z, sl, score/density); } }/*t node exists*/ else { fprintf(outf, ":::::::"); } fprintf(outf,":%s:%s:%.3f:%.3f:%.3f", getAtomName(i),(t?getAtomName(t->atomclass):""), node->loc.x,node->loc.y,node->loc.z); if(t) { fprintf(outf, ":%.2f:%.2f\n", a->bval, t->bval); } else { fprintf(outf, ":%.2f:\n", a->bval); } node = node->next; }/*while node...*/ }/*loop over NODEWIDTH*/ }/*loop over NUMATOMTYPES*/ } /*}}}writeRaw()______________________________________________________________*/ /*{{{Condense()****(SortdotNodes)***** - SJ 10/05/11 for condensing rawOutput*/ dotNode * Condense(dotNode * head, int conFlag) { // variables of the main for loop char sourceCurStart[20],sourceCurEnd[20]; dotNode *curStart = NULL, *curEnd = NULL, *prevEnd = NULL, *nextStart = NULL, *prev = NULL; atom *source = NULL; //variables for the sorting and removing duplciates loop dotNode *i = NULL, *j = NULL, *next = NULL, *iPrev = NULL, *jPrev = NULL, *minPrev = NULL, *minNode = NULL; char minTarget[20],jTarget[20],prevTarget[20],nextTarget[20]; atom *target = NULL; for(curStart=head;curStart!=NULL;curStart=nextStart) { /* main for loop which loops through the source atoms*/ source=curStart->a; sprintf(sourceCurStart,"%s%4.4s%c%s %s%c",source->r->chain,source->r->Hy36resno, source->r->resInsCode,source->r->resname,source->atomname,source->altConf); //sourceCurStart[17]='\0'; nextStart = NULL; prev = curStart; curEnd = NULL; for(curEnd=curStart->next;curEnd!=NULL;curEnd=curEnd->next) { /* for loop to determine the nodes for the same source atom*/ source=curEnd->a; sprintf(sourceCurEnd,"%s%4.4s%c%s %s%c",source->r->chain,source->r->Hy36resno, source->r->resInsCode,source->r->resname,source->atomname,source->altConf); //sourceCurEnd[17]='\0'; if(strcmp(sourceCurStart,sourceCurEnd)!=0) // if the source atom differes, then break out of the loop break; prev=curEnd; } nextStart=curEnd; // the start of the next source atoms curEnd=prev; // the last dot for the current source atom // isolating all the dots of the current source atom from the prev and next source atom if(prevEnd) prevEnd->next = NULL; curEnd->next = NULL; /* sorting code - will have a curState and a curEnd pointed to indicate the start and end of the sorted list*/ iPrev = NULL; jPrev = NULL; minPrev = NULL; i = NULL; j = NULL; next = NULL; minNode = NULL; for(i=curStart;i!=NULL;i=i->next) { /*sorting loop - sorting based on the target atom*/ target=i->t; sprintf(minTarget,"%s%4.4s%c%s %s%c",target->r->chain,target->r->Hy36resno, target->r->resInsCode,target->r->resname,target->atomname,target->altConf); //minTarget[17]='\0'; minNode = i; minPrev = iPrev; //using selection sort jPrev = i; for(j=i->next;j!=NULL; j=j->next) { target=j->t; sprintf(jTarget,"%s%4.4s%c%s %s%c",target->r->chain,target->r->Hy36resno, target->r->resInsCode,target->r->resname,target->atomname,target->altConf); //jTarget[17]='\0'; // if the target atom of j is less than that of the min then swap it: // also swap if the names are the same but the gap is less so that we // get the list sorted from smallest (most negative) to largest gap. // This will be important when we are using condensed output so that we // pick the node with the smallest gap. if ( (strcmp(jTarget,minTarget) < 0) || ((strcmp(jTarget, minTarget) == 0) && (j->gap < minNode->gap) ) ) { strcpy(minTarget,jTarget); minNode = j; minPrev = jPrev; } jPrev = j; } if(minNode != i) // if the nodes have to be swapped { next = minNode->next; if(iPrev) iPrev->next = minPrev -> next; else curStart=minNode; minPrev->next = i; minNode->next = i->next; i->next = next; i = minNode; } iPrev = i; }/* sorting loop ends */ curEnd = iPrev; //updating the end of the dots for the source atom if(conFlag) {/* code for removing duplicates - will update curEnd */ i=curStart; next = NULL; i->dotCount=1; //initialising the dot count while(i->next!=NULL) {//removing duplicates target=i->t; sprintf(prevTarget,"%s%4.4s%c%s %s%c",target->r->chain,target->r->Hy36resno, target->r->resInsCode,target->r->resname,target->atomname,target->altConf); //prevTarget[17]='\0'; target=i->next->t; sprintf(nextTarget,"%s%4.4s%c%s %s%c",target->r->chain,target->r->Hy36resno, target->r->resInsCode,target->r->resname,target->atomname,target->altConf); //nextTarget[17]='\0'; if(strcmp(prevTarget,nextTarget) == 0) //free the memory for the duplicate dot { next = i->next; i->next = next->next; free(next); next=NULL; i->dotCount++; } else { i=i->next; i->dotCount=1; } }//removing duplicate loop ends curEnd = i; //updating the end of the dots of the source atom }/*code for removing dots ends*/ //integrating the condensed list of dots of the current source atom to prev and next source atoms curEnd->next = nextStart; if(prevEnd) prevEnd->next = curStart; else head=curStart; prevEnd=curEnd; //making the current source atom as the prev source atom }/* main for loop*/ return head; } /*}}}SortdotNodes()____________________________________________________________*/ /*{{{enumerate()**************************************************************/ void enumerate(FILE *outf, char* groupname, dotNode *results[][NODEWIDTH], float probeRad, int method, int nsel, int spike, int outdots, int numSkinDots, float density) { int i, j, doit; float gs, hs, bs, hslen, bslen, tgs, ths, tbs, thslen, tbslen, psas, tsas; float dtgp, score, tGscore, tHscore, tBscore, tscore, scoreValue, a_radius; double scaledGap, slen; dotNode *node; /* char *mast[NODEWIDTH] = {"wide_contact ", "close_contact ", "small_overlap ", "bad_overlap ", "H-bond "};*/ /* 04/10/2015 SJ changed to match the new categories, see NODEWIDTH in probe.h*/ char *mast[NODEWIDTH] = {"wide_contact","close_contact","weak_H-bond","small_overlap","bad_overlap","worse_overlap","H-bond"}; /* psas and tsas are the partial and total solvent accessible surface */ fprintf(outf, " \nsubgroup: %s\n", groupname); fprintf(outf, "atoms selected: %d\npotential dots: %d\npotential area: %.1f A^2\n", nsel, numSkinDots, numSkinDots/density); if (nsel <= 0 || numSkinDots <= 0) { fprintf(outf, "empty selection\n"); return; } if (spike) { fprintf(outf, " type # %% score score/A^2 x 1000\n"); } else { fprintf(outf, " type # %%\n"); } tgs = ths = thslen = tbs = tbslen = tsas = 0.0; tGscore = tHscore = tBscore = tscore = 0.0; for (i = 0; i < NUMATOMTYPES; i++) { for (j = 0; j < NODEWIDTH; j++) { gs = hs = hslen = bs = bslen = score = psas = 0.0; node = results[i][j]; doit = (node != NULL); if (doit) { fprintf(outf, "%3s %s ", getAtomName(i), outdots?"external_dots ":mast[j]); } while(node) { if (spike)/* 04/10/2015 SJ changed the if statements in this block to reflect change in dots ordering, see NODEWIDTH in probe.h*/ { // if (j == 0 || j == 1) { /* contact dot */ if(j==0 || j ==1 || j ==2) { /* contact dot, and weak H bond*/ gs += 1.0; dtgp = node->gap; scaledGap = dtgp/GAPweight; scoreValue = exp(-scaledGap*scaledGap); score += scoreValue; tGscore += scoreValue; } // else if (j == 2 || j == 3) { /* bump */ else if (j ==3 || j ==4 || j ==5){ /* bump*/ bs += 1.0; slen = 0.5*FABS(node->gap); bslen += slen; scoreValue = - BUMPweight * slen; score += scoreValue; tBscore += scoreValue; } else { /* H-bond */ hs += 1.0; slen = 0.5*FABS(node->gap); hslen += slen; scoreValue = HBweight * slen; score += scoreValue; tHscore += scoreValue; } } else { gs += 1.0; } if (method == EXTERNALSURFACE) { a_radius = node->a->radius; psas += (a_radius + probeRad)*(a_radius + probeRad)/(a_radius * a_radius); } node = node->next; } if (doit)/* 04/10/2015 SJ changed the if statements in this block to reflect change in dots ordering, see NODEWIDTH in probe.h*/ { if (spike) { // if (j == 0 || j == 1) { /* contact dot */ if(j==0 || j ==1 || j ==2) { /* contact dot, and weak H bond*/ fprintf(outf, "%7.0f %5.1f%% %9.1f %9.2f\n", gs, 100.0*gs/numSkinDots, score/density, 1000.0*score/numSkinDots); } // else if (j == 2 || j == 3) { /* bump */ else if (j ==3 || j ==4 || j ==5){ /* bump*/ fprintf(outf, "%7.0f %5.1f%% %9.1f %9.2f\n", bs, 100.0*bs/numSkinDots, score/density, 1000.0*score/numSkinDots); } else { /* H-bond */ fprintf(outf, "%7.0f %5.1f%% %9.1f %9.2f\n", hs, 100.0*hs/numSkinDots, score/density, 1000.0*score/numSkinDots); } } else { fprintf(outf, "%7.0f %5.1f%%\n", gs, 100.0*gs/numSkinDots); } tgs += gs; ths += hs; thslen += hslen; tbs += bs; tbslen += bslen; tscore += score; if (method == EXTERNALSURFACE) { tsas += psas; /* tally the solvent accessible surface */ } } } } if (spike) { /* fprintf(outf, " score/A^2 x 1000\n");*/ fprintf(outf, "\n tot contact: %7.0f %5.1f%% %9.1f %9.2f\n", tgs, 100.0*tgs/numSkinDots, tGscore/density, 1000.0*tGscore/numSkinDots); fprintf(outf, " tot overlap: %7.0f %5.1f%% %9.1f %9.2f\n", tbs, 100.0*tbs/numSkinDots, tBscore/density, 1000.0*tBscore/numSkinDots); fprintf(outf, " tot H-bond: %7.0f %5.1f%% %9.1f %9.2f\n", ths, 100.0*ths/numSkinDots, tHscore/density, 1000.0*tHscore/numSkinDots); fprintf(outf, "\n grand tot: %7.0f %5.1f%% %9.1f %9.2f\n", (tgs+tbs+ths), 100.0*(tgs+tbs+ths)/numSkinDots, tscore/density, 1000.0*tscore/numSkinDots); fprintf(outf, "\ncontact surface area: %.1f A^2\n", (tgs+tbs+ths)/density); } else { fprintf(outf, " tot: %7.0f %5.1f%%\n\n", tgs, 100.0*tgs/numSkinDots); fprintf(outf, " contact surface area: %.1f A^2\n", tgs/density); if (method == EXTERNALSURFACE) { fprintf(outf, "accessible surface area: %.1f A^2\n\n", tsas/density); } } } /*}}}enumerate()_____________________________________________________________*/ /*{{{rawEnumerate()***********************************************************/ void rawEnumerate(FILE *outf, char* groupname, dotNode *results[][NODEWIDTH], int method, int nsel, int spike, int outdots, int numSkinDots, float density, char *namestring, char *rawname, double scoreBias) {/*rawEnumerate*/ int i, j, doit; float gs, hs, bs, hslen, bslen, tgs, ths, tbs, thslen, tbslen; float dtgp, score, tGscore, tHscore, tBscore, tscore, scoreValue; double scaledGap, slen; dotNode *node; /*autobondrot: countDots==1, rawOutput==1 */ /* *rawname holds autobondrot angle values as space_char deliniated string*/ if (nsel <= 0 || numSkinDots <= 0) { /* empty selection */ if (method == EXTERNALSURFACE) { fprintf(outf, "%9.3f", 0.0); } else if (spike) { fprintf(outf, "%9.3f", scoreBias); } if (*rawname) { fprintf(outf, " %s", rawname); } if (*namestring || *groupname) { fprintf(outf, "%s", RAW_HEADER_COMMENT); if (*namestring) { fprintf(outf, " %s", namestring); } if (*groupname) { fprintf(outf, " %s", groupname); } } fprintf(outf, "\n"); return; } tgs = ths = thslen = tbs = tbslen = 0.0; tGscore = tHscore = tBscore = tscore = 0.0; for (i = 0; i < NUMATOMTYPES; i++) { for (j = 0; j < NODEWIDTH; j++) { gs = hs = hslen = bs = bslen = score = 0.0; node = results[i][j]; doit = (node != NULL); while(node) { if (spike)/* 04/10/2015 SJ changed the if statements in this block to reflect change in dots ordering, see NODEWIDTH in probe.h*/ { // if (j == 0 || j == 1) { /* contact dot */ if(j ==0 || j ==1 || j ==2){ /* contacts, and weak H bonds*/ gs += 1.0; dtgp = node->gap; scaledGap = dtgp/GAPweight; scoreValue = exp(-scaledGap*scaledGap); score += scoreValue; tGscore += scoreValue; } // else if (j == 2 || j == 3) { /* bump */ else if(j ==3 || j ==4 || j ==5){ /* bump */ bs += 1.0; slen = 0.5*FABS(node->gap); bslen += slen; scoreValue = - BUMPweight * slen; score += scoreValue; tBscore += scoreValue; } else { /* H-bond */ hs += 1.0; slen = 0.5*FABS(node->gap); hslen += slen; scoreValue = HBweight * slen; score += scoreValue; tHscore += scoreValue; } } else { gs += 1.0; } node = node->next; } if (doit) { tgs += gs; ths += hs; thslen += hslen; tbs += bs; tbslen += bslen; tscore += score; } } } /*output one line of information */ if (method == EXTERNALSURFACE) { fprintf(outf, "%9.3f", (tgs+tbs+ths)/density); } else if (spike) /*autobondrot: spike==1 at least in all examples dcr tested*/ { fprintf(outf, "%9.3f", scoreBias + (tscore/density)); /*value # ...*/ } else { fprintf(outf, "%9.3f", scoreBias + tgs); } if (*rawname) { fprintf(outf, " %s", rawname); } /*autobondrot angle values*/ if (*namestring || *groupname) { fprintf(outf, "%s", RAW_HEADER_COMMENT); if (*namestring) { fprintf(outf, " %s", namestring); } if (*groupname) { fprintf(outf, " %s", groupname); } } fprintf(outf, "\n"); /* NOTE this is the LF for above print statements*/ /* char* rawname == buf that is filled in autobondrot.c/runnameAndTorsion()*/ /* contains the autobondrot angle values in the order in which processed*/ }/*rawEnumerate*/ /*}}}rawEnumerate()__________________________________________________________*/ /*{{{countSelected()**********************************************************/ int countSelected(atom *theAtoms, int srcFlag) { atom *a = NULL; int ns = 0; for(a = theAtoms; a; a = a->next) { if (a->flags & srcFlag) { ns++; } } return ns; } /*}}}countSelected()_________________________________________________________*/ /*{{{enumDotSkin()************************************************************/ /* NOTE: allMainAtoms/abins & allMovingAtoms/bbins must be disjoint */ /* sets of atoms (none in common) or allMovingAtoms/bbins can */ /* be NULL. */ /* allMovingAtoms refers to autobondrot set of atoms */ int enumDotSkin(atom *allMainAtoms, atomBins *abins, atom *allMovingAtoms, atomBins *bbins, pointSet dots[], int srcFlag) { atom *src = NULL, *atomList = NULL, *atomList2 = NULL; int dotTotal = 0, usesMovingAtoms = FALSE; int dummy = TRUE; /*numSkinDots used to normalize output score*/ usesMovingAtoms = ((allMovingAtoms != NULL) && (bbins != NULL)); for(src = allMainAtoms; src; src = src->next) { if (src->flags & srcFlag) { atomList = findTouchingAtoms(src, NULL, abins, 0.0, 0, &dummy); if (usesMovingAtoms) { atomList2 = findTouchingAtoms(src, atomList, bbins, 0.0, 0, &dummy); } else { atomList2 = atomList; } if (atomList2) { markBonds(src, atomList2, 1, Maxbonded); /*in enumDotSkin()*/ if (Maxbonded > 3) { fixupLongBondChains(src, atomList2, 3); } } dotTotal += countSkin(src, atomList2, dots); if(Verbose && ShowTicks) { fprintf(stderr, "%s%s \r", src->r->resname, src->r->Hy36resno); } } } if (usesMovingAtoms) { for(src = allMovingAtoms; src; src = src->next) { if (src->flags & srcFlag) { atomList = findTouchingAtoms(src, NULL, abins, 0.0, 0, &dummy); atomList2 = findTouchingAtoms(src, atomList, bbins, 0.0, 0, &dummy); if (atomList2) { markBonds(src, atomList2, 1, Maxbonded); /*in enumDotSkin()*/ if (Maxbonded > 3) { fixupLongBondChains(src, atomList2, 3); } } dotTotal += countSkin(src, atomList2, dots); if(Verbose && ShowTicks) { fprintf(stderr, "%s%s \r", src->r->resname, src->r->Hy36resno); } } } } return dotTotal; } /*}}}enumDotSkin()___________________________________________________________*/ /*{{{countSkin()**************************************************************/ int countSkin(atom *src, atom *scratch, pointSet dots[]) { atom *targ = NULL; int i = 0, within = 0, ok = 0; point3d dotloc, dotvect; pointSet *srcDots; int dotCnt = 0; if (src->elem == ignoreAtom || src->elem == atomHOd) { return 0; } /*hb-only-dummy, phantom H atom*/ srcDots = &(dots[src->elem]); if (src->elem == atomC && isCarbonylAtom(src)) { srcDots = &COdots; } for(i=0; i < srcDots->n; i++) { dotvect = srcDots->p[i]; v3add(&(src->loc), &dotvect, &dotloc); ok = TRUE; /*targ is an atom*, scratch is an atom*, atom* defined in abin.h */ /* atom* scratch is a member of atom*: which has member scratch...*/ /* so for-loop moves through all these atoms while an atom has a scratch*/ for(targ = scratch; targ; targ = targ->scratch) { if (targ->elem == ignoreAtom) {continue;} /* eliminate if within bonded atom! */ if (targ->mark && (targ->elem != atomHOd)) { #ifdef INLINE_FOR_SPEED within = INRANGEINLINE(dotloc, (targ->loc), targ->radius); #else within = inRange(&dotloc, &(targ->loc), targ->radius); #endif if (within) { ok = FALSE; break; } } } if (ok) { dotCnt++; } } return dotCnt; } /*}}}countSkin()_____________________________________________________________*/ /*{{{updateHydrogenInfo()*****************************************************/ /* Make polor hydrogens HB donors. */ /* Waters without hydrogens are protonated with "phantoms". */ /* If we are using moving atoms, all the new water H atoms go there */ /* can return a list of new atom "clones" which are just the MainAtoms waters */ /* NOTE: allMainAtoms/abins & allMovingAtoms/bbins must be disjoint */ /* sets of atoms (none in common) or allMovingAtoms/bbins can */ /* be NULL. */ /* allMovingAtoms refers to autobondrot set of atoms */ /*usual call: (i.e. NOT autobondrot) updateHydrogenInfo(outf, allMainAtoms, abins, NULL, NULL, SET1|SET2, FALSE); */ atom* updateHydrogenInfo(FILE *outf, atom *allMainAtoms, atomBins *abins, atom *allMovingAtoms, atomBins *bbins, int selectedFlag, int mustSaveMainWater) {/*updateHydrogenInfo*/ atom *src = NULL, *orig = NULL, *atomList = NULL, *a = NULL; atom *newH = NULL, *tempStaticList = NULL, *mainWaterStorage = NULL; int type = 0, newHcount = 0, i = 0, whichList = 0; int nProx = 0, usesMovingAtoms = FALSE; int dummy = FALSE; /* Information specific to making sure only one H2O Hydrogen is */ /* generated which points into ring atoms on a given aromatic ring. */ struct {atom *a; float gap;} proximity[20]; usesMovingAtoms = ((allMovingAtoms != NULL) && (bbins != NULL)); if (! usesMovingAtoms) { allMovingAtoms = NULL; bbins = NULL; } if (DumpNewHO) { fprintf(outf,"@vectorlist {water H?} color= gray\n"); } tempStaticList = NULL; /* new hydrogens for MainAtoms when usesMovingAtoms */ mainWaterStorage = NULL; /* MainAtoms waters when usesMovingAtoms */ for(whichList = 0; whichList <= 1; whichList++) { /* loop over both lists of atoms */ for(src = ((whichList==0)?allMainAtoms:allMovingAtoms); src; src = src->next) {/*loop over src atoms*/ /* Note: we have to examine ALL the atoms at this point. */ /* Not testing (src->flags & selectedFlag) for H will slow down */ /* the code but it is neccessary to get the atom types correct. */ if (isHatom(src->elem)) { /* fixup Hs in the input */ atomList = NULL; if (src->flags & selectedFlag) { atomList = findTouchingAtoms(src, NULL, abins, 0.0, 0, &dummy); if (usesMovingAtoms) { atomList = findTouchingAtoms(src, atomList, bbins, 0.0, 0, &dummy); } } if (atomList) { markBonds(src, atomList, 1, 1); /*in updateHydrogenInfo()*/ type = dotType(src, atomList, ByNABaseColor); if ( type == atomN || type == atomO || type == atomS ) {/* bgn: if (atom types) */ if(Verbose && ShowTicks) { fprintf(stderr, "%s%d \r",src->r->resname, src->r->resid); } src->props |= DONOR_PROP; /* polar hydrogens are donors */ if (UsePolarH) { src->elem = atomHpolar; /* very important */ src->radius = getRadius(src->elem, 0); } /* there are a more restrictive set of cutoffs in reduce: 0.7 & 40 */ #define H2OoccCUTTOFF 0.25 #define H2ObvalCUTTOFF 80.0 if ( ( FABS(src->occ) < H2OoccCUTTOFF || src->bval >= H2ObvalCUTTOFF ) &&(src->props & WATER_PROP) ) { src->elem = ignoreAtom; src->radius = 0.0; } /* ignore low occupancy High B water hydrogens */ for(a = atomList; a; a = a->scratch) { /* since we found a hydrogen, then we know */ /* that the heavy atom is not the donor */ /* so we can disable that property */ if( (a->mark == 1) && !(isHatom(a->elem)) && (FABS(src->occ) >= H2OoccCUTTOFF) && (src->bval < H2ObvalCUTTOFF) ) { a->props &= ~DONOR_PROP; /*set to one's complement of ...*/ } } }/* end: if (atom types) */ } }/* fixup Hs in the input */ else if ( (src->flags & selectedFlag) /* fixup water oxygen */ && (src->props & WATER_PROP) && (src->elem == atomO) ) {/* else if water */ /* fixup water oxygen */ /*if connected H already found, the DONOR prop will not be set*/ /*if it is, we look further ahead in this residue for Hs */ if (src->props & DONOR_PROP) { for(a = src->next; a && (src->r == a->r); a = a->next) { if ( isHatom(a->elem) && (src->altConf == a->altConf) && (a->occ > 0.1)) { src->props &= ~DONOR_PROP; /* turn off donor property */ } } } /* if donor flag still set... */ if (src->props & DONOR_PROP) {/* water O without occupied Hs */ /* begin clone procedure */ /* special case for when we want to remember just these waters */ /* in this state */ if (mustSaveMainWater && usesMovingAtoms && (whichList==0)) { atom* waterClone = NULL; waterClone = (atom *)malloc(sizeof(atom)); if (waterClone) { *waterClone = *src; waterClone->next = mainWaterStorage; /* add clone to list */ mainWaterStorage = waterClone; } } /* end clone procedure */ src->props |= AMBIGWATER_PROP; /* this is a water which needs H? added */ src->props |= ACCEPTOR_PROP; src->props &= ~DONOR_PROP; /* acceptor, not donor */ orig = src; /* keep track of the original Oxygen atom */ atomList = findTouchingAtoms(src, NULL, abins, 0.25, 0, &dummy); if (usesMovingAtoms) { atomList = findTouchingAtoms(src, atomList, bbins, 0.25, 0, &dummy); } /* the 0.25 comes from: OxygenVDW(1.4) + 2*probeRad(0.3) == HpolVDW(1.0) + OHBondLen(1.0) and we want a minimum overlap of 0.1 before we consider possible Hs */ /* for each nearby HB acceptor, keep only 1 per aromatic ring */ nProx = 0; for(a = atomList; a; a = a->scratch) { if (a->props & ACCEPTOR_PROP) {/* bgn: if acceptor_prop */ float hbgap = gapSize(&(a->loc), &(orig->loc), a->radius + 2.0f); int foundAromRingAtom = FALSE; if (a->props & AROMATIC_PROP) { for (i=0; i < nProx; i++) { if ( a->r == proximity[i].a->r && proximity[i].a->props & AROMATIC_PROP) { if (hbgap < proximity[i].gap) { proximity[i].a = a; proximity[i].gap = hbgap; } foundAromRingAtom = TRUE; break; } } } if ( (!foundAromRingAtom) && nProx < sizeof(proximity)/sizeof(proximity[0])) { proximity[nProx].a = a; proximity[nProx].gap = hbgap; ++nProx; } }/* end: if acceptor_prop */ } /* build new HOH hydrogens in the direction of each acceptor */ newHcount = 0; for (i=0; i < nProx; i++) {/*scan over neighbors*/ a = proximity[i].a; if (a->props & ACCEPTOR_PROP) {/*close acceptor! Make a phantom H in the direction of each Acceptor*/ newH = (atom *)malloc(sizeof(atom)); if (newH) {/*newH allocated, set its parameters*/ point3d o2aVec; double waterOHbondLen; *newH = *orig; if (usesMovingAtoms) { newH->next = tempStaticList; /* hook into temp list */ tempStaticList = newH; } else { newH->next = src->next; /* hook into list of atoms */ src->next = newH; } v3sub(&(a->loc), &(orig->loc), &o2aVec); #define BEST_HBOND_OVERLAP 0.6 waterOHbondLen = 1.0 + MAX(-1.0, MIN(0.0, proximity[i].gap + BEST_HBOND_OVERLAP)); v3scale(&o2aVec, waterOHbondLen); v3add(&(orig->loc), &o2aVec, &(newH->loc)); newH->nextInBin= NULL; /* added to bins below */ newH->scratch = NULL; newH->elem = atomHOd; /*hb-only-dummy, phantom H atom*/ newH->radius = getRadius(newH->elem, 0); newH->covRad = getCovRad(newH->elem); /*atomprops.h/atomProp AtomTbl covRad*/ strcpy(newH->atomname, " H? "); newH->props |= DONOR_PROP; newH->props &= ~ACCEPTOR_PROP; /* donor, not acceptor */ /* *** adding an atom not in input! *** */ addNeighbor(newH, (usesMovingAtoms ? bbins : abins)); if (DumpNewHO) { if (DebugLevel>3) { fprintf(stderr, "HETATM%5d H%-2d%c%3.3s%s%4.4s%c %8.3f%8.3f%8.3f%6.2f%6.2f new H\n", 0, ++newHcount, newH->altConf, newH->r->resname, newH->r->chain, newH->r->Hy36resno, newH->r->resInsCode, newH->loc.x, newH->loc.y, newH->loc.z, newH->occ, newH->bval); } fprintf(outf, "{%4.4s%c%3.3s%s%4.4s%c}P %8.3f%8.3f%8.3f\n", orig->atomname, orig->altConf, orig->r->resname, orig->r->chain, orig->r->Hy36resno, orig->r->resInsCode, orig->loc.x, orig->loc.y, orig->loc.z); fprintf(outf, "{%4.4s%c%3.3s%s%4s%c}L %8.3f%8.3f%8.3f\n", newH->atomname, newH->altConf, newH->r->resname, newH->r->chain, newH->r->Hy36resno, newH->r->resInsCode, newH->loc.x, newH->loc.y, newH->loc.z); } if (! usesMovingAtoms) { src = src->next; /* bypass new atom when we iterate */ /* when inserted after src atom */ } }/*newH allocated, set its parameters*/ }/*close acceptor! Make a phantom H in the direction of each Acceptor*/ }/*scan over neighbors*/ }/* water O without occupied Hs */ }/* else if water */ /* fixup water oxygen */ /* non-water polar (N/O/S) atoms in selection list */ else if ( (src->flags & selectedFlag) && !(src->props & WATER_PROP) && (src->elem == atomN || src->elem == atomO || src->elem == atomS)) { /* Remove any DONOR property for non-water non-H */ /* polar atoms since the Hydrogen atoms will be the DONORs. */ /* This was assigned in select.c to facilitate selections. */ src->props &= ~DONOR_PROP; } }/*loop over src atoms */ /* end loop over list of src atoms*/ }/* loop over both lists of atoms */ if (tempStaticList) {/* new hydrogens for MainAtoms when usesMovingAtoms */ for(src = allMovingAtoms; src; src = src->next) { if (src->next == NULL) { src->next = tempStaticList; break; /* make sure we don't continually loop */ } } } return mainWaterStorage; }/*updateHydrogenInfo*/ /*}}}updateHydrogenInfo()____________________________________________________*/ /*{{{dump_changes()***********************************************************/ void dump_changes(FILE *outf) { fprintf(outf, "Probe change log:\n"); fprintf(outf, "\n"); fprintf(outf, "Note: Not captured prior to Aug 1998.\n"); fprintf(outf, " JMW = J Michael Word\n"); fprintf(outf, "\n"); fprintf(outf, " 8/ 7/98 - JMW - added segID parsing for XPLOR compatibility\n"); fprintf(outf, " 8/10/98 - JMW - added formatting for O and XtalView output\n"); fprintf(outf, "11/12/98 - JMW - use table to parent H atoms in std. residues,\n"); fprintf(outf, " made GAPweight expand with expanded radius,\n"); fprintf(outf, " and reorganized O format output in chunks\n"); fprintf(outf, " and changed xview format colors and modified\n"); fprintf(outf, " hydrogen name parsing for HG## and HE##\n"); fprintf(outf, "11/15/98 - JMW - added -stdbonds and -noparent flags and built\n"); fprintf(outf, " a complete std AA and NA connectivity table\n"); fprintf(outf, "11/23/98 - JMW - extended std bond table for truncated H names\n"); fprintf(outf, " 3/ 8/99 - JMW - atom names can include * and ' characters\n"); fprintf(outf, " 3/11/99 - JMW - fixed HBond problem with implicit hydrogens\n"); fprintf(outf, " and altered the skipping of MCMC to skip them\n"); fprintf(outf, " only if within a single chain and not a HET\n"); fprintf(outf, " 3/16/99 - JMW - Made water donor/acceptor insensitive to names\n"); fprintf(outf, " and extended atom name parsing for wierd HETS\n"); fprintf(outf, " 4/05/99 - JMW - Fixed typo in error checking in parse.c for\n"); fprintf(outf, " processing within _ of _,_,_\n"); fprintf(outf, " 4/06/99 - JMW - Updated stdbond table for NT and OXT\n"); fprintf(outf, " 4/07/99 - JMW - Cleaned up compiler warnings on initialization\n"); fprintf(outf, " and unused variables\n"); fprintf(outf, " 5/10/99 - JMW - Added AIB,ABU and MSE to mc hbond NH & O list\n"); fprintf(outf, " and added HN to the list of aliases for H\n"); fprintf(outf, " 7/12/99 - JMW - Begin changes for version 2\n"); fprintf(outf, " 7/17/99 - JMW - Added support for autobondrot, display ref.\n"); fprintf(outf, " 8/ 3/99 - JMW - Fixed typo in UpdateHydrogenInfo about bins\n"); fprintf(outf, " 8/ 5/99 - JMW - Sort autobondrot atoms in residue order,\n"); fprintf(outf, " add include file processing to autobondrot,\n"); fprintf(outf, " add -notick and -kinemage probe options\n"); fprintf(outf, " list process Nt/Ct/res for moving atoms\n"); fprintf(outf, " 8/16/99 - JMW - Flag -kinemage forces kin format (for -auto)\n"); fprintf(outf, " added CA-CB for gly in stdbonds for mutants\n"); fprintf(outf, " 9/ 2/99 - JMW - Added stdbond info to support alternate name\n"); fprintf(outf, " C5A for the methyl on a THY nucleic acid base\n"); fprintf(outf, "10/ 5/99 - JMW - Dropped unused variables flagged by compiler\n"); fprintf(outf, "12/ 6/99 - JMW - -OUT shows the contact surface area and the\n"); fprintf(outf, " solvent accessible surface area\n"); fprintf(outf, " -ADDVDW for radius offset(-SCALEVDW for scale)\n"); fprintf(outf, " -OUTColor sets the point color for -OUT\n"); fprintf(outf, "12/13/99 - JMW - added shortcuts: -SCSurface, -CONTACT,\n"); fprintf(outf, " -SASurface, -ACCESS\n"); fprintf(outf, "12/13/99 - JMW - renamed: -SASurface to -ASurface\n"); fprintf(outf, " 1/ 3/00 - JMW - added names to surface shortcuts\n"); fprintf(outf, " 1/13/00 - JMW - renamed: -CONTACT to -EXPOsed\n"); fprintf(outf, " 1/14/00 - JMW - fixed TRANS keyword in autobondrot\n"); fprintf(outf, " 5/ 4/00 - JMW - changed -u output to print spike coords and\n"); fprintf(outf, " added -oldu flag to give the old output\n"); fprintf(outf, " added -ignore flag to drop data from input\n"); fprintf(outf, " 5/31/00 - JMW - Added modified bases used in tRNAs\n"); fprintf(outf, " Included new properties TEST_ACCEPT_ANGLE and\n"); fprintf(outf, " CH_DONOR\n"); fprintf(outf, " 6/ 9/00 - JMW - constructed chain of residues view of atomlist\n"); fprintf(outf, " in order to support ring normals\n"); fprintf(outf, " 7/ 5/00 - JMW - inlined inRange for speed, remembered inosine,\n"); fprintf(outf, " dropped dummy atoms, added -segid flag,\n"); fprintf(outf, " fixed bug when water was used with -auto\n"); fprintf(outf, " 7/27/00 - JMW - added code to free extra memory at end of run\n"); fprintf(outf, "10/31/00 - JMW - added O2* & O2' to NABackbone list (oversite)\n"); fprintf(outf, "11/ 6/00 - JMW - changed malloc.h and memory.h to stdlib.h\n"); fprintf(outf, "11/14/00 - JMW - added -basecolor & -colorbase to color DNA/RNA\n"); fprintf(outf, "11/19/00 - JMW - fixed naBase typo for 7MG and added -nofaceHB\n"); fprintf(outf, "11/27/00 - JMW - updated properties of odd bases in select.c\n"); fprintf(outf, "11/28/00 - JMW - added H2U to base class table and moved PSU\n"); fprintf(outf, "12/01/00 - JMW - fixed H atomclass bug with -colorbase\n"); fprintf(outf, " 2/14/01 - JMW - changed LowGoodCut from -0.2 to -0.4\n"); fprintf(outf, " 7/11/01 -(v2.4?) JMW - dots no longer broken out by element in .kin,\n"); fprintf(outf, " added -element to generate old style output\n"); fprintf(outf, " which made kinemage buttons for elements. \n"); fprintf(outf, " 7/25/01 -(v2.5?) JMW - negative residue numbers in patterns \n"); fprintf(outf, " (Sometimes need to include extra parentheses or a \n"); fprintf(outf, " space to prevent the selection from being treated as\n"); fprintf(outf, " a command line option.) \n"); fprintf(outf, "10/ 1/01 - JMW - allow Nterminal fragnents to Hbond\n"); fprintf(outf, "10/ 9/01 - JMW - take out Nterm code and make CHX HBs an option\n"); fprintf(outf, " 1/23/03 - 2.09 - JMW - Run default now:-3 -mc -het -self \"altA ogt33\"\n"); fprintf(outf, " was: -het -self \"altA\"\n"); fprintf(outf, " Default for -ignore is now \"olt1\"\n"); fprintf(outf, " Recognise insert codes like 123a, 2B-3, & insC\n"); fprintf(outf, " Fixed bounding box water-H bug in autobondrot!\n"); fprintf(outf, "10/14/03 - 2.10 - JMW - Fixed bug (found by Patrick Reardon) caused by\n"); fprintf(outf, " filtering out non-target contacts too early \n"); fprintf(outf, " Also made sure atoms in different models\n"); fprintf(outf, " can not interact in any way.\n"); fprintf(outf, " Also changed how alternate conformations in\n"); fprintf(outf, " in different residues interact. Now altA and\n"); fprintf(outf, " altB never see one another. Previously, they\n"); fprintf(outf, " could if alts A and B weren't the same res.\n"); fprintf(outf, " Added -minocc flag with a default of 0.02\n"); fprintf(outf, " as the minimum non-zero occupancy.\n"); fprintf(outf, " Added -changes flag to dump change log.\n"); fprintf(outf, "Oct 2004: DCR code annotations, fixups, default modifiations\n"); fprintf(outf, "10/07/04 - 2.10 dcr041007 now presumes atom named _SE_ is Selenium\n"); fprintf(outf, " i.e. refmac or cns missadjusted name.\n"); fprintf(outf, "halides flagged as negative ions and H-bond acceptors\n"); fprintf(outf, "default: probe infile > outfile \n"); fprintf(outf, "same as: probe -3 -mc -het -self \"altA ogt33\" infile > outfile \n"); fprintf(outf, "try : probe -4H -mc -het -self \"altA ogt33\" infile > outfile \n"); fprintf(outf,"dcr041009 pointmasters M==McMc, S==ScSc, O==Other(McSc,Het--)\n"); fprintf(outf,"dcr041010 -onlybadout option for Bad clashes only\n"); fprintf(outf,"dcr041017 default -4H, pointtmaster P==McSc, report counts\n"); fprintf(outf,"dcr041020 master= {surface}, more work on count reports\n"); fprintf(outf,"dcr041020 annotated (NOT implemented) writeHbonds JACKnANDREA\n"); fprintf(outf,"dcr041023 2D array of count information\n"); fprintf(outf,"dcr041026 NOT report count IntersectBothWays re probe update\n"); fprintf(outf,"dcr041101 -summary colon deliminated multi-count report\n"); fprintf(outf,"dcr041101 -oneline summary all on oneline\n"); fprintf(outf,"dcr041101 -DEFAULTs same as: <>, \n"); fprintf(outf," but allows flags like -summary or -oneline\n"); /*jmw & dcr agreement of 041110 on version name and maintenance by dcr*/ fprintf(outf,"2.11.041112 by agreement of 041110: maintained now by dcr\n"); fprintf(outf,"041112 -nomodeltest for mage/probe (Lmodeltest) see 050121\n"); fprintf(outf,"041114 more fussing with NMR models\n"); fprintf(outf,"050111,17,18,19 flow, esp. autobondrot, annotations...\n"); fprintf(outf,"050119 very low occ atoms have presence but not show contacts\n"); fprintf(outf,"050121 remove -nomodeltest stuff, mage now sends model # \n"); fprintf(outf,"060129 jEdit type fold-comments on each subroutine \n"); fprintf(outf," single vdw contact button replaces both wide and small \n"); fprintf(outf," -mastername flags extra master={name} (default: dots)\n"); fprintf(outf,"060212 something about hets also marked prot,rna,dna... \n"); fprintf(outf,"060831 -NOOCC atoms NOT filtered by occ value \n"); fprintf(outf,"061018 not treat HIS as aromatic ACCEPTOR \n"); fprintf(outf,"2.12.0708xx ... version number change (approx date) \n"); fprintf(outf,"070801 rmi Updated for wwPDB remediation and PDB v3.0 \n"); fprintf(outf,"070801 Recognizes both new and old heavy atom names \n"); fprintf(outf," Recognizes both new and old residue names for DNA \n"); fprintf(outf," Builds atom connectivity and is able to generate contact dots on \n"); fprintf(outf," Xplor, PDBv2.3, PDBv3.0, or files with mixed format even within a \n"); fprintf(outf," single residue \n"); fprintf(outf,"070821 add strcasecmp to utility.c rwgk \n"); fprintf(outf,"070905 add alternate names of thymine methyl \n"); fprintf(outf,"070913 add detail on probe unformatted to help \n"); fprintf(outf,"070913 fixed handling of mixed RNA files \n"); fprintf(outf,"071010 added support for Hybrid36 atom and residue numbers \n"); fprintf(outf,"071025 added support of Coot style RNA names \n"); fprintf(outf,"071101 added support for two character chain names \n"); fprintf(outf,"071107 added support for OT1 and OT2 as C-term carboxylate oxygens \n"); fprintf(outf,"071128 bug fix in parsing of command line chainIds \n"); fprintf(outf,"110413 bug fix onlybadout outputs only bo, not cc gjk\n"); fprintf(outf,"110830 MacOSX10.6 GCC finicky re. no format string in \n"); fprintf(outf," writeOutput function for color \n"); fprintf(outf,"110909 changed the elseif portion of coloroutput to have no extra space \n"); fprintf(outf,"111028 (v2.6?) Swati Jain (SJ) -condense flag, used with -u option \n"); fprintf(outf," gives one line per source atom - target atom pair, also \n"); fprintf(outf," gives number of dots in that interaction, as an additional column.\n"); fprintf(outf,"2.14.111204 DCR OneDotEach for atoms across gap on line of sight.\n"); fprintf(outf," Logic is quite different so separate examineDots subroutines.\n"); fprintf(outf," LOneDotEach option calls examineOneDotEach() subroutine.\n"); fprintf(outf," General cleanup of notes and rearrangements for code readability.\n"); fprintf(outf,"111205-120109 DCR tuning ODE, methylene defined in select.c\n"); fprintf(outf,"111028 - 1303xx tweaks and new param to trunk, mainly Jeff \n"); fprintf(outf,"130427 dcr reconciliation of trunk <- probeVector incl ODE\n"); fprintf(outf,"2.15.130427 version number change for merged code! \n"); fprintf(outf,"2.16.130509 jjh added support for segid instead of chaind\n"); fprintf(outf,"2.16.130520 jjh fixed bug in segid handling\n"); fprintf(outf,"04/16/2015 - SJ added the -sepworse flag, if true will seperate the overlaps of >= 0.4 and overlaps of >=0.5. This is default by false. Had to change NODEWIDTH value (see probe.h)\n One can change the 0.5 cutoff for worse overlap by specifying the flag DIVWorse\n"); fprintf(outf,"09/24/2021 - RMT Version 2.17 fixed crash when finding ambiguous Oxygens\n"); fprintf(outf,"10/05/2021 - RMT Version 2.18 makes default C=O radius scale depend on table value\n"); fprintf(outf,"12/09/2021 - RMT Version 2.19 adds commend-line argument to dump atom info\n"); fprintf(outf,"12/21/2021 - RMT Version 2.20 dumps atom information after hydrogens have been updated.\n"); fprintf(outf,"12/21/2021 - RMT Version 2.21 does not dump Phantom Hydrogen information.\n"); fprintf(outf,"03/31/2022 - RMT Version 2.22 does not mark Nitrogens in the HIS ring or atoms in the\n"); fprintf(outf," TRP 5-membered ring as acceptors. Too many side contacts are made to these rings.\n"); fprintf(outf,"04/04/2022 - RMT Version 2.23 fixes copy/paste bug to allow F and G pointmasters.\n"); fprintf(outf,"09/29/2022 - RMT Version 2.24 sorts raw dots by increasing gap, providing smallest for condensed output.\n"); fprintf(outf,"02/08/2023 - RMT Version 2.25 makes the Nitrogens in TRP and HIS 5-membered rings acceptors.\n"); fprintf(outf,"02/11/2023 - RMT Version 2.26 makes the Nitrogen in TRP not an acceptor.\n"); exit(0); }/*dump_changes()*/ /*}}}dump_changes()__________________________________________________________*/ /*{{{countsummary()***********************************************************/ /*count arrays indecies: Nth pass, contact width (wide to overlap, Hbond,...)*/ /*N either 1 or 2 depending on direction of relationship */ /* Old: contact width 20111215dcr*/ /* 0 1 2 3 4 5 mcmc wide: mcmc close: mcmc small : mcmc bad: mcmc H-bond: mcmc H-weak scsc wide: scsc close: scsc small: scsc bad: scsc H-bond: scsc H-weak mcsc wide: mcsc close: mcsc small: mcsc bad: mcsc H-bond: mcsc H-weak other wide: other close: other small: other bad: other H-bond other H-weak wide sum: close sum: small sum: bad sum: H-bond sum H-weak sum */ /*5 contact width criteria include optional 6th index 5 weakHbond 20111215dcr*/ /* New: contact width 04/15/2015 SJ*/ /* 0 1 2 3 4 5 6 mcmc wide: mcmc close: mcmc H-weak: mcmc small: mcmc bad: mcmc worse: mcmc H-bond: scsc wide: scsc close: scsc H-weak: scsc small: scsc bad: scsc worse: scsc H-bond: mcsc wide: mcsc close: mcsc H-weak: mcsc small: mcsc bad: mcsc worse: mcsc H-bond: other wide: other close: other H-weak: other small: other bad: other worse: other H-bond wide sum: close sum: H-weak sum: small sum: bad sum: worse sum: H-bond sum: */ void countsummary(FILE *outf, char* modestr, int Lines, int Pass) /*dcr041101*/ {/* 04/21/2015 SJ: don't think this needs to check flags before printing, it would just print 0 is LworseOverlap, and LweakHbond is false*/ char message[200]; int N=0; /*index Number*/ int k=0; if(Pass == 1 || Pass == 0) {N = 0;} else {N = 1;} summcont[N][0] = mcmccont[N][0]+scsccont[N][0]+mcsccont[N][0]+othrcont[N][0]; summcont[N][1] = mcmccont[N][1]+scsccont[N][1]+mcsccont[N][1]+othrcont[N][1]; summcont[N][2] = mcmccont[N][2]+scsccont[N][2]+mcsccont[N][2]+othrcont[N][2]; summcont[N][3] = mcmccont[N][3]+scsccont[N][3]+mcsccont[N][3]+othrcont[N][3]; summcont[N][4] = mcmccont[N][4]+scsccont[N][4]+mcsccont[N][4]+othrcont[N][4]; summcont[N][5] = mcmccont[N][5]+scsccont[N][5]+mcsccont[N][5]+othrcont[N][5]; summcont[N][6] = mcmccont[N][6]+scsccont[N][6]+mcsccont[N][6]+othrcont[N][6]; /* added SJ 04/15/2015*/ /*20111215dcr weakHbond optional Nodewith criteria*/ if(Pass == 0) {/*expect another pass, so save these values, and clear arrays for next pass*/ for(k=0 ; k<=6 ; k++) /* changed from k<=5 to k<=6 SJ 04/15/2015*/ {/*copy values*/ mcmccont[1][k] = mcmccont[0][k]; scsccont[1][k] = scsccont[0][k]; mcsccont[1][k] = mcsccont[0][k]; othrcont[1][k] = othrcont[0][k]; summcont[1][k] = summcont[0][k]; } for(k=0 ; k<=6 ; k++)/* changed from k<=5 to k<=6 SJ 04/15/2015*/ {/*clear arrays*/ mcmccont[0][k] = 0; scsccont[0][k] = 0; mcsccont[0][k] = 0; othrcont[0][k] = 0; summcont[0][k] = 0; } } if(Pass == 2) {/*should have been a previous pass with values saved*/ for(k=0 ; k<=6 ; k++)/* changed from k<=5 to k<=6 SJ 04/15/2015*/ { mcmccont[0][k] = mcmccont[0][k] + mcmccont[1][k]; scsccont[0][k] = scsccont[0][k] + scsccont[1][k]; mcsccont[0][k] = mcsccont[0][k] + mcsccont[1][k]; othrcont[0][k] = othrcont[0][k] + othrcont[1][k]; summcont[0][k] = summcont[0][k] + summcont[1][k]; } } if(Pass == 1 || Pass == 2) {/*output*/ if(Lines == 1) /*for accummulating oneliners for comparisons dcr041101*/ { /*Old: 0 wide contact,1 close contact,2 small overlap,3 bad overlap,4 H-bonds*/ /*Old: INDEX 5 NOW USED FOR OPTIONAL 6TH CONTACT CRITERIA OF WEAKHBONDS*/ /* New:0 wide contact,1 close contact,2 weak H-bonds, 3 small overlap,4 bad overlap,5 worse overlap, 6 H-bonds - SJ 04/15/2015*/ // added the [0][6] to all SJ 04/15/2015 fprintf(outf,": %s " ":%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld " ":%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld " ":%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld " ":%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld " ":%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :\n" ,inputfilename ,mcmccont[0][0],mcmccont[0][1],mcmccont[0][2],mcmccont[0][3],mcmccont[0][4] ,mcmccont[0][5],mcmccont[0][6] ,scsccont[0][0],scsccont[0][1],scsccont[0][2],scsccont[0][3],scsccont[0][4] ,scsccont[0][5],scsccont[0][6] ,mcsccont[0][0],mcsccont[0][1],mcsccont[0][2],mcsccont[0][3],mcsccont[0][4] ,mcsccont[0][5],mcsccont[0][6] ,othrcont[0][0],othrcont[0][1],othrcont[0][2],othrcont[0][3],othrcont[0][4] ,othrcont[0][5],othrcont[0][6] ,summcont[0][0],summcont[0][1],summcont[0][2],summcont[0][3],summcont[0][4] ,summcont[0][5],summcont[0][6] ); } else /*for kinemage*/ { fprintf(outf,"@text\n"); fprintf(outf,"probe: %s\n",modestr); fprintf(outf,"%s\n",inputfilename); /*Old: 0 wide contact,1 close contact,2 small overlap,3 bad overlap,4 H-bonds*/ //fprintf(outf,":CONTACT: WIDE : CLOSE : SMALL : BAD : H-BOND : weakHbond :\n"); /*New: SJ 04/15/2015 0 wide contact,1 close contact,2 weak H-bonds, 3 small overlap,4 bad overlap,5 worse overlap, 6 H-bonds*/ /* added [0][6] for all SJ 04/21/2015*/ fprintf(outf,":CONTACT: WIDE : CLOSE : weak H-bonds : SMALL : BAD : WORSE : H-BOND :\n"); fprintf(outf,":MCMC :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :\n" ,mcmccont[0][0],mcmccont[0][1],mcmccont[0][2],mcmccont[0][3],mcmccont[0][4] ,mcmccont[0][5],mcmccont[0][6]); fprintf(outf,":SCSC :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :\n" ,scsccont[0][0],scsccont[0][1],scsccont[0][2],scsccont[0][3],scsccont[0][4] ,scsccont[0][5],scsccont[0][6]); fprintf(outf,":MCSC :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :\n" ,mcsccont[0][0],mcsccont[0][1],mcsccont[0][2],mcsccont[0][3],mcsccont[0][4] ,mcsccont[0][5],mcsccont[0][6]); fprintf(outf,":OTHER :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :\n" ,othrcont[0][0],othrcont[0][1],othrcont[0][2],othrcont[0][3],othrcont[0][4] ,othrcont[0][5],othrcont[0][6]); fprintf(outf,":SUM :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :\n" ,summcont[0][0],summcont[0][1],summcont[0][2],summcont[0][3],summcont[0][4] ,summcont[0][5],summcont[0][6]); } if (Verbose) { sprintf(message,"%s",inputfilename); note(message); /*Old: 0 wide contact,1 close contact,2 small overlap,3 bad overlap,4 H-bonds*/ //sprintf(message,":CONTACT: WIDE : CLOSE : SMALL : BAD : H-BOND : weakHbond :"); /*New: SJ 04/15/2015 0 wide contact,1 close contact,2 weak H-bonds, 3 small overlap,4 bad overlap,5 worse overlap, 6 H-bonds*/ /* added [0][6] for all SJ 04/21/2015*/ sprintf(message,":CONTACT: WIDE : CLOSE : weak H-bonds : SMALL : BAD : WORSE : H-BOND :\n"); note(message); sprintf(message,":MCMC :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :" ,mcmccont[0][0],mcmccont[0][1],mcmccont[0][2],mcmccont[0][3],mcmccont[0][4] ,mcmccont[0][5],mcmccont[0][6]); note(message); sprintf(message,":SCSC :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :" ,scsccont[0][0],scsccont[0][1],scsccont[0][2],scsccont[0][3],scsccont[0][4] ,scsccont[0][5],scsccont[0][6]); note(message); sprintf(message,":MCSC :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :" ,mcsccont[0][0],mcsccont[0][1],mcsccont[0][2],mcsccont[0][3],mcsccont[0][4] ,mcsccont[0][5],mcsccont[0][6]); note(message); sprintf(message,":OTHER :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :" ,othrcont[0][0],othrcont[0][1],othrcont[0][2],othrcont[0][3],othrcont[0][4] ,othrcont[0][5],othrcont[0][6]); note(message); sprintf(message,":SUM :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :%9ld :" ,summcont[0][0],summcont[0][1],summcont[0][2],summcont[0][3],summcont[0][4] ,summcont[0][5],summcont[0][6]); note(message); } }/*output*/ } /*}}}countsummary()__________________________________________________________*/