program o2mol c c ... compile/link with: f77 -o O2MOL o2mol.f ; strip O2MOL c execute with: O2MOL < gsreal.odb c c extract display matrix from o datablock .gs_real c write molscript transformation commands c 5/93 Janet Smith c c .gs_real looks something like: c c.GS_REAL R 20 (5f13.5) c 0.22011 138.92300 35.86800 3.29500 -0.14422 c -0.11191 0.12297 0.00000 0.06680 -0.18808 c -0.09281 0.00000 0.15226 -0.02349 0.15720 c 0.00000 17.01891 23.03784 -14.44577 1.00000 c c gs(1) is sometimes the scale factor c gs(2-4) is the current screen center in Angstroms c gs(5-16) the 4x4 display matrix, i.e. rotation of the c display from stored coordinates c gs(17-19) is the scaled translation vector c c It appears that the displayed coordinates are obtained from c stored Angstrom coordinates by the following transformation: c c gs(5)*x + gs(9)*y + gs(13)*z + gs(17) = x(display) c gs(6)*x + gs(10)*y + gs(14)*z + gs(18) = y(display) c gs(7)*x + gs(11)*y + gs(15)*z + gs(19) = z(display) c c The current scale factor (zoom) is multiplied into the c elements of the transformation (gs(5-19)) c c There are two problems with the display matrix that require it c to be orthogonalized: c c 1) first element of the datablock is the scale factor at c the time of the last recentering, and may not be the scale c factor applied to the current matrix, which is updated c whenever the zoom knob is changed c c 2) beginning with o version 5.9, the default format for writing c floating point numbers (f8.3) is too imprecise for storing c elements of a rotation matrix, and users will generally not c remember to specify a format when writing out this datablock c c the matrix output by this program is a kosher orthogonal c matrix, thanks to subroutine matmgs by Morten Kjelgaard c c Eulerian angles in the Rossmann & Blow convention (modified c to transform atoms rather than axes) are also computed c c datablock for o transformation (Angstrom to Angstrom) via c .LSQ_RT_GS is computed c real gsreal(20),r(3,3) character title*26,type*2,fmt*20,line*72 c pi = acos(-1.0) rtodeg = 180.0 / pi read (5,'(a)') line write(6,*) line if (line(1:8) .ne. '.GS_REAL') then write (6,*) ' Not a .GS_REAL datablock' call exit end if read (5,*) gsreal r(1,1) = gsreal(5) r(2,1) = gsreal(6) r(3,1) = gsreal(7) r(1,2) = gsreal(9) r(2,2) = gsreal(10) r(3,2) = gsreal(11) r(1,3) = gsreal(13) r(2,3) = gsreal(14) r(3,3) = gsreal(15) x = r(1,1) + r(1,2) + r(1,3) y = r(2,1) + r(2,2) + r(2,3) z = r(3,1) + r(3,2) + r(3,3) scale = sqrt ( (x*x + y*y + z*z) / 3.0 ) ierr = 0 call matmgs (r,3,ierr) if (ierr .ne. 0) then write (6,*) ' Singular matrix' call exit end if call euler (r,th1,th2,th3) c Eulerian angles are nontranposed Rossmann/Blow c Signs must be reversed because they will be applied to c coordinates as successive rotations about z, x and z th1 = -th1 * rtodeg th2 = -th2 * rtodeg th3 = -th3 * rtodeg write (6,1) ((r(i,j),j=1,3),i=1,3),th1,th2,th3 1 format (/' For MolScript:' $ //' transform atom * by centre position in molecule _____' $ /' by rotation ',3f10.6/20x,3f10.6/20x,3f10.6,' ;' $ //' OR' $ //' transform atom * by centre position in molecule _____' $ /' by rotation z',f8.2 $ /' by rotation x',f8.2 $ /' by rotation z',f8.2,' ;') title = '.LSQ_RT_GS' number = 12 fmt = '(x,3f12.6)' tx = gsreal(17) / scale ty = gsreal(18) / scale tz = gsreal(19) / scale write (6,2) title,type,number,fmt, $ ((r(i,j),i=1,3),j=1,3),tx,ty,tz 2 format ( /' For coordinate transforms in o:' $ //2a,i10,1x,a $ /1x,3f12.6/1x,3f12.6/1x,3f12.6/1x,3f12.6) call exit end C TRANSFORM MATRIX R TO EULERIAN ANGLES C Rossmann/Blow convention, no transposition C Do all calculations in double precision C I/O single precision SUBROUTINE EULER (RR,TT1,TT2,TT3) REAL RR(3,3) DOUBLE PRECISION R(3,3),T1,T2,T3,r3 DO 5 I=1,3 DO 5 J=1,3 R(J,I) = RR(J,I) 5 CONTINUE IF (ABS(R(3,3)) .LT. 0.999999) THEN T1 = DATAN2(R(3,1),-R(3,2)) T3 = DATAN2(R(1,3), R(2,3)) IF (ABS(T3).LT.0.10 .OR. ABS(T3).GT.3.00) THEN T2 = DATAN2( R(2,3)/DCOS(T3) , R(3,3) ) ELSE T2 = DATAN2( R(1,3)/DSIN(T3) , R(3,3) ) END IF ELSE T1 = DATAN2(R(1,2),R(1,1)) T3 = 0.0 R3 = DMIN1( DABS(R(3,3)) , DBLE(1.0) ) R3 = DSIGN( R3 , R(3,3) ) T2 = DACOS(R3) END IF TT1 = T1 TT2 = T2 TT3 = T3 RETURN END c======================================================================= subroutine matmgs (A, n, errcod) * m a t m g s -- Modified Gram-Schmidt Orthogonalization by Columns. * This routine orthogonalizes a n by n square matrix using the modified * Gram-Schmidt procedure. If the input matrix has a singularity, errcod * is returned with a non-zero value. implicit none integer n, errcod real A(n,n) * When ------- Who ---------------- What ------------------------------- * 28-Mar-1990 Morten Kjeldgaard Modified some old code of mine. c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ integer i, ii, ik, ido, j real d, dot, dnorm do 50 i = 1,n-1 dnorm = 0.0 do 10 ii = 1,3 dnorm = dnorm + A(ii,i)*A(ii,i) 10 continue if (dnorm .lt. 1.0e-12) then errcod = -1 return endif ido = i + 1 do 40 j = ido, n dot = 0.0 do 20 ik = 1,n dot = dot + A(ik,i)*A(ik,j) 20 continue do 30 ik = 1,n A(ik,j) = A(ik,j) - A(ik,i) * dot/dnorm 30 continue 40 continue 50 continue * --- Normalize the n x n matrix by columns. do 80 j = 1,n d = 0.0 do 60 i = 1,n d = d + A(i,j)**2 60 continue d = sqrt (d) do 70 i = 1,n A(i,j) = A(i,j) / d 70 continue 80 continue end