C PROTEIN DATA BANK SOURCE CODE BREAKRING C AUTHOR. C.MARZEC,L.DAY C ENTRY DATE. 1/94 UNSUPPORTED C LAST REVISION. 1/94 C PURPOSE. GENERATE PARAMETERS FROM 5-MEMBERED C PURPOSE. RING C LANGUAGE. FORTRAN C c Version 1.1 --- 6/25/91 --- by Christopher Marzec c c BREAKRING accepts five sets of (x,y,z) Cartesian coordinates of c the members of a ribose ring. It outputs the four internal c ring parameters [q, P, s, Gamma] and the bond lengths: c b(j) = |R(j+1)-R(j)| , where R(j) is atom position vector of atom j, c labelled so j=1 for C1; j=2 for C2; etc; j=5 for O4. c c BREAKRING is the inverse of MAKERING. The programs were written c on the basis of a full parameterization of five membered rings c described in C.J Marzec and L.A. Day (1993) "An exact description c of five-membered ring configurations. I. Parameterization via an c amplitude S, an angle Gamma, the pseudorotation amplitude, q, and phase c angle, P, and bond lengths." J. Biomolec. Struct. Dynam. 10, 1091-1123. c c BREAKRING calls S/R break, which performs the calculation. c Input is from keyboard or from file "makering.output" created by c program makering.f . Output is on screen and into file c "breakring.output". c implicit real*8 (a-h,o-z) dimension x(5),y(5),z(5),r(3,5),tauval(5),r1(3),r2(3),r3(3),r4(3) dimension r4(3),bond(5) character*1 answer character*2 atom(5) character*10 vertex(5) character*13 dihed(5) data atom/'C1','C2','C3','C4','O4'/ data dihed/'','','', 1 '',''/ data vertex/'','','', 1 '',''/ open(11,file='breakring.output') write(6,10) 10 format('Breakring accepts the x,y,z coordinates (in Angstroms) of') write(6,11) 11 format('five ribose ring members and outputs [q,P,S,Gamma; bondlengths].'/) c Read the (x,y,z) coords for each ring: C1, C2, C3, C4, O4 c Note: O4 is labelled with j = 5; C1 with j = 1; C2 with j = 2; etc. do 20 j = 1,5 write(6,22) atom(j) 22 format('Enter x,y,z coordinates for atom ',a2,':') read(5,*) x(j),y(j),z(j) r(1,j) = x(j) r(2,j) = y(j) r(3,j) = z(j) 20 write(6,25) x(j),y(j),z(j) 25 format(' (entry is ',3(1pd12.4),')') c S/R break makes the [q,P,s,gamma] calculation: call break(x,y,z,qval,pval,sval,gamma) write(6,72) qval,pval,sval,gamma 72 format(/'[q,P,S,Gamma] = ',/,'[',1pd11.4,' A,',1pd12.4,' deg.,', 1 1pd12.4,' A,',1pd12.4,' deg. ]') c c calculate bond lengths: write(6,202) do 120 j = 1,5 jp = j+1 if(jp.eq.6) jp = 1 bond(j) = dsqrt((r(1,jp)-r(1,j))**2+(r(2,jp)-r(2,j))**2 1 +(r(3,jp)-r(3,j))**2) 120 continue c print bond lengths: write(6,111) bond(1) write(6,112) bond(2) write(6,113) bond(3) write(6,114) bond(4) write(6,115) bond(5) 111 format(' bond length b1 = ',1pd12.4,' A.') 112 format(' bond length b2 = ',1pd12.4,' A.') 113 format(' bond length b3 = ',1pd12.4,' A.') 114 format(' bond length b4 = ',1pd12.4,' A.') 115 format(' bond length b0 = b5 = ',1pd12.4,' A.') c store [q,P,S,Gamma; bonds] in file "breakring.output": write(11,*) qval,pval,sval,gamma write(11,31) 31 format('bond lengths (b1...b5):') do 30 j = 1,5 30 write(11,*) bond(j) c write(6,200) 200 format(/'Do you wish to calculate the dihedrals? Enter y or n .') read(5,201) answer 201 format(a1) if(answer.ne.'y'.and.answer.ne.'Y') go to 299 202 format(/) write(6,202) c calculate ring atom dihedrals: do 210 j = 1,5 c calculate dihedral , where atom j = atom j+5. j1 = j+1 j2 = j+2 j3 = j+3 if(j1.gt.5) j1 = j1-5 if(j2.gt.5) j2 = j2-5 if(j3.gt.5) j3 = j3-5 do 220 i = 1,3 r1(i) = r(i,j) r2(i) = r(i,j1) r3(i) = r(i,j2) 220 r4(i) = r(i,j3) call dihedral(r1,r2,r3,r4,tauval(j1)) write(6,230) dihed(j1),tauval(j1) 230 format('Dihedral ',a13,' =',1pd12.4,' degrees.') 210 continue c compare to Pas pseudorotation calculation according to Altona c and Sundaralingam (1972) : tanpnum = tauval(4)+tauval(1)-tauval(3)-tauval(5) tanpden = 2.d0*tauval(2)*(dsin(4.d0*.314159265d0)+dsin(2.d0*.314159265d0)) pval = 180.d0*datan2(tanpnum,tanpden)/3.1415926530 if(pval.lt.0.d0) pval = pval + 360.d0 write(6,250) pval 250 format(/'For comparison, the Altona and Sundaralingam ' 1 'pseudorotation value, Pas, equals',(1pd12.4),' degrees.') 299 continue c write(6,400) 400 format(/'Do you wish to calculate bond angles? Enter y or n .') read(5,201) answer if(answer.ne.'y'.and.answer.ne.'Y') go to 499 write(6,202) do 410 j = 1,5 j1 = j+1 j2 = j+2 if(j1.gt.5) j1 = j1-5 if(j2.gt.5) j2 = j2-5 do 420 i = 1,3 r1(i) = r(i,j) r2(i) = r(i,j1) 420 r3(i) = r(i,j2) call angle(r1,r2,r3,ang) write(6,430) vertex(j1),ang 410 continue 430 format('Bond angle ',a10,' =',1pd12.4' degrees.') 499 continue write(6,999) 999 format(/'The calculated [q,P,S,Gamma; bondlengths] are 1 found in file "breakring.output".') c stop end subroutine break(x,y,z,qval,pcp,s,gamma) c breakring.f calculates the [q,P,s and gamma] values for a five-membered c ring, from input data consisting of [x(j),y(j),z(j)] Cartesian triples, c where j = 1..4 for C1...C4; and j=5 for O4 c implicit real*8 (a-h,o-z) dimension x(5),y(5),z(5),rp(3),rpp(3) dimension sinj(5),cosj(5),zunit(3),znew(5),xnew(5),ynew(5) dimension vecj(3),cos2j(5),sin2j(5),xunit(3),yunit(3),eta(5) dimension rcm(3,5),cm(3) complex z1,z2,z3,im c x(j,n) contains x val for atom j for ring n; j=5 for O4; c j=1 for C1; j=2 for C2; etc c data im/(0.d0,1.d0)/ twopi = 2.d0*3.14159265358979323846d0 d2r = twopi/360.d0 fac = dsqrt(2.d0/5.d0) c fill matrices used for Cremer/Pople analysis: do 14 j = 1,5 sinj(j) = dsin(twopi*j/5) cosj(j) = dcos(twopi*j/5) cos2j(j) = dcos(twopi*2*j/5) sin2j(j) = dsin(twopi*2*j/5) 14 continue c calculate the position of the center of mass (unit mases) c vector cm: do 2011 i = 1,3 2011 cm(i) = 0.d0 do 2010 j = 1,5 cm(1) = cm(1) + x(j)/5.d0 cm(2) = cm(2) + y(j)/5.d0 cm(3) = cm(3) + z(j)/5.d0 2010 continue c calculate the position vectors rcm in center-of-mass c coordinate system: rcm(j) = r(j)-cm(j) do 40 j = 1,5 rcm(1,j) = x(j) - cm(1) rcm(2,j) = y(j) - cm(2) rcm(3,j) = z(j) - cm(3) 40 continue c calculate rp (= A) and rpp (= B) according to Cremer and Pople: do 55 i = 1,3 rp(i) = 0.d0 rpp(i) = 0.d0 do 55 j = 1,5 rp(i) = rp(i) + rcm(i,j)*sinj(j) rpp(i) = rpp(i) + rcm(i,j)*cosj(j) 55 continue c now calculate the z unit vector: zunit = rp X rpp call cross (rp,rpp,zunit) c normalize zunit: call dot(zunit,zunit,zmag2) do 57 i = 1,3 c Use yunit = rp c yunit(i) = rp(i) yunit(i) = rpp(i) 57 zunit(i) = zunit(i)/dsqrt(zmag2) c normalize yunit: call dot(yunit,yunit,ymag2) do 570 i = 1,3 570 yunit(i) = yunit(i)/dsqrt(ymag2) c set xunit = yunit X zunit call cross(yunit,zunit,xunit) c check mag: call dot(zunit,zunit,zmag) call dot(zunit,rp,zxrp) call dot(zunit,rpp,zxrpp) call dot(yunit,yunit,ymag) call dot(xunit,xunit,xmag) c Now project out z coord along zunit direction: do 60 j = 1,5 do 61 i = 1,3 61 vecj(i) = rcm(i,j) call dot(vecj,zunit,zval) 60 znew(j) = zval c now calculate q,P values: qcosp = 0.d0 qsinp = 0.d0 do 70 j = 1,5 qcosp = qcosp + fac*znew(j)*cos2j(j) qsinp = qsinp - fac*znew(j)*sin2j(j) 70 continue phicp = datan2(qsinp,qcosp)/d2r qval = dsqrt(qcosp**2+qsinp**2) c convert Cremer/Pople angle phicp into the usual c pseudorotation angle pcp: pcp = phicp + 90.d0 if(pcp.gt.360.d0) pcp = pcp - 360.d0 if(pcp.lt.0.d0) pcp = pcp + 360.d0 30 format(f15.4,',',f15.4,',',i5) c Reconstruct the zj coord from q and P, as a check do 85 j = 1,5 zdif = znew(j)-fac*qval*dcos(phicp*d2r+.4d0*twopi*j) if(dabs(zdif).gt.1.d-12) write(6,86) j,zdif 86 format('Problem here: j,zdif =',i5,1pd12.4) 85 continue 888 format(/) c Project rcm(i,j) onto the mean C/P plane: do 200 j = 1,5 do 210 i = 1,3 210 vecj(i) = rcm(i,j) call dot(xunit,vecj,xval) call dot(yunit,vecj,yval) xnew(j) = xval ynew(j) = yval c Now ring atom coordinates are in (xnew(j),ynew(j),znew(j)) for j=1 to 5 200 continue c Rotate xnew and ynew st O4 lies on y-axis: ang5 = -twopi/4.d0 + datan2(ynew(5),xnew(5)) cosrot = dcos(ang5) sinrot = dsin(ang5) do 215 j = 1,5 xval = xnew(j) yval = ynew(j) xnew(j) = xval*cosrot + yval*sinrot 215 ynew(j) = -xval*sinrot + yval*cosrot c Ccalculate s and Gamma: c Use s**2 exp(2i*Gamma) = exp(-2i*eta(3))*SUMj bj**2 exp(2i*etaj) = z2 do 300 j = 1,5 c calculate etaj: 300 eta(j) = datan2(ynew(j),xnew(j)) z2 = cmplx(0.d0,0.d0) t11 = 0.d0 t22 = 0.d0 t12 = 0.d0 sumr = 0.d0 sumi = 0.d0 do 310 j = 1,5 sumr = sumr + (xnew(j)**2-ynew(j)**2) sumi = sumi + 2.d0*ynew(j)*xnew(j) bj2 = xnew(j)**2 + ynew(j)**2 z2 = z2 + bj2*cdexp(2.d0*im*eta(j)) t12 = t12 + xnew(j)*ynew(j) t11 = t11 + xnew(j)*xnew(j) t22 = t22 + ynew(j)*ynew(j) 310 continue c calculate s and gamma from sumr and sumi; these c refer to the coord system in which the plots are calculated. s4 = sumr**2 + sumi**2 sorig = s4**(.25d0) gammaorig = datan2(sumi,sumr) if(gammaorig.gt.twopi) gammaorig = gammaorig-twopi if(gammaorig.lt.0.d0) gammaorig = gammaorig + twopi gammaorig = .5d0*gammaorig/d2r write(6,888) 886 format('z2,sumr,sumi =',4(1pd12.4)) x2 = z2 y2 = -im*z2 y2 = 2.d0*t12 x2 = t11 - t22 rho4 = x2**2 + y2**2 s = rho4**(.25d0) twogamma = datan2(y2,x2)/d2r if(twogamma.lt.0.d0) twogamma = twogamma + 360.d0 if(twogamma.gt.360.d0) twogamma = twogamma - 360.d0 gamma = twogamma/2.d0 gammadif = dabs(gamma-gammaorig) sdif = dabs(s-sorig) if(gammadif.gt.1.d-3.or.sdif.gt.1.d-3) write(6,4555) gamma,gammaorig, 1 s,sorig 4555 format('Problem: gamma,gammaorig,s,sorig =',4(1pd12.4)) c check complex arithmetic: z1 = s**2*cdexp(2.d0*im*gamma*d2r) z3 = z1-z2 if(cabs(z3).gt.1.d-6) write(6,6534) z1,z2,z3 6534 format('Problem: z1,z2,z3 =',3(1pd12.4,2x,1pd12.4)) 100 continue return end subroutine cross(a,b,c) c c = a X b implicit real*8 (a-h,o-z) dimension a(3),b(3),c(3) c(3) = a(1)*b(2) - a(2)*b(1) c(1) = a(2)*b(3) - a(3)*b(2) c(2) = a(3)*b(1) - a(1)*b(3) return end subroutine dot (a,b,dotprod) c dotprod = a dot b implicit real*8 (a-h,o-z) dimension a(3),b(3) dotprod = a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end subroutine dihedral(r1,r2,r3,r4,tauvalnorm) implicit real*8 (a-h,o-z) dimension r1(3),r2(3),r3(3),r4(3),b1(3),b2(3),b3(3) real*8 norm1(3),norm2(3),nmag1,nmag2 c S/R dihedral accepts position vectors r1, r2, r3, and r4 c and calculates dihedral tauval = <1-2-3-4> in degrees. c First calculate bond vectors bj = r(j+1)-r(j) do 10 i = 1,3 b1(i) = r2(i) - r1(i) b2(i) = r3(i) - r2(i) 10 b3(i) = r4(i) - r3(i) c c Find tauvalnorm via cos(tau) = norm1 X norm2, where norm1 is c the normal to plane containing atoms 1, 2, and 3; and norm2 c is the normal to plane containing atoms 2, 3, and 4. call cross(b1,b2,norm1) call cross(b2,b3,norm2) c normalize: call dot(norm1,norm1,nmag1) nmag1 = dsqrt(nmag1) call dot(norm2,norm2,nmag2) nmag2 = dsqrt(nmag2) do 60 i = 1,3 norm1(i) = norm1(i)/nmag1 60 norm2(i) = norm2(i)/nmag2 c dot product of normal vectors: call dot(norm1,norm2,dotnorm) c sign of dihedral angle: call dot(norm1,b3,b3z) tauvalnorm = (b3z/dabs(b3z))* 1 (180.d0/3.141592653589793d0)*dacos(dotnorm) return end subroutine angle(r1,r2,r3,ang) c S/R angle calculates the angle (degrees) between vectors c r1-r2 and r3-r2 implicit real*8 (a-h,o-z) dimension r1(3),r2(3),r3(3),b1(3),b2(3) do 10 i = 1,3 b1(i) = r1(i) - r2(i) 10 b2(i) = r3(i) - r2(i) call dot(b1,b1,a1mag2) call dot(b2,b2,a2mag2) call dot(b1,b2,dotprod) dotprod = dotprod/dsqrt(a1mag2*a2mag2) ang = 180.d0*dacos(dotprod)/3.141592653589d0 return end