#define MAXSTRANDS 500 #define MAXRES 1000 molecule m; residue r, r1, rholder; int ac, sn, rn; string str; string strands[ MAXSTRANDS ]; int nstr; int nres[ MAXSTRANDS ]; residue res[ MAXSTRANDS, MAXRES ]; int getcanonicalstackdist( file f, residue r, residue r1 ) { atom a, a1; float d, theta, rvl, r1vl; int number_of_violations; point rN1,rC4,rC2,rN3,rmid,rxmid,rxdir,rydir,rzdir,rv,r1v; point r1N1,r1C4,r1C2,r1N3,r1mid,r1xmid,r1xdir,r1ydir,r1zdir; point rP,r1P,rO3,r1O3; number_of_violations = 0; for( a in r ){ if (a.atomname =~ "P" ) rP = a.pos; if (a.atomname =~ "O3?" ) rO3 = a.pos; } for( a1 in r1 ){ if (a1.atomname =~ "P" ) r1P = a.pos; if (a1.atomname =~ "O3?" ) r1O3 = a.pos; } if( distp( rP, r1O3) <= distp( r1P, rO3 ) ){ rholder = r1; r1 = r; r = rholder; } for( a in r ){ if (a.atomname =~ "[HP'*]" ) continue; if ( a.atomname == "N1" ) rN1 = a.pos; else if ( a.atomname == "C2" ) rC2 = a.pos; else if ( a.atomname == "C4" ) rC4 = a.pos; else if ( a.atomname == "N3" ) rN3 = a.pos; for( a1 in r1 ){ if ( a1.atomname =~ "[HP'*]" ) continue; if ( a1.atomname == "N1" ) r1N1 = a1.pos; else if ( a1.atomname == "C2" ) r1C2 = a1.pos; else if ( a1.atomname == "C4" ) r1C4 = a1.pos; else if ( a1.atomname == "N3" ) r1N3 = a1.pos; d = distp( a.pos, a1.pos ); if ((d > 9.0 ) || ( d < 2.0 )){ number_of_violations++; } } } rmid = ( rN1 + rC4 )/2.0; rxmid = ( rC2 + rN3 )/2.0; if ( (r.resname == "ADE" ) || (r.resname == "A") || (r.resname == "GUA" ) || (r.resname == "G") ) rxdir = rxmid - rmid; else rxdir = rmid - rxmid; rydir = rC4 - rmid; rzdir = rxdir ^ rydir; rv = rzdir; rvl = sqrt( rzdir.x*rzdir.x + rzdir.y*rzdir.y + rzdir.z*rzdir.z ); r1mid = ( r1N1 + r1C4 )/2.0; r1xmid = ( r1C2 + r1N3 )/2.0; if ( (r1.resname == "ADE" ) || (r1.resname == "A") || (r1.resname == "GUA" ) || (r1.resname == "G") ) r1xdir = r1xmid - r1mid; else r1xdir = r1mid - r1xmid; r1ydir = r1C4 - r1mid; r1zdir = r1xdir ^ r1ydir; r1v = r1zdir; r1vl = sqrt( r1zdir.x*r1zdir.x + r1zdir.y*r1zdir.y + r1zdir.z*r1zdir.z ); theta = acos ( ( rv @ r1v ) / ( rvl * r1vl ) ); if ((number_of_violations < 10 ) && (theta < 20.0)){ for( a in r ){ for( a1 in r1 ){ d = distp( a.pos, a1.pos ); fprintf( f, "%s:%s-%s:%s %8.3f\n", r.resname, a.atomname, r1.resname, a1.atomname, d ); } } } }; for( ac = 2; ac <= argc; ac = ac + 1 ){ m = getpdb( argv[ ac ] ); str = ""; nstr = 0; rn = 0; for( r in m ){ if( str != r.strandname ){ nstr = nstr + 1; strands[ nstr ] = r.strandname; nres[ nstr ] = 0; } nres[ nstr ] = nres[ nstr ] + 1; res[ nstr, nres[ nstr ] ] = r; str = r.strandname; } for( sn = 1; sn <= nstr; sn = sn + 1 ){ for( rn = 1; rn < nres[ sn ]; rn = rn + 1 ){ r = res[ sn, rn ]; r1 = res[ sn, rn+1 ]; getcanonicalstackdist( stdout, r, r1 ); } } }