PROGRAM pdbtodiam

c     Make Diamond files from PDB files. Specific for I2 setting of PpPCD

      CHARACTER*80      line
      CHARACTER*132     filnam     
      CHARACTER*4       label
      CHARACTER*3       resnam
      REAL              x(3), y(3)
      REAL              cell(6)/1., 1., 1., 90., 97.64, 90./

c     Get file name

      WRITE (6, 1000) ' Enter PDB mnemonic:'
1000  FORMAT (6a)
      READ (5, 1000) filnam
      len = length(filnam)
      WRITE (6, 1000) filnam(1:len)

c     Open files

      OPEN (unit=1, form='FORMATTED', status='OLD', readonly,
     1     file=filnam(1:len)//'.pdb')
      OPEN (unit=2, form='FORMATTED', status='UNKNOWN',
     1     file=filnam(1:len)//'.diam')
c     OPEN (unit=3, form='FORMATTED', status='UNKNOWN',
c    1     file=filnam(1:len)//'ca.diam')

c     Loop through remark cards looking for unit cell parameters

      j = 0
      DO WHILE (j .EQ. 0)
        READ (1, 1000, end=45) line
        IF (line(1:6) .EQ. 'REMARK') THEN
          j = INDEX(line, 'beta=')
        endIF
      endDO
      READ (line(j+5:j+10), 1001) cell(5)
45    REWIND 1
      tanbet = TAND(cell(5)-90.)
      cosbet = COSD(cell(5)-90.)

c     Write Diamond file headers

      WRITE (2, 1001) cell
1001  FORMAT (6f10.3)
c     WRITE (3, 1001) cell
      WRITE (2, 1000) ' Conversion of PDB file -- ', filnam
c     WRITE (3, 1000) ' Conversion of PDB file -- ', filnam
      WRITE (2, 1000) '      X         Y         Z          B',
     1     '    TYPE IRES ENTRY WEIGHT SEQ RES   ATOM'
c     WRITE (3, 1000) '      X         Y         Z          B',
c    1     '    TYPE IRES ENTRY WEIGHT SEQ RES   ATOM'

c     Loop through files

      nresold = -1
      mres = 0
      DO WHILE (.TRUE.)
5       READ (1, 1000, end=9) line
        IF (line(1:4) .NE. 'ATOM') goto 5

        READ (line, 1002, err=39) natom, label, resnam, nres, x, q, b
1002    FORMAT (6x, i5, 2x, a4, a3, 2x, i4, 4x, 3f8.3, 2f6.2)

        IF (label(1:1) .EQ. 'C') THEN
          itype = 6
        ELSE IF (label(1:1) .EQ. 'N') THEN
          itype = 7
        ELSE IF (label(1:1) .EQ. 'O') THEN
          itype = 8
        ELSE IF (label(1:1) .EQ. 'S') THEN
         itype = 16
        ELSE IF (label(1:2) .EQ. 'CL') THEN
          itype = 17
        ELSE IF (label(1:2) .EQ. 'FE') THEN
          itype = 23
        ELSE IF (label(1:2) .EQ. 'BR') THEN
          itype = 35
        ELSE IF (label(1:1) .EQ. 'I') THEN
          itype = 53
        ELSE IF (q .EQ. 0.0) THEN
          itype = 0 
        endIF

        IF (nres .NE. nresold) THEN
          mres = mres + 1
          nresold = nres
        endIF

c     Transform coordinates

        y(1) = x(1) + tanbet*x(3)
        y(2) = x(2)
        y(3) = x(3)/cosbet

        WRITE (2, 1003) y, b, itype, mres, natom, q, resnam, nres, 
     1     line(73:73), label
1003    FORMAT (4f10.3, 3i5, f8.2, 2x, a3, i4, a1, 2x, a4)
c       IF (label .EQ. 'CA  ') THEN
c         WRITE (3, 1003) x, b, itype, nres, natom, q, resnam, 
c    1      nres, line(73:73), label     
c       endIF
      endDO

c     Close files

39    WRITE (6, 1000) ' Error in reading following line:'
      WRITE (6, 1000) line
9     CLOSE (unit=1)
      CLOSE (unit=2)
c     CLOSE (unit=3)


      end
 
      INTEGER FUNCTION length(line)

      CHARACTER*132   line

      DO i = 132,1,-1
        IF (line(i:i) .NE. ' ') goto 9
      endDO

9     length = i

      return
      end