Logo Search packages:      
Sourcecode: raster3d version File versions  Download package

rastep.f

      PROGRAM RASTEP
********************************************************************************
*
* Usage: 
*    rastep [-h] [-iso] [-Bcolor Bmin Bmax] [-prob xx] [-radius r] [-fancy[0-9]]
*           [-tabulate histogram.file] [-by_atomtype] [-suv_check] [-cn_check]
*
*
*     -auto       auto-orientation of viewpoint
*     -h          suppresses header records in output
*     -iso        forces isotropic B values (spheres rather than
*                 ellipsoids) even if ANISOU cards present
*     -aniso            converts isotropic B values into Uij terms so that they
*                 can be included in the analysis
*     -Bcolor Bmin Bmax 
*                 color by Biso values; Bmin = dark blue, Bmax = light red
*     -Acolor           color by anisotropy;  red < white (A=0.5) < green
*     -prob xx    draws ellipsoids to enclose this 
*                 probability level (default = 0.50)
*     -radius           draws bonds with this radius in Angstroms
*                 (default = 0.10)
*     -fancy[0-9] increasingly complex rendition of ellipsoids
*                 fancy0  [default] solid surface only
*                 fancy1  principle axes + transparent bounding ellipsoid
*                 fancy2      equatorial planes only
*                 fancy3  equatorial planes + transparent ellipsoid
*                 fancy4  longest principle axis only
*                 fancy5      for ORTEP lovers - one octant missing
*                 fancy6  same as fancy5, with missing octant colored grey
*
*===============================================================================
* The following options are used by parvati scripts
*
*     -tabulate [file]instead of creating a Raster3D input file, 
*                 list all atoms with principle axes and anisotropy.
*                 Optionally write a histogram of anisotropy to speficied 
*                 output file; otherwise output is to stderr 
*
*                       output from -tabulate for this version of rastep
*                 ATOM RESNAME RESNUM  EIGEN1 EIGEN2 EIGEN3 ANISOTROPY Uiso
*
*     -com [file] find <anisotropy> in shells from center of mass
*     -nohydrogens      don't plot hydrogens even if present
*     -mini       small size plot (176x208) with auto-orientation
*     -suv_check  use Suv to validate similarity of bonded ellipsoids
*     -cn_check   use CCuij to check similarity across C-N bonds
*
********************************************************************************
*
* EAM Jul 97      - initial version
* EAM Dec 97      - version 2.4b release
* EAM Jan 98      - add tabulation option
* EAM May 98      - integrate with PARVATI script
* EAM Jun 99      - fix bug (lack of sqrt) in -iso processing, trap read errors
*             -Acolor flag to color by anisotropy (not yet entirely satisfactory)
* EAM Jul 99      - V2.4l -tab output revised slightly, 
*             NPD ellipsoids colored magenta
*             -nohydro flag to suppress drawing hydrogens
*             -mini    flag to generate smaller pictures
*             don't draw bonds between atoms with different ALT flags
* EAM Aug 99      - V2.4m 
*             work harder at suppressing hydrogens
*             add auto-orientation (NB: scaling is wrong in this case!)
* EAM Dec 99      - V2.5
*             clean up output formats a little
* EAM Jun 2000    - additional error reporting
* EAM Jul 2000    - apply Bcolor to bonds as well as atoms
* EAM Sep 2000    - Suv similarity test
* EAM Apr 2001    - V2.6 
*             ORTEP_LIKE ellipsoids (one octant missing)
*             error count
* EAM Feb 2002    - rework ANISOU and Suv processing to gain back some speed
*             maybe I should have an array of iso/aniso flags to save
*             time during testing?
*             the rest of CARD() array can go too?
*             all the tests on IF ATOM(I)(1:).eq.'ATOM' are now unneeded
* EAM May 2003    - expand default color table to allow for off-by-one atom names
* EAM Oct 2003    - trap and report 0 values for axial Uij components in input
* EAM May 2006    - initial arrays to 0 for gfortran
* EAM Feb 2008    - 2.7s initialize even more arrays for gfortran
* EAM Aug 2009    - slightly revised output format for -tabulate
* EAM Sep 2009    - Explicitly test CCuij for C-N bonds
* EAM Oct 2009    - report CCuij = if C or N is NPD
*             -aniso flag
* EAM Mar 2010    - CCuij for phosphate backbone also
* 
*     I/O units for colour/co-ordinate input, specs output, user output
*
      INCLUDE 'VERSION.incl'
*
      INTEGER INPUT, OUTPUT, NOISE
      PARAMETER (INPUT=5, OUTPUT=6, NOISE=0)
      PARAMETER (MAXCOL=5000, MAXATM=300000)
      REAL RGB(3,MAXCOL), VDW(MAXCOL)
      REAL SPAM(5,MAXATM)
      REAL UIJ(6,MAXATM)
      real center(3)
      CHARACTER*24 MASK(MAXCOL),TEST
      CHARACTER*80 ATOM(MAXATM),CARD
      character*3  resname
      character*1  spacer
      LOGICAL MATCH
      logical           hflag, ellipses, bcflag, tflag, atflag, comflag
      logical           acflag, nohydro, mini, auto, aniflag
      integer           fancy
      character*132     flags
c
c     Data structures used for auto-orientation
      real  Rr(3,3), U(4,4), Xmom(5)
c
      COMMON /ASSCOM/ assout, verbose
      integer         assout
      logical         verbose
c
      real  quadric(10), anisou(6)
      real      eigens(4), evecs(4,4), evecinv(4,4), evecit(4,4)
      real      qq(4,4), qp(4,4), temp(4,4)
c
      external    anitoquad
      integer     anitoquad
c
      real  problevel(50)
c
      real start(3),end(3)
      real MARGIN
      parameter (MARGIN = 1.15)
c
      real  Uprin(3), Umean, Usigma, anisotropy, ellipticity
      integer     histogram(20), hislun
      real  anisi(MAXATM), sum_a, sum_a2, anis_mean, anis_sigma
      real  Biso(MAXATM), sum_b, sum_b2, sum_ab
      integer     nanis, niso, nhyd, nonpos
      logical     hwhacky
      integer     nerrors
c
      character*2 atomtype
      integer     natype
c
      real*8      wsum, xsum, ysum, zsum
      real*8      xcom, ycom, zcom
      real  adist(110)
      integer     hdist(110), comlun, dshells
c
c     Support for validation of similarity of bonded atoms
      logical     suvflag, cnflag
      integer     suvlun, suvbad, cnlun, cnbad, ccuijpair
      character*1  prevchain
      character*10 name_i, name_j
      real  anisov(6)
      real      cn_min_ccuij
c
c     Default to CPK colors and VDW radii
      integer DEFCOLS
      parameter (DEFCOLS = 17)
      character*60 defcol(DEFCOLS)
      data defcol /
     & 'COLOUR###### CA ##############   0.175   0.175   0.175  1.70',
     & 'COLOUR###### C  ##############   0.175   0.175   0.175  1.70',
     & 'COLOUR#######C################   0.625   0.625   0.625  1.70',
     & 'COLOUR#######N################   0.125   0.125   1.000  1.60',
     & 'COLOUR#######O################   0.750   0.050   0.050  1.50',
     & 'COLOUR#######S################   1.000   1.000   0.025  1.85',
     & 'COLOUR#######H################   1.000   1.000   1.000  1.20',
     & 'COLOUR#######P################   0.050   0.750   0.050  1.80',
     & 'COLOUR##### CA ###############   0.175   0.175   0.175  1.70',
     & 'COLOUR##### C  ###############   0.175   0.175   0.175  1.70',
     & 'COLOUR######C#################   0.625   0.625   0.625  1.70',
     & 'COLOUR######N#################   0.125   0.125   1.000  1.60',
     & 'COLOUR######O#################   0.750   0.050   0.050  1.50',
     & 'COLOUR######S#################   1.000   1.000   0.025  1.85',
     & 'COLOUR######H#################   1.000   1.000   1.000  1.20',
     & 'COLOUR######P#################   0.050   0.750   0.050  1.80',
     & 'COLOUR########################   1.000   0.000   1.000  2.00'
     &            /
c
c     Critical values for probability ellipsoids of a trivariate normal
c     distribution. From Table 6.1 of ORTEP-III manual (Oak Ridge National
c     Laboratory Report ORNL-6895, 1996). Tabulated below in increments of
c     2% in probability.  Default contours enclose a probability level of
c     50% (critical value 1.5382).
c
      data  problevel /     0.4299, 0.5479, 0.6334, 0.7035, 0.7644, 
     &                        0.8192, 0.8694, 0.9162, 0.9605, 1.0026,
     &                        1.0430, 1.0821, 1.1200, 1.1570, 1.1932,
     &                        1.2288, 1.2638, 1.2985, 1.3330, 1.3672,
     &                        1.4013, 1.4354, 1.4695, 1.5037, 1.5382,
     &                        1.5729, 1.6080, 1.6436, 1.6797, 1.7164,
     &                        1.7540, 1.7924, 1.8318, 1.8724, 1.9144,
     &                        1.9580, 2.0034, 2.0510, 2.1012, 2.1544,
     &                        2.2114, 2.2730, 2.3404, 2.4153, 2.5003,
     &                        2.5997, 2.7216, 2.8829, 3.1365, 6.0000 /
c
c
      assout   = noise
      verbose  = .false.
      hflag    = .false.
      acflag   = .false.
      bcflag   = .false.
      tflag    = .false.
      atflag   = .false.
      comflag  = .false.
      suvflag  = .false.
      cnflag   = .false.
      ellipses = .true.
      hwhacky  = .false.
      nohydro  = .false.
      mini     = .false.
      auto     = .false.
      aniflag  = .false.
      fancy    = 0
      prob     = 0.50
      radius   = 0.10
      nerrors  = 0
c
      prevchain = '@'
c
c     Gfortran is nuts
      wsum = 0
      xsum = 0
      ysum = 0
      zsum = 0
      sum_a = 0
      sum_b = 0
      do i=1,20
          histogram(i) = 0
      enddo
      do i=1,110
          adist(i) = 0
          hdist(i) = 0
      enddo
      do i = 1,3
          do j = 1,3
              Rr(i,j) = 0
          enddo
      enddo
      do i = 1,4
          eigens(i) = 0
          do j = 1,4
              evecs(i,j) = 0
            evecinv(i,j) = 0
            evecit(i,j) = 0
            qq(i,j) = 0
            qp(i,j) = 0
            temp(i,j) = 0
            U(i,j) = 0
          enddo
      enddo
      do i = 1,5
          Xmom(i) = 0
      enddo
c
      narg = iargc()
      i = 1
    5 continue
          call getarg( i, flags )
          if (flags(1:5) .eq. '-help') goto 701
          if (flags(1:6) .eq. '-debug') verbose = .true.
          if (flags(1:2) .eq. '-h') hflag = .true.
          if (flags(1:4) .eq. '-iso') ellipses = .false.
          if (flags(1:6) .eq. '-aniso') aniflag = .true. 
          if (flags(1:4) .eq. '-rad') then
            i = i + 1
            if (i.gt.narg) goto 701
            call getarg( i, flags )
            read (flags,*,err=701) radius
            if (radius.lt.0) radius = 0.0
          end if
          if (flags(1:4) .eq. '-pro') then
            i = i + 1
            if (i.gt.narg) goto 701
            call getarg( i, flags )
            read (flags,*,err=701) prob
            if (prob.le.0.) stop 'illegal probability level'
*           If prob > 1 assume they meant it in percent
            if (prob.gt.1.) prob = prob / 100.
          end if
          if (flags(1:5) .eq. '-Acol') then
            acflag = .true.
            bcflag = .false.
          end if
          if (flags(1:5) .eq. '-Bcol') then
            bcflag = .true.
            acflag = .false.
            i = i + 1
            if (i.gt.narg) goto 701
            call getarg( i, flags )
            read (flags,*,err=701) Bmin
            i = i + 1
            if (i.gt.narg) goto 701
            call getarg( i, flags )
            read (flags,*,err=701) Bmax
          endif
          if (flags(1:6) .eq. '-fancy') then
            if (flags(7:7).eq.'0') then
                fancy = 0
            else if (flags(7:7).eq.'1') then
                fancy = 1
            else if (flags(7:7).eq.'2') then
                fancy = 2
            else if (flags(7:7).eq.'3') then
                fancy = 3
            else if (flags(7:7).eq.'4') then
                fancy = 4
            else if (flags(7:7).eq.'5') then
                fancy = 5
            else if (flags(7:7).eq.'6') then
                fancy = 6
            else
                fancy = 1
            endif
          endif
          if (flags(1:4) .eq. '-tab') then
            tflag  = .true.
            hflag  = .true.
            acflag = .false.
            bcflag = .false.
            fancy  = 0
            do j=1,MAXATM
                anisi(j) = 0.0
            end do
            hislun = NOISE
            if (i.ge.narg) goto 799
            call getarg(i+1,flags)
            if (flags(1:1) .ne. '-') then
                hislun = 1
                open(unit=hislun,file=flags,status='UNKNOWN'
     &              )
                i = i + 1
            endif
          endif
          if (flags(1:4) .eq. '-com') then
            comflag = .true.
            comlun = NOISE
            if (i.ge.narg) goto 799
            call getarg(i+1,flags)
            if (flags(1:1) .ne. '-') then
                comlun = 2
                open(unit=comlun,file=flags,status='UNKNOWN'
     &              )
                i = i + 1
            endif
          endif
          if (flags(1:4) .eq. '-suv') then
            suvflag = .true.
            suvlun = NOISE
            suvlimit = 0.975
            if (i.ge.narg) goto 799
            call getarg(i+1,flags)
            if (flags(1:1) .ne. '-') then
                read (flags,*,err=701) suvlimit
                i = i + 1
            endif
          endif
          if (flags(1:9) .eq. '-cn_check') then
            cnflag = .true.
            cnlun = NOISE
            cn_min_ccuij = 0.95
            if (i.ge.narg) goto 799
            call getarg(i+1,flags)
            if (flags(1:1) .ne. '-') then
                cnlun = 3 
                open(unit=cnlun,file=flags,status='UNKNOWN')
                i = i + 1
            endif
          endif
          if (flags(1:8) .eq. '-by_atom') then
            atflag = .true. 
          endif
          if (flags(1:8) .eq. '-nohydro') then
            nohydro = .true. 
          endif
          if (flags(1:8) .eq. '-mini') then
            mini = .true. 
            auto = .true. 
          endif
          if (flags(1:8) .eq. '-auto') then
            auto = .true. 
          endif
c
      i = i + 1
      if (i.le.narg) goto 5
      goto 799
  701 continue
      write (noise,*) 'Raster3D Thermal Ellipsoid Program ',
     &                  VERSION
      write (noise,'(/,A)') 'syntax:'
      write (noise,'(A)')
     &      'rastep     [-h] [-iso] [-Bcolor Bmin Bmax] [-prob Plevel]'
      write (noise,'(A)')
     &  '   [-fancy[0-6]] [-radius R] [-auto]'
      write (noise,'(A,A)')
     &  '   [-nohydrogens] [-suv [suv_limit]] [-cn_check [cnfile]]'
      write (noise,'(A,A)')
     &  '   [-tabulate [tabfile]] [-by_atomtype] [-com [comfile]]'
      call exit(-1)
  799 continue

c
c Critical values for the radius corresponding to a sphere
c enclosing the requested probability level are taken from
c Table 6.1 of the ORTEP manual
      iprob = (prob+0.01)*50.
      pradius = problevel(iprob)

c
      write (noise,800)
      write (noise,*) 'Raster3D Thermal Ellipsoid Program ',
     &                  VERSION
      write (noise,*) 'E A Merritt -  11 Mar 2010'
      write (noise,800)
  800 format('************************************************')
c
      if (.not.ellipses) then
        write (noise,801) float(iprob)/50.
  801   format(' Spheres will bound Biso probability level', f5.2)
      else
        write (noise,802) float(iprob)/50.
  802   format(' Ellipsoids will bound probability level', f5.2)
      endif
      write (noise,803) pradius
  803 format(' Corresponding critical value           ', f7.4)
      if (aniflag) write (noise,*) 
     &     ' Isotropic atoms will be included in the analysis'
c
      if (acflag) then
            write (noise,*) 'Atoms will be colored based on Anisotropy'
      endif
c
      if (bcflag) then
      write (noise,*) 'Atom colors will be assigned based on Biso'
      write (noise,*) '    from dark blue = Bmin =', Bmin
      write (noise,*) '      to light red = Bmax =', Bmax
      Umin = Bmin / (8. * 3.14159*3.14159)
      Umax = Bmax / (8. * 3.14159*3.14159)
      Umin = Umin
      Umax = Umax
      endif
c
      if (.not. hflag) then
      WRITE(OUTPUT,'(A,A,I5,A)') 
     &     'Raster3D thermal ellipsoid program ',VERSION,
     &     INT(prob*100.+0.5), '% probability bounds'
      if (mini) then
        WRITE(OUTPUT,'(A)') '22  26    tiles in x,y'
      else
        WRITE(OUTPUT,'(A)') '80  64    tiles in x,y'
      endif
      WRITE(OUTPUT,'(A)') ' 8   8    pixels (x,y) per tile'
      WRITE(OUTPUT,'(A)') '4         3x3 virtual pixels -> 2x2 pixels'
      WRITE(OUTPUT,'(A)') '1 1 1     white background'
      WRITE(OUTPUT,'(A)') 'F         no, shadows are dorky'
      WRITE(OUTPUT,'(A)') '25        Phong power'
      WRITE(OUTPUT,'(A)') '0.15      secondary light contribution'
      WRITE(OUTPUT,'(A)') '0.05      ambient light contribution'
      WRITE(OUTPUT,'(A)') '0.25      specular reflection component'
      WRITE(OUTPUT,'(A)') '0.0       No perspective'
      WRITE(OUTPUT,'(A)') '1 1 1     main light source position'
      end if

c
      if (auto) then
        ASPECT = 22./26.
      else
        ASPECT = 80./64.
      endif
c
      NCOL = 0
      NATM = 0
      NANI = 0
      nanis = 0
      niso  = 0
      nhyd  = 0
      nonpos = 0
10    CONTINUE
        READ(INPUT,'(A80)',END=50) CARD
        IF (CARD(1:4).EQ.'COLO') THEN
          NCOL = NCOL + 1
          IF (NCOL.GT.MAXCOL) THEN
            WRITE(NOISE,*) 'Colour table overflow.  Increase ',
     &                     'MAXCOL and recompile.'
            STOP 10
          ENDIF
          READ(CARD,'(6X,A24,3F8.3,F6.2)') MASK(NCOL),
     &          (RGB(I,NCOL),I=1,3), VDW(NCOL)
        ELSEIF (nohydro .AND. CARD(77:78).EQ.' H') THEN
        goto 10
      ELSEIF (CARD(1:6).EQ.'ANISOU') THEN
        nani = nani + 1
        if (card(13:27).ne.atom(natm)(13:27)) goto 14
        read (card(29:70),*,err=12,end=12) (uij(i,natm),i=1,6)
        do i=1,6
            uij(i,natm) = uij(i,natm) * 0.0001
        enddo
        noerr = 1
        do i=1,3
            if (uij(i,natm).eq.0.0) then
                uij(i,natm) = 0.0001
                noerr = 0
            endif
        enddo
        if (noerr.eq.0) goto 15
        ELSEIF (CARD(1:4).EQ.'ATOM'.OR.CARD(1:4).EQ.'HETA') THEN
          NATM = NATM + 1
          IF (NATM.GT.MAXATM) THEN
            WRITE(NOISE,*) 'Atom array overflow.  Increase ',
     &                     'MAXATM and recompile.'
            STOP 20
          ENDIF
          ATOM(NATM) = CARD
        uij(1,natm) = -1.0
        ELSEIF (CARD(1:3).EQ.'END') THEN
          GO TO 50
        ENDIF
        GO TO 10
12    write(noise,*) '*** Format problem - ', card(13:70)
      nerrors = nerrors + 1
      goto 10
14    write(noise,*) '*** ANISOU record out of order - ', card(13:70)
      nerrors = nerrors + 1
      goto 10
15    write(noise,*) '*** Illegal ANISOU values - ', card(13:70)
      nerrors = nerrors + 1
      goto 10
*     Come here when EOF or 'END' record is reached
50    CONTINUE
      IF (NATM.EQ.0) THEN
        WRITE(NOISE,*) 'No atoms in input.'
        STOP 30
      ENDIF
*     Load default colors after any that were read in
      IF (NCOL.LE.MAXCOL-DEFCOLS) THEN
        DO i = 1,DEFCOLS
        NCOL = NCOL + 1
          READ(defcol(i),'(6X,A24,3F8.3,F6.2)') MASK(NCOL),
     &          (RGB(J,NCOL),J=1,3), VDW(NCOL)
        ENDDO
      ENDIF
*
      IF (NCOL.EQ.0) THEN
        WRITE(NOISE,*) 'No colours in input.'
        STOP 40
      ENDIF
      XMAX = -1E20
      XMIN =  1E20
      YMAX = -1E20
      YMIN =  1E20
      ZMAX = -1E20
      ZMIN =  1E20
      DO 100 IATM=1,NATM
        CARD = ATOM(IATM)
c
c     Do a little pre-processing to make later bookkeeping easier
c     At least screen out non-conformant PDB files that contain
c     something obviously not an element type in columns 77:78
c     Hydrogen naming conventions are totally messed up.
c     
      if (atflag) then
          resname = card(18:20)
c         if (resname.eq.'MET' .and. card(14:15).eq.'SD') then
c           atom(iatm)(77:78) = 'SD'
c         end if
          if ( resname.eq.'HOH' .or. resname.eq.'H2O'
     &          .or. resname.eq.'WAT') then
            atom(iatm)(77:78) = 'OW'
          end if
          if (atom(iatm)(78:78).ge.'0' .and. atom(iatm)(78:78).le.'9')
     &            atom(iatm)(77:78) = '  '
      end if
      if (nohydro .and. atom(iatm)(77:78).eq.'  ') then
          do 70 i = 13,16
            if (atom(iatm)(i:i).eq.' ') goto 70
            if (atom(iatm)(i:i).ge.'1' 
     &                .and. atom(iatm)(i:i).le.'4') goto 70
            if (atom(iatm)(i:i).ne.'H') goto 71
            atom(iatm)(77:78) = ' H'
   70     continue
   71     continue
      end if
          
c
        TEST = CARD(7:30)
        DO 80 ICOL=1,NCOL
          IF (MATCH(TEST,MASK(ICOL))) THEN
            READ(CARD,'(30X,3F8.3,6X,F6.2)',end=82,err=82) 
     &              X,Y,Z, Biso(IATM)
          IF (Biso(IATM).LE.0.0) THEN
            nerrors = nerrors + 1
            write(noise,*) '*** Illegal Biso ',Biso(IATM),' - ',
     &                        atom(iatm)(13:27)
            Biso(IATM) = 0.0
          ENDIF
            Uiso = Biso(IATM) / (8. * 3.14159*3.14159)
            SPAM(1,IATM) = X
            SPAM(2,IATM) = Y
            SPAM(3,IATM) = Z
          SPAM(4,IATM) = Uiso
            SPAM(5,IATM) = ICOL
          RAD  = sqrt(Uiso) * PRADIUS
            XMAX = MAX(XMAX,X+RAD)
            XMIN = MIN(XMIN,X-RAD)
            YMAX = MAX(YMAX,Y+RAD)
            YMIN = MIN(YMIN,Y-RAD)
            ZMAX = MAX(ZMAX,Z+RAD)
            ZMIN = MIN(ZMIN,Z-RAD)
c         atomtype = CARD(77:78)
c         if (atomtype.ne.'  ') then
c           weight = amass(atomtype)
c         else
            weight = 13.4
          wsum = wsum + weight
          xsum = xsum + weight * X
          ysum = ysum + weight * Y
          zsum = zsum + weight * Z
          if (uij(1,iatm).lt.0 .and. aniflag) then
            uij(1,iatm) = Uiso
            uij(2,iatm) = Uiso
            uij(3,iatm) = Uiso
            uij(4,iatm) = 0.0
            uij(5,iatm) = 0.0
            uij(6,iatm) = 0.0
c           write(0,*) 'Setting Uiso for atom ',card
          endif
            GO TO 100
          ENDIF
80      CONTINUE
        WRITE(NOISE,*) 'No colour table mask matches this atom:'
        WRITE(NOISE,*) ATOM(IATM)
        STOP 90
82    continue
      write(noise,*) 'Input format problem in record'
      write(noise,*) CARD
      STOP 90
100   CONTINUE
      XMID = (XMAX+XMIN)/2.
      YMID = (YMAX+YMIN)/2.
      ZMID = (ZMAX+ZMIN)/2.
      TX = -XMID
      TY = -YMID
      TZ = -ZMID
      xcom = xsum / wsum
      ycom = ysum / wsum
      zcom = zsum / wsum
      IF (ASPECT.GE.1.) THEN
*       The X direction is wider than the Y
        XROOM = ASPECT
        YROOM = 1.
        ZROOM = 2.
      ELSE
        XROOM = 1.
        YROOM = ASPECT
        ZROOM = 2.
      ENDIF
      XSPAN = XMAX-XMIN
      YSPAN = YMAX-YMIN
      ZSPAN = ZMAX-ZMIN
      SCALE = MAX(XSPAN/XROOM,YSPAN/YROOM,ZSPAN/ZROOM)
*     Leave a little extra room as a border:
      SCALE = SCALE / 0.90
      if (mini) then
          scale = sqrt(xspan**2+yspan**2+zspan**2) * aspect
          if (scale .lt. 9.) scale = 9.
      end if
*     
*     These are for the center-of-mass table
      DMAX  = MAX( ABS(XMAX-XCOM), ABS(XMIN-XCOM) )**2
     &      + MAX( ABS(YMAX-YCOM), ABS(YMIN-YCOM) )**2
     &      + MAX( ABS(ZMAX-ZCOM), ABS(ZMIN-ZCOM) )**2
      DMAX  = SQRT(DMAX)
      dshells = 100
*
      if (.not. hflag) then
      WRITE(OUTPUT,120) TX,TY,TZ,SCALE
120   FORMAT('1 0 0 0   input co-ordinate + radius transformation'/
     &       '0 1 0 0'/
     &       '0 0 1 0'/
     &       4F10.3)
      WRITE(OUTPUT,'(A)') '3         mixed object types'
      WRITE(OUTPUT,'(A)') '*'
      WRITE(OUTPUT,'(A)') '*'
      WRITE(OUTPUT,'(A)') '*'
      end if
c
c Auto-orientation
c 25-Aug-1999 
c     Find Eigenvectors of moment of inertia tensor.
c     Arrange smallest Eigenvalue along Y, largest along Z.
c Problems:
c     - Could emphasize side-chain over backbone by ignoring 
c       atoms O and N, but in practice this doesn't seem to help much.
c     - Scaling is wrong, because it was done before rotation.
c
      if (auto) then
            do 125 iatm = 1, natm
        if   (atom(iatm)(1:4).ne.'ATOM' 
     &        .and. atom(iatm)(1:4).ne.'HETA') goto 125
C         if (atom(iatm)(13:15).eq.' O ')  goto 125
C         if (atom(iatm)(13:15).eq.' N ')  goto 125
        if (atom(iatm)(77:78).eq.' H' .and. nohydro) goto 125
        x = spam(1,iatm) - xcom
        y = spam(2,iatm) - ycom
        z = spam(3,iatm) - zcom
        Rq = (x*x + y*y + z*z)
        Rv = sqrt(Rq)
        Rr(1,1) = Rr(1,1) + x*x
        Rr(1,2) = Rr(1,2) + x*y
        Rr(1,3) = Rr(1,3) + x*z
        Rr(2,1) = Rr(2,1) + y*x
        Rr(2,2) = Rr(2,2) + y*y
        Rr(2,3) = Rr(2,3) + y*z
        Rr(3,1) = Rr(3,1) + z*x
        Rr(3,2) = Rr(3,2) + z*y
        Rr(3,3) = Rr(3,3) + z*z
        Xmom(1) = Xmom(1) + 1.
        Xmom(2) = Xmom(2) + Rv
        Xmom(3) = Xmom(3) + Rq
        Xmom(4) = Xmom(4) + Rv*Rq
        Xmom(5) = Xmom(5) + Rq*Rq
  125 continue
      if (verbose) then
        write (NOISE,'(/,A,5G13.5)') ' Radial moments:',Xmom
        write (NOISE,'(A)')        ' Moment of inertia tensor'
      end if
      Rg = sqrt(Xmom(3)/Xmom(1))
      U(1,1) = Xmom(3)
      U(2,2) = Xmom(3)
      U(3,3) = Xmom(3)
      do k = 1,3
        do l = 1,3
          U(k,l) = (U(k,l) - Rr(k,l)) / Xmom(1)
        end do
        if (verbose) write (NOISE,'(3G13.6)') (U(k,l),l=1,3)
      end do
      call jacobi( U, 3, 4, Eigens, Evecs )
c     Re-order so that long axis is vertical, short axis towards viewer
        kmax = 1
        if (Eigens(2).gt.Eigens(kmax)) kmax = 2
        if (Eigens(3).gt.Eigens(kmax)) kmax = 3
        kmin = 1
        if (Eigens(2).le.Eigens(kmin)) kmin = 2
        if (Eigens(3).le.Eigens(kmin)) kmin = 3
        kmid = 6 - (kmin + kmax)
      if (verbose) then
        write (NOISE,'(A,/,3F13.5,10X,3i2)') ' Eigenvalues:', 
     &            Eigens(kmin),Eigens(kmid),Eigens(kmax)
     &            ,kmin,kmid,kmax
        write (NOISE,'(A)') ' Eigenvectors:'
        do k = 1,3
          write (NOISE,'(3F13.5)') 
     &            Evecs(k,kmin),Evecs(k,kmid),Evecs(k,kmax)
        enddo
      end if
      do k = 1,3
        U(k,1) = Evecs(k,kmid)
        U(k,2) = Evecs(k,kmin)
        U(k,3) = Evecs(k,kmax)
        U(k,4) = 0.0
        U(4,k) = 0.0
      enddo
      U(4,4) = 1.0
c     Beware! may be left-handed at this point!
      det = tinv4( evecinv, U )
      if (det .lt. 0.0) then
        do k = 1,3
          U(k,3) = -U(k,3)
        enddo
      endif
c     OK, now it should be right-handed
      WRITE(OUTPUT,'(A)') '# Auto-orientation matrix'
      WRITE(OUTPUT,'(A)') '16'
      WRITE(OUTPUT,'(A)') 'ROTATION'
      do k = 1,3
        write (OUTPUT,'(3F13.5)') U(1,k),U(2,k),U(3,k)
      enddo
      WRITE(OUTPUT,'(A)') '# End auto-orientation'
      end if
c
c Label output records
c
      if (.not. tflag) then
      WRITE(OUTPUT,'(A,A)') 
     &            '# Thermal ellipsoids from Rastep Version ',VERSION
      WRITE(OUTPUT,'(A,F5.2)') '# Probability level',float(iprob)/50.
      end if
c
c Write ellipsoids to input file for render
c
      IF (fancy.eq.0 .and. .not.tflag) GOTO 139
c
c First, optional pass, to write fancy stuff associated with ellipsoids
c
      IATM = 1
  130 CONTINUE
      IF (ATOM(IATM)(1:4).EQ.'ATOM' .OR.
     &    ATOM(IATM)(1:4).EQ.'HETA') THEN
      X = SPAM(1,IATM)
      Y = SPAM(2,IATM)
      Z = SPAM(3,IATM)
      ICOL = SPAM(5,IATM)
      if (bcflag) then
          call U2RGB( SPAM(4,IATM), Umin, Umax, RED, GREEN, BLUE )
          RED   = RED*RED
          GREEN = GREEN*GREEN
          BLUE  = BLUE*BLUE
      else if (acflag) then
          call A2RGB( 1.0, RED, GREEN, BLUE )
          RED   = RED*RED
          GREEN = GREEN*GREEN
          BLUE  = BLUE*BLUE
      else 
          RED   = RGB(1,ICOL)
          GREEN = RGB(2,ICOL)
          BLUE  = RGB(3,ICOL)
      endif
      IF (ellipses .and. uij(1,iatm).ge.0.) THEN
          do i=1,6
            anisou(i) = uij(i,iatm)
          enddo
          if (anitoquad(anisou,pradius,quadric,eigens,evecs).lt.0)then
              write(noise,*) '*** Non-positive definite ellipsoid - ',
     &                        atom(iatm)(13:27)
            nonpos = nonpos + 1
            nerrors = nerrors + 1
            Biso(iatm) = 0.0
            goto 138
          endif
          goto 132
  132     continue
          radlim = pradius * max( eigens(1),eigens(2),eigens(3) )
          radlim = radlim * MARGIN
c
c     Only for debugging ellipsoids
          if (verbose) then
            write (noise,901) 'ANISOU ',X,Y,Z,ANISOU
            write (noise,902) 'QUADRIC',QUADRIC
            write (noise,903) 'Eigenvalues', (EIGENS(i),i=1,3),
     &                 'prob', prob,'limiting radius', radlim
            write (noise,904) 'Evecs ',((evecs(i,j),i=1,3),j=1,3)
          endif
901       format(a,3f8.3,6f8.4)
902       format(a,10f8.3)
903       format(a,3f8.3,4x,a,f8.3,4x,a,f8.3)
904       format(a,9f7.3)
c
c     Tabulate principal axes of ellipsoid for each atom
        if (tflag) then
          do i=1,3
            Uprin(i) = eigens(i)**2
          enddo
          if (Uprin(2).gt.Uprin(1)) then
            Umean    = Uprin(1)
            Uprin(1) = Uprin(2)
            Uprin(2) = Umean
          endif
          if (Uprin(3).gt.Uprin(1)) then
            Umean    = Uprin(1)
            Uprin(1) = Uprin(3)
            Uprin(3) = Umean
          endif
          if (Uprin(3).gt.Uprin(2)) then
            Umean    = Uprin(2)
            Uprin(2) = Uprin(3)
            Uprin(3) = Umean
          endif
c
c       Anisotropy we define as Umin / Umax
c       as in shelxpro output
          anisotropy  = min(Uprin(1),Uprin(2),Uprin(3))
     &                  / max(Uprin(1),Uprin(2),Uprin(3))
c
c       But don't count atoms which are perfectly isotropic
          if (atom(iatm)(77:78).eq.' H') then
            nhyd  = nhyd + 1
            if (anisotropy .ne. 1.0) hwhacky = .true.
          else if (anisotropy .eq. 1.0) then
            niso  = niso + 1
          else
            anisi(iatm) = anisotropy
            sum_A = sum_A + anisotropy
            sum_B = sum_B + Biso(iatm)
            nanis = nanis + 1
          end if
c
c       Ellipticity we define as 1 / anisotropy
          if (anisotropy.eq.0) then
            ellipticity = 0
          else
            ellipticity = 1. / anisotropy
          end if
c
c       Longhi et al (1997) JMB 268, 779-799.
c       proposed another measure A = sigU / meanU
          Umean  = (Uprin(1) + Uprin(2) + Uprin(3)) / 3.0
          Usigma = (Uprin(1)-Umean)**2 
     &               + (Uprin(2)-Umean)**2 + (Uprin(3)-Umean)**2
          Usigma = sqrt( Usigma ) / 3.0
          alonghi = Usigma / Umean
c
c       Might want to check correlation with Uiso
          Uiso = SPAM(4,iatm)
c
c       Cosmetic changes to atom identifier for the sake of sorting
c       We will force there to be exactly three entities printed.
c       PDB format is just a mess:
c       cols 13:16      atom
c       col     17      alternate conf
c       cols 18:20      residue
c       col     22      chain
c       cols 23:27      resnum
c
          do i = 16, 13, -1
            if (ATOM(iatm)(i:i) .ne. ' ') j = i
          enddo
          do i = 13, 17
            if (ATOM(iatm)(i:i) .ne. ' ') k = i
          enddo
          do i = j, k
            if (ATOM(iatm)(i:i) .eq. ' ') ATOM(iatm)(i:i) = '_'
          enddo
c         if (ATOM(iatm)(17:17) .ne. ' ') then
c           do i = 18,19
c           if (ATOM(iatm)(i:i) .eq. ' ') ATOM(iatm)(i:i) = '_'
c           enddo
c         endif
          spacer = ' '
          if (ATOM(iatm)(22:22) .ne. ' ') then
            spacer = '_'
            do i = 23,25
            if (ATOM(iatm)(i:i) .eq. ' ') ATOM(iatm)(i:i) = '_'
            enddo
          endif
          write (output,905) ATOM(iatm)(13:17),
     &            ATOM(iatm)(18:22),spacer,ATOM(iatm)(23:27),
     &            Uprin(1),Uprin(2),Uprin(3),anisotropy,Uiso
905     format(A5,1X,A5,A1,A5,3F9.4,2X,F9.4,F9.4,F9.4)
c
c       Also make a histogram of anisotropies
          i = anisotropy * 20. + 1
          histogram(i) = histogram(i) + 1
c
c       And a table of <anis> by distance from center of mass
          if (comflag .and. anisotropy.lt.1.0) then
            dist = sqrt( (SPAM(1,iatm)-XCOM)**2 
     &                   + (SPAM(2,iatm)-YCOM)**2 
     &                       + (SPAM(3,iatm)-ZCOM)**2 )
            i = (dist/dmax) * float(dshells) + 1
            adist(i) = adist(i) + anisotropy
            hdist(i) = hdist(i) + 1
          endif
c
c       End tabulation code
        endif

c
c     Skip hydrogens if requested
        if (nohydro .and. ATOM(IATM)(77:78).eq.' H') goto 138

c
c     Draw principal axes inside bounding ellipsoid
        if (fancy.eq.1) then
          do i=1,3
            size = eigens(i) * pradius - 0.02
            start(1) = x - size*evecs(1,i)
            start(2) = y - size*evecs(2,i)
            start(3) = z - size*evecs(3,i)
            end(1)   = x + size*evecs(1,i)
            end(2)   = y + size*evecs(2,i)
            end(3)   = z + size*evecs(3,i)
            write (output,907) start, end
907         format(' 3',/,
     &             3f9.3,' 0.02',3f9.3,' 0.02','  0.5 1.0 0.3')
          enddo
        endif
c
c     Draw longest principle axis only
c     (experimental use only - not supported or documented)
        if (fancy.eq.4) then
          imax = 1
          if (eigens(2).gt.eigens(imax)) imax = 2
          if (eigens(3).gt.eigens(imax)) imax = 3
          size = eigens(imax) * pradius 
          imin = 1
          if (eigens(2).lt.eigens(imin)) imin = 2
          if (eigens(3).lt.eigens(imin)) imin = 3
          imed = 6 - (imax+imin)
          size = size * eigens(imax)/eigens(imed)
          start(1) = x - size*evecs(1,imax)
          start(2) = y - size*evecs(2,imax)
          start(3) = z - size*evecs(3,imax)
          end(1)   = x + size*evecs(1,imax)
          end(2)   = y + size*evecs(2,imax)
          end(3)   = z + size*evecs(3,imax)
          write (output,907) start, end
        endif
c
c     Construct 3 quadrics corresponding to the 3 orthogonal planes
c     through the center of our ellipsoid
      if (fancy.eq.2 .or. fancy.eq.3) then
        eigens(4) = 1.0
        evecs(4,4)= 1.0
        det = tinv4( evecinv, evecs )
        call trnsp4( evecit, evecinv )
        do k = 1,3
          do i = 1,4
          do j = 1,4
            QQ(i,j) = 0.0
          enddo
          QQ(i,i) = 1. / (eigens(i)*eigens(i))
          enddo
        QQ(k,k) = 1000.
        QQ(4,4) = -pradius*pradius
        call tmul4( TEMP, QQ, evecinv )
        call tmul4( QP, evecit, TEMP )
        if (acflag) then
          call A2RGB( anisotropy, red, green, blue )
          red   = red*red
          green = green*green
          blue  = blue*blue
        endif
        write (output,151) 14, X,Y,Z, radlim, red, green, blue
        write (output,152) QP(1,1),QP(2,2),QP(3,3),QP(1,2),QP(2,3),
     &                       QP(1,3),QP(1,4),QP(2,4),QP(3,4),QP(4,4)
        enddo
      endif
c
      endif
      ENDIF
  138 continue
      IATM = IATM + 1
      IF (IATM.LE.NATM) GOTO 130
c
c     Set transparency for enclosing ellipoids
      if (fancy.eq.1 .or. fancy.eq.3) then
      write (output,'(A,/,A)') '9 Begin transparent ellipsoids','8 '
      write (output,'(A)') ' 15.  0.6   1.0 1.0 1.0   0.6   0 0 0 0'
      else if (fancy.eq.2 .or. fancy.eq.4) then
      goto 160
      endif
c
  139 CONTINUE
c
c If we're just tabulating ellipticity, start making tables
c
      if (tflag) then
      total = 0
      do i = 1,20
          total = total + histogram(i)
      end do
      if (total.eq.0) then
          write (noise,*)  'No ANISOU records found'
          call exit(-1)
      end if
      if (nanis.eq.0) then
         write (noise,*)   'No anisotropic atoms found'
         call exit(-1)
      end if
      if (hwhacky) then
         write(noise,*)
     &                       'You seem to have anisotropic hydrogens',
     &                   ' - is this some kind of joke?'
         nerrors = nerrors + 1
      end if

      write (hislun,'(A)') '# Anisotropy  Fraction   Number'
      write (hislun,'(A)') '#   range     of atoms of atoms'
      do i = 1,20
          write (hislun,'(2F5.2,3X,F8.3,I10)') 
     &            (float(i)-1.)/20., float(i)/20., 
     &            float(histogram(i))/total, histogram(i)
      end do
c       Calculate mean and sigma of distribution
      sum_a2 = 0.0
      sum_b2 = 0.0
      sum_ab = 0.0
      anis_mean  = sum_A / float(nanis)
      anis_sigma = 0.0
      biso_mean  = sum_B / float(nanis)
      ccoef      = 0.0
      do i = 1,NATM
          if (anisi(i).ne.0) then
            sum_a2 = sum_a2 + (anisi(i) - anis_mean)**2
            sum_b2 = sum_b2 + (Biso(i)  - biso_mean)**2
            sum_ab = sum_ab 
     &                   + (anisi(i)-anis_mean) * (Biso(i)-biso_mean)
          end if
      end do
      if (nanis.gt.1) then
          anis_sigma = sqrt( sum_a2 / float(nanis-1) )
          ccoef = sum_ab / sqrt(sum_a2 * sum_b2)
      end if

      write (hislun,'(A)')    '#'
      write (hislun,'(A,I10)')'#  number of ANISOU records:',nani
      write (hislun,'(A,I10)')'#       non-isotropic atoms:',nanis
      write (hislun,'(A,I10)')'#           isotropic atoms:',niso
      if (nonpos.gt.0)
     &  write (hislun,'(A,I10)')'# nonpositive-definite APDs:',nonpos
      if (nhyd.gt.0)
     &  write (hislun,'(A,I10)')'#          ANISOU hydrogens:',nhyd
      write (hislun,'(A)') '#'
      write (hislun,'(A,F7.3)')
     &       '# correlation of anisotropy with B_iso:', ccoef
      write (hislun,'(A)') '#'
      write (hislun,'(A)') '#              Anisotropy  B_iso'
      write (hislun,'(A)') '#  AtomType   mean  sigma   mean  number'
      write (hislun,'(A)') '# ---------   ----------- ------  ------'
      write (hislun,'(A,F7.3,F7.3,F7.2,I8)')    
     &       '#|    Total',anis_mean,anis_sigma,biso_mean,nanis

c
c     If we're tabulating ellipticity, but not by atom_type, then we're done
c
      if (.not. atflag) goto 145
  140 continue
c
C     Find an atom type we haven't done yet, exit if none left
c     This only works if the PDB file contains the atom type in
c     columns 77:78 (standard since sometime in 1997, but many
c     files do not conform to this standard)
c
      do i = 1, NATM
          if (anisi(i).ne.0.0 .and. atom(i)(77:78).ne.'  ') then
            atomtype = atom(i)(77:78)
            goto 141
          end if
      end do
      goto 145
  141 continue

      sum_A     = 0.0
      sum_B     = 0.0
      sum_a2    = 0.0
      biso_mean = 0.0
      natype    = 0
c
c     Loop over atoms looking for the right ones
c
      do i = 1, NATM
          if (anisi(i).ne.0.0 .and. atomtype.eq.atom(i)(77:78)) then
            natype = natype + 1
            sum_A = sum_A + anisi(i)
          end if
      end do
      anis_mean  = sum_A / float(natype)
      anis_sigma = 0.0
      do i = 1, NATM
          if (anisi(i).ne.0.0 .and. atomtype.eq.atom(i)(77:78)) then
            sum_a2 = sum_a2 + (anisi(i)-anis_mean)**2
            biso_mean = biso_mean + Biso(i)
            atom(i)(77:78) = '  '
          end if
      end do
      if (natype.gt.1) anis_sigma = sqrt( sum_a2 / float(natype-1) )
      biso_mean  = biso_mean / float(natype)

      write (hislun,'(A,4X,A,F7.3,F7.3,F7.2,I8)')    
     &       '#|   ',atomtype,anis_mean,anis_sigma,biso_mean,natype
      goto 140
      end if
      IATM = 1
      goto 150
c
c Write out plot of mean anisotropy as a function of distance from c.o.m.
c The hisclean routine is not strictly necessary; it collapses shells at
c the extreme ranges of distance that contain only a few atoms.
c
  145 continue
      if (comflag) then
            write (comlun,'(A/A,A/A/A)') '#',
     &          '# Mean anisotropy as a function of distance',
     &      ' from center of mass',
     &          '#',
     &          '#  Distance   <anis>   #atoms'
      call hisclean( adist, hdist, dshells, nanis )
      do i = 1, dshells
          if (hdist(i).gt.0) then
            adist(i) = adist(i) / float(hdist(i))
            dmid = (dmax/float(dshells))*(float(i)-0.5)
            write (comlun,146) dmid, adist(i), hdist(i)
          endif
      enddo
  146 format(F10.3,F10.3,I10)
      endif
      if (suvflag.or.cnflag) goto 160
      call exit(0)

c
c Second pass write a single sphere or ellipsoid for each atom
c
  150 CONTINUE
      IF (ATOM(IATM)(1:4).EQ.'ATOM' .OR.
     &    ATOM(IATM)(1:4).EQ.'HETA') THEN
      if (nohydro .and. ATOM(IATM)(77:78).eq.' H') goto 154
      X = SPAM(1,IATM)
      Y = SPAM(2,IATM)
      Z = SPAM(3,IATM)
      ICOL = SPAM(5,IATM)
      if (bcflag) then
          call U2RGB( SPAM(4,IATM), Umin, Umax, RED, GREEN, BLUE )
          RED   = RED*RED
          GREEN = GREEN*GREEN
          BLUE  = BLUE*BLUE
      else if (acflag) then
          call A2RGB( 1.0, RED, GREEN, BLUE )
          RED   = RED*RED
          GREEN = GREEN*GREEN
          BLUE  = BLUE*BLUE
      else
          RED   = RGB(1,ICOL)
          GREEN = RGB(2,ICOL)
          BLUE  = RGB(3,ICOL)
      endif
      IF (.not. ellipses) THEN
          RAD  = sqrt(SPAM(4,IATM)) * PRADIUS
            WRITE(OUTPUT,151) 2, X,Y,Z,RAD,RED,GREEN,BLUE
      ELSE IF (uij(1,iatm).gt.0.) THEN
          do i=1,6
            anisou(i) = uij(i,iatm)
          enddo
          if (anitoquad(anisou,pradius,quadric,eigens,evecs).lt.0)then
              write(noise,*) '*** Non-positive definite ellipsoid - ',
     &                        atom(iatm)(13:26)
            nerrors = nerrors + 1
            Biso(iatm) = 0.0
            red   = 1.0
            green = 0.0
            blue  = 1.0
          endif
          radlim = pradius * max( eigens(1),eigens(2),eigens(3) )
          radlim = radlim * MARGIN
          if (acflag) then
            anisotropy = min( eigens(1),eigens(2),eigens(3) )
     &                       / max( eigens(1),eigens(2),eigens(3) )
            anisotropy = anisotropy * anisotropy
            call A2RGB( anisotropy, red, green, blue )
            red   = red*red
            green = green*green
            blue  = blue*blue
          endif
          if (fancy.eq.5 .or. fancy.eq.6) then
            write (output,'(I2,/,A,I2,/,A)') 8,
     &                ' -1.0 -1.0  -1.0 -1.0 -1.0  0.0  0 0 0',
     &                  fancy-1, 'ORTEP_LIKE'
            if (fancy.eq.6) write (output,171) 0.5, 0.5, 0.5
            write (output,172) 0, x,y,z, (evecs(i,1),i=1,3)
            write (output,172) 0, x,y,z, (evecs(i,2),i=1,3)
            write (output,172) 0, x,y,z, (evecs(i,3),i=1,3)
  171       format('BOUNDING_COLOR ',3F6.3)
  172       format('BOUNDING_PLANE ',I2,6F10.4)
          endif
          write (output,151) 14, x,y,z,radlim,red,green,blue
          write (output,152) (quadric(i),i=1,10)
      ELSE
          RAD  = sqrt(SPAM(4,IATM)) * PRADIUS
            WRITE(OUTPUT,151) 2, X,Y,Z,RAD,RED,GREEN,BLUE
  151 FORMAT(I2,/,3(1X,F8.3),4F8.3)
  152 FORMAT(10F12.4)
      ENDIF
      goto 154
      ENDIF
  154 continue
      IATM = IATM + 1
      IF (IATM.LE.NATM) GOTO 150
      IF (fancy.eq.1 .or. fancy.eq.3) then
      write (output,'(A)') '9 end transparent ellipsoids'
      endif
      IF (fancy.eq.5 .or. fancy.eq.6) then
      write (output,'(A)') '9 end ortep ellipsoids'
      endif
  160 continue
c
c Write bonds to file also. Atoms are considered bonded if they lie
c closer to each other than 0.6 * sum of VDW radii.
C If two atoms of different colors are bonded, make half-bond
C cylinders with each color.
C
      if (nerrors.eq.0) write(noise,*) 
     &       '... no errors found in input file'
      if (radius.eq.0.0 .and. .not.suvflag .and. .not.cnflag) goto 210
      if (suvflag) then
         write (suvlun,'(A,A,F5.3,A)') 
     &      'Checking for neighboring atoms with dissimilar Uij ',
     &          '(Suv < ',suvlimit,')...'
         suvbad = 0
      endif
      if (cnflag) then
         write (NOISE,'(A,F5.3)') 
     &      'Checking for C-N linkages with CCuij < ', cn_min_ccuij
         cnbad = 0
      endif
c
      DO 202 IATM=1,NATM
      IF (nohydro .AND. ATOM(IATM)(77:78).EQ.' H') GOTO 202
      if (suvflag .or. cnflag) then
        if (Biso(IATM).eq.0.0) goto 202
        if (uij(1,IATM).lt.0.) goto 202
      endif
      XI = SPAM(1,IATM)
      YI = SPAM(2,IATM)
      ZI = SPAM(3,IATM)
      ICOL = SPAM(5,IATM)
      VDWI = VDW(ICOL)
      DO 201 JATM=IATM+1,NATM
      CLOSE2 = 4.537
      DX2 = (XI - SPAM(1,JATM))**2
      if (dx2 .gt. close2) goto 201
      DY2 = (YI - SPAM(2,JATM))**2
      if (dy2 .gt. close2) goto 201
      DZ2 = (ZI - SPAM(3,JATM))**2
      DIST2 = DX2 + DY2 + DZ2
c
c       Checking for bonded atoms with dissimilar Uij
        if (suvflag .or. cnflag) then
          IF (DIST2 .GT. CLOSE2) GOTO 201
          IF (Biso(JATM).eq.0.0) goto 180
          if (uij(1,jatm).lt.0.) goto 180
CDEBUG Version 2.6c used explicit test on 0.6*VdW distance
CDEBUG but this is very slow so I have fixed CLOSE2 at the C-S bond length
CDEBUG It might be worth doing a first cut using, say, 2.25A (>S-S bond)
CDEBUG and then only check the actual VdW distance for pairs making the cut
C         JCOL = SPAM(5,JATM)
C         VDWJ = VDW(JCOL)
C         CLOSE2 = (0.6 * (VDWI + VDWJ)) **2
C         IF (DIST2 .GT. CLOSE2) GOTO 201
CDEBUG Add this section back to see if results match V2.6c numbers
          similarity = Suv( uij(1,iatm), uij(1,jatm) )
          if (suvflag .and. (similarity .lt. suvlimit)) then
            write (suvlun,161) 
     &                ATOM(IATM)(13:17),ATOM(IATM)(18:27),
     &                  ATOM(JATM)(13:17),ATOM(JATM)(18:27),
     &                similarity
            suvbad = suvbad + 1
  161     format(1X,A5,1X,A10,8X,A5,1X,A10,4X,'Suv = ',F8.4)
          endif
c
c         Check in particular for C-N links that may be TLS group boundaries
  180     continue
          if (cnflag) then
              ccuijpair = 0
              if ((ATOM(IATM)(13:15) .eq. ' C ') .and. (ATOM(JATM)(13:15) .eq. ' N ')) then
                  ccuijpair = 1
            endif
              if ((ATOM(IATM)(13:15) .eq. ' O3') .and. (ATOM(JATM)(13:15) .eq. ' P ')) then
                  ccuijpair = 2
            endif

              if (ccuijpair.ne.0) then
                  cc = CCuv( uij(1,iatm), uij(1,jatm) )
                if (ATOM(IATM)(22:22).ne.prevchain) then
                      prevchain = ATOM(IATM)(22:22)
                    write (cnlun,'(1x)')
                    write (cnlun,'(1x)')
                endif
                name_i = ATOM(IATM)(18:27)
                name_j = ATOM(JATM)(18:27)
                if (name_i(5:5) .eq. ' ') name_i(5:5) = 'X'
                if (name_j(5:5) .eq. ' ') name_j(5:5) = 'X'
                do ichar = 6,9
                    if (name_i(ichar:ichar).eq.'_') name_i(ichar:ichar) = ' '
                    if (name_j(ichar:ichar).eq.'_') name_j(ichar:ichar) = ' '
                enddo
                write (cnlun,182) name_i, name_j, cc
                if (cc .lt. cn_min_ccuij) then
                    if (ccuijpair.eq.1) write (NOISE,183) name_i, name_j, cc
                    if (ccuijpair.eq.2) write (NOISE,184) name_i, name_j, cc
                    cnbad = cnbad + 1
                endif
  182         format(1X,A10,8X,A10,4X,'CCuij = ',F8.4)
  183         format(1X,A15,8X,A15,4X,'C-N CCuij = ',F8.4)
  184         format(1X,A15,8X,A15,4X,'O3''-P CCuij = ',F8.4)
              endif
          endif

          if (tflag) goto 201
        endif

c     More stringent distance test when drawing bonds
        IF (nohydro .AND. ATOM(JATM)(77:78).EQ.' H') GOTO 201
        JCOL = SPAM(5,JATM)
        VDWJ = VDW(JCOL)
        CLOSE2 = (0.6 * (VDWI + VDWJ)) **2
        IF (DIST2  .GT. CLOSE2) GOTO 201

c       Don't draw bonds between alternate conformers of same residue
        IF (ATOM(IATM)(17:17).ne.' '.and.ATOM(JATM)(17:17).ne.' '
     &        .and. ATOM(IATM)(17:17).ne.ATOM(JATM)(17:17)) goto 201
c
c       Atoms coloured by B value
        if (bcflag) then
           write(output,211)
     1            SPAM(1,IATM),SPAM(2,IATM),SPAM(3,IATM),radius,
     2            SPAM(1,JATM),SPAM(2,JATM),SPAM(3,JATM),radius,
     3            1.0, 1.0, 1.0
           call U2RGB( SPAM(4,IATM), Umin, Umax, RED1, GREEN1, BLUE1 )
           call U2RGB( SPAM(4,JATM), Umin, Umax, RED2, GREEN2, BLUE2 )
           write(output,212)
     1            RED1,GREEN1,BLUE1, RED2,GREEN2,BLUE2, 0., 0., 0.
c       Same color atoms
        elseif (RGB(1,ICOL) .EQ. RGB(1,JCOL) .AND.
     1       RGB(2,ICOL) .EQ. RGB(2,JCOL) .AND.
     2       RGB(3,ICOL) .EQ. RGB(3,JCOL)) THEN
          WRITE(OUTPUT,211)
     1         SPAM(1,IATM),SPAM(2,IATM),SPAM(3,IATM),radius,
     2         SPAM(1,JATM),SPAM(2,JATM),SPAM(3,JATM),radius,
     3         RGB(1,ICOL),RGB(2,ICOL),RGB(3,ICOL)
        ELSE
          DO K=1,3
            center(K) = (SPAM(K,IATM)+SPAM(K,JATM))/2
          ENDDO
          WRITE(OUTPUT,211)
     1         SPAM(1,IATM),SPAM(2,IATM),SPAM(3,IATM),radius,
     2         center(1),center(2),center(3),radius,
     3         RGB(1,ICOL),RGB(2,ICOL),RGB(3,ICOL)
          WRITE(OUTPUT,211)
     1         center(1),center(2),center(3),radius,
     2         SPAM(1,JATM),SPAM(2,JATM),SPAM(3,JATM),radius,
     3         RGB(1,JCOL),RGB(2,JCOL),RGB(3,JCOL)
        ENDIF
  201 CONTINUE
  202 CONTINUE
  210 CONTINUE

C211  FORMAT(1H3,/,11(f8.3))
211   FORMAT(1H3,/,3(1x,f8.3),f7.3,3(1x,f8.3),f7.3,3(f6.3))
212   FORMAT(2H17,/,9f8.3)
c
      if (suvflag) then
          if (suvbad.eq.0) write (suvlun,*) '... No Suv outliers'
      endif
      if (cnflag) then
          if (cnbad.eq.0) write (suvlun,*) '... No CCuij outliers'
      endif
      if (suvflag .or. cnflag) call exit(0)
c
      write (noise,'(/)')
      write (noise,156) 'X  min max center-of-mass:', XMIN, XMAX, xcom
      write (noise,156) 'Y  min max center-of-mass:', YMIN, YMAX, ycom
      write (noise,156) 'Z  min max center-of-mass:', ZMIN, ZMAX, zcom
      write (noise,156) '     scale:', SCALE
  156 format(1x,a,2f8.2,f10.3)
      END


      LOGICAL FUNCTION MATCH (SUBJ, MASK)
      CHARACTER*24 SUBJ,MASK
      MATCH = .FALSE.
      DO 10 I = 1, 24
        IF (SUBJ(I:I).NE.MASK(I:I) .AND. MASK(I:I).NE.'#') RETURN
10    CONTINUE
      MATCH = .TRUE.
      RETURN
      END

C     Dummy routine to make linker happy (called by QINP in quadric.f)
      SUBROUTINE TRANSF (X,Y,Z)
      RETURN
      END

      SUBROUTINE ASSERT (LOGIC, DAMMIT)
      LOGICAL LOGIC
      CHARACTER*(*) DAMMIT
      COMMON /ASSOUT/ NOISE
      IF (LOGIC) RETURN
      WRITE (NOISE,*) 'ERROR >>>>>> ',DAMMIT
C     STOP 1234
      CALL EXIT(-1)
      END


CCC   Return RGB triple that runs from dark blue at Bmin 
CC    to light red at Bmax
C
      subroutine U2rgb( Uiso, Umin, Umax, r, g, b )
      real              Uiso, Umin, Umax, r, g, b
c
      real    fraction, h, s, v
c
      fraction = (Uiso-Umin) / (Umax-Umin)
      if (fraction.lt.0.) fraction = 0.
      if (fraction.gt.1.) fraction = 1.
      h = 240. * (1.-fraction)
      s = 0.95
      v = 0.75 + fraction/4.
      call hsv2rgb( h, s, v, r, g, b )
      return
      end


CCC   Return RGB triple that runs from 
CC    red for A=0.0 -> white for A=0.5 -> green for A=1.0
C
      subroutine A2rgb( A, r, g, b )
      real              A, r, g, b
c
      real    fraction, h, s, v
c
      if (A .lt. 0.5) h = 0.0
      if (A .ge. 0.5) h = 120.
      fraction = abs( (A-0.5)*2.0 )
c     s = fraction**2
      s = fraction
      v = 1.0 - s/2.0
      if (A .le. 0.0) then
          h = 300.
          v = 1.0
      endif
      call hsv2rgb( h, s, v, r, g, b )
      return
      end


CCC   Color format conversion from Hue/Saturation/Value to Red/Green/Blue
CC    minimal (i.e. NO) error checking
C
      subroutine hsv2rgb( h, s, v, r, g, b )
      real                h, s, v, r, g, b
c
      real    f, p, q, t
      integer i
c
      i = h /60.
      f = h/60. - float(i)
      p = v * (1. - s)
      q = v * (1. - s*f)
      t = v * (1. - s*(1. - f))
      if (i.eq.5) then
          r = v
          g = p
          b = q
      else if (i.eq.4) then
          r = t
          g = p
          b = v
      else if (i.eq.3) then
          r = p
          g = q
          b = v
      else if (i.eq.2) then
          r = p
          g = v
          b = t
      else if (i.eq.1) then
          r = q
          g = v
          b = p
      else
          r = v
          g = t
          b = p
      endif
      return 
      end

CCC   This subroutine is not strictly necessary.
CC    It smooths the histogram/curve of Anisotropy by 
C     distance from center of mass

      subroutine hisclean( value, number, shells, nanis )
c
      integer shells
      real    value(shells)
      integer number(shells)
      integer nanis
c
      integer minshell, maxshell, min10, max10
      integer nsum
      real    vsum
c
      if (nanis .lt. 1200) then
          shells = shells / 2
          do i = 1, shells
            number(i) = number(2*i-1) + number(2*i)
            value(i)  = value(2*i-1) + value(2*i)
          enddo
      endif
c
      do i = shells, 1, -1
          if (number(i).gt.0) minshell = i
      enddo
      do i = 1, shells
          if (number(i).gt.0) maxshell = i
      enddo
      nsum = 0
      min10 = minshell
      do i = minshell, shells
          if (nsum + number(i) .gt. 10) then
            goto 101
          else
            nsum  = nsum + number(i)
            min10 = i
          endif
      enddo
101   continue
      nsum = 0
      max10 = maxshell
      do i = maxshell, 1, -1
          if (nsum + number(i) .gt. 10) then
            goto 102
          else
            nsum  = nsum + number(i)
            max10 = i
          endif
      enddo
102   continue
c
      nsum = 0
      vsum = 0.0
      do i = minshell, min10
          nsum = nsum + number(i)
          vsum = vsum + value(i)
          number(i) = 0
          value(i)  = 0.0
      enddo
      i = (minshell + min10) / 2
      number(i) = nsum
      value(i)  = vsum
c
      nsum = 0
      vsum = 0.0
      do i = maxshell, max10, -1
          nsum = nsum + number(i)
          vsum = vsum + value(i)
          number(i) = 0
          value(i)  = 0.0
      enddo
      i = (maxshell + max10) / 2
      number(i) = nsum
      value(i)  = vsum
c
      return
      end


Generated by  Doxygen 1.6.0   Back to index