C PROTEIN DATA BANK SOURCE CODE MAKERING C AUTHOR. C.MARZEC,L.DAY C ENTRY DATE. 1/94 UNSUPPORTED C LAST REVISION. 1/94 C PURPOSE. CALCULATE 5-MEMBERED RING FROM C PURPOSE. PARAMETERS C LANGUAGE. FORTRAN C c Version 1.1 --- 6/25/91 --- by Christopher Marzec c c MAKERING calculates a five-membered ring structure from parameters c [q,P,s,Gamma;{bj}]; q and P are the same parameters defined by c Cremer and Pople; S and Gamma are as defined by C.J. Marzec c and L.A. Day (1993) "An exact description of five-membered ring c configurations. I. Parameterization via an amplitude S, an c angle Gamma, the pseudorotation amplitude, q, and phase angle, P, and c bond lengths." J. Biomolec. Struct. Dynam. 10, 1091-1123. c MAKERING is the inverse of BREAKRING. c c j = (1,2,3,4,5) corr. to (C1,C2,C3,C4,O4) ; this is the same as c using j = 0 for O4 in the equations. c c b(j) = |R(j+1)-R(j)| ; so b(1) = C2-C1 bond length; b(4) = O4-C4 ; c b(5) = b(0) = C1-O4 c c b(j) are the true 3-D bond lengths; bp(j) are their projections c into the CP plane. S/R make computes matrix r(i,j), which c contains x, y, z (i=1,2,3) coordinates for atom j. c c makering.f is the inverse of breakring.f, modulo solid body c translations and rotations away from "standard position", c which are projected out by breakring.f . c c Input is from keyboard or from file "breakring.output" created by breakring. c Output is on screen and into file "makering.output". c implicit real*8 (a-h,o-z) dimension b(5),r(3,5),bav(5) character*1 answer character*2 atname(5) character*7 bondname(5) c The bond length data below are empirical, averages from the Brookhaven c data sets: data bav/1.5190d0,1.5234d0,1.5226d0,1.4320d0,1.4298d0/ data bondname/'','','','',''/ data atname/'C1','C2','C3','C4','O4'/ twopi = 2.d0*3.14159265358979323d0 open(20,file='makering.output') pi = twopi/2.d0 d2r = twopi/360.d0 do 8 j = 1,5 8 b(j) = bav(j) write(6,9) 9 format('Makering accepts internal ring parameters [q,P,S,Gamma]'/, 1 'and bond lengths, and it outputs Cartesian coordinates',/ 2 'for the ring atoms.'/) c read input values: write(6,10) 10 format('Input q (A), P (deg.), S (A), and Gamma (deg)') read(5,*) qval,pval,sval,gamma write(6,21) qval,pval,sval,gamma 21 format(/'The [q,P,s,gamma] values input are [', 1 1pd10.4,2x,1pd10.4,2x,1pd10.4,2x,1pd10.4']'/) write(6,22) 22 format('Do you wish to supply your own bond lengths? Enter y or n.'/) read(5,23) answer 23 format(a1) do 12 j = 1,5 if(answer.ne.'y'.and.answer.ne.'Y') go to 15 write(6,11) bondname(j),j 11 format('Input ',a7,' bond length b',i1,':') read(5,*) b(j) 15 write(6,13) bondname(j),b(j) 13 format(a7,' bond length value is',1pd12.4) 12 continue if(sval.gt.1.d0) write(6,41) sval 41 format(/'The largest observed S value is less than 1.0 A.', 1 ' The value',1pd12.4,' may be too large;' 2 ' S/R make might not converge.') if(qval.gt.0.8d0) write(6,42) qval 42 format(/'The largest observed q value is less than 0.8 A.', 1 ' The value',1pd12.4,' may be too large,' 2 ' and S/R make might not converge.') c dmax is the maximum change in Angstroms allowed in S/R make dmax = 1. c if(sval.gt.0.8.or.qval.gt.1.0) dmax = 0.5d0 init = 1 c If init = 1, S/R sosolve initializes v values to those of a c maximally symmetrical plane-ring. If you wish to calculate a ring c from a neighboring ring (ie., without initializing the v values), c set init = 0. c S/R make requires input: [q,P,s,gamma; bonds] ; parameters dmax and c init are explained above. Output is position matrix r(i,j) c for 1=1 to 3 and j = 1 to 5. call make(qval,pval,sval,gamma,b,r,dmax,init) c Write (x,y,z) coords into output file coords.mpso write(6,1013) do 1005 j = 1,5 write(20,*) r(1,j),r(2,j),r(3,j) write(6,1004) atname(j),r(1,j),r(2,j),r(3,j) 1004 format('Atom ',a2,' has x,y,z =',3(1pd12.4)) 1005 continue write(6,1013) write(20,1013) 1013 format(/) write(20,1012) write(20,*) qval,pval,sval,gamma 1012 format('These coordinates arise from [q,p,s,Gamma] = ') write(20,1011) 1011 format('where:') do 1010 j = 1,5 1010 write(20,13) bondname(j),b(j) write(6,1015) 1015 format('The Cartesian coordinates of the constructed ring'/, 1 'are found in file "makering.output".') stop end subroutine leq(a,b,nn,mm,na,nb,err) implicit real*8(a-h,o-z) dimension a(na,nn),b(nb,mm) c leq is a linear equations subroutine to solve matrix eqn ax=b c n is dimension of a c m is number of right hand sides (columns in b) c on return, soln matrix x is in b c err = -1.d0 for singular matrix, = 1.d0 for non-singular matrix n = nn m = mm n1 = n-1 c find maximum element in each row, and divide do 19 i = 1,n r = dabs(a(i,1)) do 16 j = 2,n 16 r = dmax1(r,dabs(a(i,j))) do 17 j = 1,n 17 a(i,j) = a(i,j)/r do 18 j = 1,m 18 b(i,j) = b(i,j)/r 19 continue do 9 j = 1,n1 c find maximum element of j'th column imax = j r = dabs(a(j,j)) l = j+1 do 2 i = l,n if(r-dabs(a(i,j))) 1,2,2 1 imax = i r = dabs (a(i,j)) 2 continue if(imax-j) 6,6,3 3 do 4 k = j,n r = a(j,k) a(j,k) = a(imax,k) 4 a(imax,k) = r do 5 k = 1,m r = b(j,k) b(j,k) = b(imax,k) 5 b(imax,k) = r 6 do 9 i = l,n r = -a(i,j)/a(j,j) do 7 k = l,n 7 a(i,k) = a(i,k) + r*a(j,k) do 8 k = 1,m 8 b(i,k) = b(i,k) + r*b(j,k) 9 continue c the matrix is now in triangular form err = 1.d0 do 11 i = 1,n if(a(i,i)) 11,10,11 c tests if diagonal is zero 10 err = -1.d0 go to 15 11 continue c start the back substitution and find solution do 14 k = 1,m b(n,k) = b(n,k)/a(n,n) do 13 l = 1,n1 i = n-l r = 0.d0 imax = i+1 do 12 j = imax,n jj = i+n+1-j 12 r = r + a(i,jj)*b(jj,k) c r = r + a(i,j)*b(j,k) 13 b(i,k) = (b(i,k)-r)/a(i,i) 14 continue 15 return end subroutine make(qval,pval,sval,gamma,b,r,dmax,init) c Input values qval and sval in Angstroms; pval and gamma in degrees c c Vector v contains plane-ring vertex position variables: c v(1)= deltax(1);...v(5)=deltax(5); v(6)=deltay(1)...v(10)=deltay(5) ; c where x and y are Cartesian atom position coordinates defined in the c "standard" coordinate system. c implicit real*8 (a-h,o-z) dimension w(10,10),delta(10),f(10),bp(5),ws(10,10) dimension v(10),r(3,5),b(5) twopi = 2.d0*3.14159265358979323d0 pi = twopi/2.d0 d2r = twopi/360.d0 pvalr = pval*d2r gammar = gamma*d2r c Project the b(j) values onto the CP plane to find projected length bp(j): c use zj = q*sqrt(2/5)*cos(P-pi/2+4*pi*j/5) do 25 j = 1,5 zjp=qval*dsqrt(2.d0/5.d0)*dcos(pvalr-pi/2.d0+.4d0*twopi*(j+1)) zj=qval*dsqrt(2.d0/5.d0)*dcos(pvalr-pi/2.d0+.4d0*twopi*j) r(3,j) = zj sinaj = (zjp-zj)/b(j) bp(j) = b(j)*dsqrt(1.d0-sinaj**2) 25 continue c initialize v matrix if init = 1: if(init.ne.1) go to 11 c rp(j) is the magnitude of the radius vector to atom j. c For the initial guess, take a symmetrical plane-ring whose c sides are 1.5 A long. This guess will converge to the physically c acceptable ring structure. So rp*sin(36deg)=1.5/2., and rp = 1.27A rp = 1.27d0 do 10 j = 1,5 etaj = pi/2.d0 - 2.d0*pi*j/5.d0 c set x and y components of vector rpvec, whose length is rp v(j) = rp*dcos(etaj) v(j+5) = rp*dsin(etaj) 10 continue c initial guess for unphysical star solns (5-pointed star with O4 c at top) ******************************************* istar = 0 if(istar.eq.0) go to 1234 v(5) = 0. v(10) = .788 v(1) = .463 v(6) = -.638 v(2) = -.75 v(7) = .243 v(3) = .75 v(8) = .243 v(4) = -.463 v(9) = -.638 1234 continue c A third, "tooth" shaped initial guess generates the second family of c unphysical solutions. c **************************************************** 11 continue c We wish to solve the 10 equations for the 10 plane-ring c variables v0(j). We linearize the solution around current guess v(j), c writing v0(j) = v(j) + delta(j), c which gives 10 equations f(i) = w(i,j)*delta(j). c Begin loop to find v(1-10) iteratively in terms of s and gamma: do 100 iloop = 1,100 c first zero w and f: do 301 i = 1,10 f(i) = 0.d0 do 301 j = 1,10 w(i,j) = 0.d0 301 continue c fill w matrix: w(5,5) = 1.d0 do 30 j = 1,5 w(1,j) = 1.d0 w(2,j+5) = 1.d0 w(3,j) = 2.d0*v(j) w(3,j+5) = -2.d0*v(j+5) w(4,j) = 2.d0*v(j+5) w(4,j+5) = 2.d0*v(j) jp = j+1 if(j.eq.5) jp = 1 w(j+5,j) = -2.d0*(v(jp)-v(j)) w(j+5,jp) = 2.d0*(v(jp)-v(j)) w(j+5,j+5) = -2.d0*(v(jp+5)-v(j+5)) w(j+5,jp+5) = 2.d0*(v(jp+5)-v(j+5)) 30 continue do 31 i = 1,10 do 31 j = 1,10 31 ws(i,j) = w(i,j) c fill f vector (LHS) f(5) = -v(5) f(3) = (sval**2)*dcos(2.d0*gammar) f(4) = (sval**2)*dsin(2.d0*gammar) do 35 j = 1,5 jp = j+1 if(jp.eq.6) jp = 1 f(1) = f(1) - v(j) f(2) = f(2) - v(j+5) f(3) = f(3) - v(j)**2 + v(j+5)**2 f(4) = f(4) -2.d0*v(j)*v(j+5) f(j+5) = bp(j)**2 - ((v(jp)-v(j))**2 + (v(jp+5)-v(j+5))**2) 35 continue c store the f vector temporarily in delta: do 40 j = 1,10 40 delta(j) = f(j) c Solve for the correction vector delta: call leq(w,delta,10,1,10,10,err) c check soln by reassembling equations: do 45 i = 1,10 fcheck = 0.d0 do 44 j = 1,10 44 fcheck = fcheck + ws(i,j)*delta(j) difi = dabs(fcheck-f(i)) if(dif1.gt.1.d-8) write(6,47) i,f(i),fcheck,difi 47 format('Problem with eqn(',i2,'); f(i),fcheck,difi =', 1 3(1pd21.4)) 45 continue deltamax = 0.d0 c calculate largest delta value do 463 j = 1,10 463 if(dabs(delta(j)).gt.deltamax) deltamax = dabs(delta(j)) c adjust delta's st largest is dmax degrees = dmax*twopi/180 rad c zlam adjusts the amount of the correction vector applied. zlam = dmax/deltamax c if(zlam.gt.1.d0.or.sval.le.3.d0) zlam = 1.d0 if(zlam.gt.1.d0) zlam = 1.d0 c correct the v vector do 50 j = 1,10 50 v(j) = v(j) + delta(j)*zlam c check to see if done: c corsum is the sum of the absolute corrections. corsum = 0.d0 do 55 j = 1,10 55 corsum = corsum + dabs(delta(j)) c Uncomment next line for info about bad convergence; c if(iloop.gt.20) write(6,48) iloop,zlam,corsum 48 format('iloop, zlam, corsum =',i5,2(1pd12.4)) if(corsum.lt.1.d-5) go to 101 100 continue write(6,777) 777 format('Loop 100 of S/R make did not exit properly, 1 after convergence. Try a more moderate choice of q and/or S') 101 continue c Save the coordinates of the planering: do 110 j = 1,5 r(1,j) = v(j) r(2,j) = v(j+5) 110 continue c Now the r(i,j) contain the coordinates for the new ring, c written in the "standard" coordinate system. return end