C **********************************************************************
C *  PDB
C=**********************************************************************
C=* From: SPIDER - MODULAR IMAGE PROCESSING SYSTEM                     *
C=* Copyright (C)2004, P. A. Penczek                                   *
C=*                                                                    *
C=* University of Texas - Houston Medical School                       *
C=* Email:  pawel.a.penczek@uth.tmc.edu                                *
C=*                                                                    *
C=* This program is free software; you can redistribute it and/or      *
C=* modify it under the terms of the GNU General Public License as     *
C=* published by the Free Software Foundation; either version 2 of the *
C=* License, or (at your option) any later version.                    *
C=*                                                                    *
C=* This program is distributed in the hope that it will be useful,    *
C=* but WITHOUT ANY WARRANTY; without even the implied warranty of     *
C=* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
C=* General Public License for more details.                           *
C=*                                                                    *
C=* You should have received a copy of the GNU General Public License  *
C=* along with this program; if not, write to the                      *
C=* Free Software Foundation, Inc.,                                    *
C=* 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.      *
C=*                                                                    *
C=**********************************************************************

	SUBROUTINE PDB

        INCLUDE 'CMBLOCK.INC'

        IF (FCHAR(4:6) .EQ. 'CG3') THEN
C          ------------------------- 'PDB GRAVITY'
           CALL PDBCG3

        ELSEIF (FCHAR(4:5) .EQ. 'IF') THEN
C          ------------------------- 'GET PDB FILE PARAMETERS'
	   CALL PDBIF

	ELSEIF (FCHAR(4:5) .EQ. 'SH') THEN
C          ------------------------- 'PDB SHIFT
	   CALL PDBSH

	ELSEIF (FCHAR(4:7) .EQ. 'RT3A') THEN
C          ------------------------- 'PDB ROTATE AROUND AN ARBITRARY POINT
	   CALL PDBRT3A

	ELSEIF (FCHAR(4:7) .EQ. 'RT3L') THEN
C          ------------------------- 'PDB ROTATE AROUND A LINE DEFINED by Two points
	   CALL PDBRT3L

	ELSEIF (FCHAR(4:6) .EQ. 'RT3') THEN
C          ------------------------- 'PDB ROTATE AROUND AN ARBITRARY POINT
	   CALL PDBRT3

	ELSEIF (FCHAR(4:6) .EQ. 'STP') THEN
C          ------------------------- 'COPY SPIDER COORDINATES TO PDB FORMAT'
	   CALL SPTOPDB

	ELSEIF (FCHAR(4:6) .EQ. 'PTS') THEN
C          ------------------------- 'COPY PDB FILE TO SPIDER FORMAT'
	   CALL PDBTOSP
	ENDIF
	END

C **********************************************************************

	SUBROUTINE PDBSH

        INCLUDE 'CMLIMIT.INC'
	INCLUDE 'CMBLOCK.INC'
        CHARACTER (LEN=MAXNAM) :: PDBFILE,TRFILE
	CHARACTER *80  RECLIN,HEAD*10, NULL*1,ATOM*4,RESIDUE*3
	LOGICAL        EX,FLAGCELLDIM,FLAGATOM,ZNUM,COMBIND

	DATA  LUN2,LUN3/25,26/

	NATOM  = 1
	LENREC = 0
        CALL OPAUXFILE(.TRUE.,PDBFILE,NULL,LUN2,LENREC,'O',
     &                 'PDB INPUT',.TRUE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

        CALL OPAUXFILE(.TRUE.,PDBFILE,NULL,LUN3,LENREC,'N',
     &                 'PDB OUTPUT',.TRUE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN
	
	CALL RDPRM3S(SXT,SYT,SZT,NOT_USED,
     &               'SHIFT IN ANGSTROMS SXT, SYT, SZT: ',IRTFLG)

C       CHANGE COORDINATE SYSTEM 
 	SX =  SYT
 	SY =  SXT
 	SZ = -SZT

70	READ(LUN2,30) RECLIN
30	FORMAT(A80)

	IF (RECLIN(1:4) .NE. 'ATOM' .AND.
     &      RECLIN(1:4) .NE. 'END'  .AND.
     & 	    RECLIN(1:4) .NE. 'TER' ) THEN
	    WRITE(LUN3,30) RECLIN
	    GOTO 70 

	ELSEIF (RECLIN(1:4) .EQ. 'ATOM') THEN
	   READ(RECLIN,75) HEAD,N,ATOM,NULL,RESIDUE,NULL,NR2,
     &                     NULL,XO,YO,ZO,
     &                     OCCUPANCY,TEMPERATURE,N
75	   FORMAT(A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,1X,I3)	

	   X = XO + SX
	   Y = YO + SY
	   Z = ZO + SZ

C REMOVED ALTERATION DEC 07 al
C          ALTERED NOV 07 al AFTER RECEIVING BUG REPORT
C           TEMP = X
C           X    = Y
C           Y    = TEMP
C           Z    = -Z 

	   WRITE(RECLIN(7:11),80) NATOM
80	   FORMAT(I5)

	   NATOM = NATOM + 1
	   WRITE(RECLIN(31:54),85) X,Y,Z
85	   FORMAT(3F8.3)

	   WRITE(LUN3,30) RECLIN
	   GOTO 70

	ELSEIF (RECLIN(1:3) .EQ. 'TER') THEN
	   WRITE(RECLIN(7:11),80)NATOM
	   WRITE(LUN3,30)RECLIN
	   NATOM = NATOM+1
	   GOTO 70

        ELSEIF (RECLIN(1:3) .EQ. 'END') THEN
	   WRITE(LUN3,30)RECLIN
	   GOTO 999
	ELSE
	   GOTO 70
	ENDIF

999	CLOSE(LUN2)
	CLOSE(LUN3)
	END


C **********************************************************************
C*
C*Rotate PDB file around an arbitrary center
C*
C********************************************

	SUBROUTINE PDBRT3A

	INCLUDE 'CMBLOCK.INC'
	INCLUDE 'CMLIMIT.INC'
        CHARACTER (LEN=MAXNAM) :: PDBFILE,TRFILE  
      
        CHARACTER *80 RECLIN,HEAD*10, NULL*1,ATOM*4,RESIDUE*3
	LOGICAL       EX,FLAGCELLDIM,FLAGATOM,ZNUM,COMBIND
	DOUBLE PRECISION RM(3,3)

	DATA  LUN2,LUN3/25,26/

	NATOM=1
	LENREC=0
	CALL OPAUXFILE(.TRUE.,PDBFILE,NULL,LUN2,LENREC,'O',
     &                 'PDB INPUT',.TRUE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN
 
	CALL OPAUXFILE(.TRUE.,PDBFILE,NULL,LUN3,LENREC,'N',
     &                 'PDB OUTPUT',.TRUE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

	CALL RDPRM3S(PHI,THETA,PSI,NOT_USED,
     &     'ROTATION ANGLE IN DEGREES PHI, THETA, PSI: ',IRTFLG)

	CALL  BLDR(RM,PSI,THETA,PHI)

	CALL RDPRM3S(XC,YC,ZC,NOT_USED,
     &     'CENTER OF ROTATION IN ANGSTROMS X,Y,Z: ',IRTFLG)

70	READ(LUN2,30) RECLIN
30	FORMAT(A80)

	IF (RECLIN(1:4) .NE. 'ATOM' .AND.
     &      RECLIN(1:4).NE.'END'.AND.
     & 	    RECLIN(1:4).NE.'TER' ) THEN
	   WRITE(LUN3,30) RECLIN
	   GOTO 70 

	ELSEIF (RECLIN(1:4) .EQ. 'ATOM') THEN
	   READ(RECLIN,75) HEAD,N,ATOM,NULL,RESIDUE,NULL,NR2,
     &                     NULL,XO,YO,ZO,
     &                     OCCUPANCY,TEMPERATURE,N
75	   FORMAT(A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,1X,I3)	
C	   WRITE(*,*) XO,YO,ZO

C          CHANGE COORDINATE SYSTEM
	   W  = XO
	   XO = YO
	   YO = W
	   ZO = -ZO
	   XO = XO-XC
	   YO = YO-YC
	   ZO = ZO-ZC

C          AFTER ROTATION CHANGE IT BACK, THAT'S WHY THE ORDER OF 
C          X,Y IS CHANGED AND Z HAS INVERTED SIGN

	   Y = RM(1,1)*XO+RM(1,2)*YO+RM(1,3)*ZO+XC
	   X = RM(2,1)*XO+RM(2,2)*YO+RM(2,3)*ZO+YC
	   Z = -(RM(3,1)*XO+RM(3,2)*YO+RM(3,3)*ZO+ZC)
	   WRITE(RECLIN(7:11),80) NATOM
80	   FORMAT(I5)
	   NATOM = NATOM+1
	   WRITE(RECLIN(31:54),85)  X,Y,Z
85	   FORMAT(3F8.3)
	   WRITE(LUN3,30)RECLIN
	   GOTO 70

	ELSEIF (RECLIN(1:3) .EQ. 'TER') THEN
	   WRITE(RECLIN(7:11),80)NATOM
	   WRITE(LUN3,30)RECLIN
	   NATOM = NATOM + 1
	   GOTO 70

        ELSEIF (RECLIN(1:3) .EQ. 'END') THEN
	   WRITE(LUN3,30)RECLIN
	   GOTO 999

	ELSE
	   GOTO 70
	ENDIF

999	CLOSE(LUN2)
	CLOSE(LUN3)
	END


C **********************************************************************
C*
C*Rotate PDB file around a line defined by two points in 3D space
C*
C********************************************

	SUBROUTINE PDBRT3L

	INCLUDE 'CMBLOCK.INC'
	INCLUDE 'CMLIMIT.INC'
        CHARACTER (LEN=MAXNAM) :: PDBFILE,TRFILE
        CHARACTER *80 RECLIN,HEAD*10, NULL*1,ATOM*4,RESIDUE*3
	LOGICAL       EX,FLAGCELLDIM,FLAGATOM,ZNUM,COMBIND
	DOUBLE PRECISION RM(3,3),R1(3,3),R2(3,3),R3(3,3)

	DATA  LUN2,LUN3/25,26/

	NATOM  = 1
	LENREC = 0
	CALL OPAUXFILE(.TRUE.,PDBFILE,NULL,LUN2,LENREC,'O',
     &                 'PDB INPUT',.TRUE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

30	FORMAT(A80)
 
	CALL OPAUXFILE(.TRUE.,PDBFILE,NULL,LUN3,LENREC,'N',
     &                 'PDB OUTPUT',.TRUE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

	CALL RDPRM1S(ALPHA,NOT_USED,'ROTATION ANGLE ALPHA',IRTFLG)

	CALL RDPRM3S(X1,Y1,Z1,NOT_USED,
     &        'FIRST POINT DEFINING ROTATION AXIS IN ANGSTROMS X,Y,Z',
     &        IRTFLG)

	CALL RDPRM3S(X2,Y2,Z2,NOT_USED,
     &       'SECOND POINT DEFINING ROTATION AXIS IN ANGSTROMS X,Y,Z',
     &       IRTFLG)

	XX    = X2-X1
	YY    = Y2-Y1
	ZZ    = Z2-Z1
	PSI   = -ATAN2D(YY,XX)
	THETA = ATAN2D(ZZ,SQRT(XX*XX+YY*YY))

	CALL  BLDR(R1,PSI,THETA,90.0)
	CALL  BLDR(R2,0.0,ALPHA,0.0)

C       -1
C       R   U(ALPHA) R

         DO  11  I=1,3
            DO  11  J=1,3
               R3(J,I)=0.0
               DO  11  K=1,3
11       R3(J,I)=R3(J,I)+R2(K,I)*R1(J,K)

         DO  12  I=1,3
            DO  12  J=1,3
               RM(J,I)=0.0
               DO  12  K=1,3
12       RM(J,I)=RM(J,I)+R1(I,K)*R3(J,K)

70	READ(LUN2,30) RECLIN
	IF (RECLIN(1:4) .NE. 'ATOM' .AND.
     &      RECLIN(1:4) .NE. 'END'  .AND.
     &      RECLIN(1:4) .NE. 'TER' ) THEN
	    WRITE(LUN3,30) RECLIN
	    GOTO 70 

	ELSEIF (RECLIN(1:4) .EQ. 'ATOM') THEN
	   READ(RECLIN,75) HEAD,N,ATOM,NULL,RESIDUE,NULL,NR2,
     &                     NULL,XO,YO,ZO,
     &                     OCCUPANCY,TEMPERATURE,N
75	   FORMAT(A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,1X,I3)	
C	   WRITE(*,*) XO,YO,ZO

C          CHANGE COORDINATE SYSTEM
	   W  = XO
	   XO = YO
	   YO = W
	   ZO = -ZO
	   XO = XO-X1
	   YO = YO-Y1
	   ZO = ZO-Z1

C          AFTER ROTATION CHANGE IT BACK, THAT'S WHY THE ORDER OF X,Y IS CHANGED
C          AND Z HAS INVERTED SIGN

	   Y=RM(1,1)*XO+RM(1,2)*YO+RM(1,3)*ZO+X1
	   X=RM(2,1)*XO+RM(2,2)*YO+RM(2,3)*ZO+Y1
	   Z=-(RM(3,1)*XO+RM(3,2)*YO+RM(3,3)*ZO+Z1)
	   WRITE(RECLIN(7:11),80) NATOM
80	   FORMAT(I5)
	   NATOM = NATOM+1
	   WRITE(RECLIN(31:54),85)  X,Y,Z
85	   FORMAT(3F8.3)
	   WRITE(LUN3,30)RECLIN
	   GOTO 70

	ELSEIF (RECLIN(1:3) .EQ. 'TER') THEN
	   WRITE(RECLIN(7:11),80)NATOM
	   WRITE(LUN3,30)RECLIN
	   NATOM = NATOM + 1
	   GOTO 70

        ELSEIF (RECLIN(1:3) .EQ. 'END') THEN
	   WRITE(LUN3,30)RECLIN
	   GOTO 999

	ELSE
	   GOTO 70
	ENDIF

999	CLOSE(LUN2)
	CLOSE(LUN3)
	END

C **********************************************************************
C*
C* Copy SPIDER document file to a PDB file
C*
C********************************************

	SUBROUTINE SPTOPDB

	INCLUDE 'CMBLOCK.INC'
	INCLUDE 'CMLIMIT.INC'
        CHARACTER (LEN=MAXNAM) :: PDBFILE,TRFILE
 	CHARACTER *95 RECLIN,HEAD*10, NULL*1,ATOM*4
	CHARACTER *95 ANAM*4,RESIDUE*3,PHEAD*5,SENQ*1
	LOGICAL EX,FLAGCELLDIM,FLAGATOM,ZNUM,COMBIND

	DATA  LUN2,LUN3/25,26/

	OPEN(99,FILE='ttt')

	PHEAD='ATOM '
	LENREC = 0
30	FORMAT(A80)
        CALL OPAUXFILE(.TRUE.,PDBFILE,NULL,LUN2,LENREC,'O',
     &                 'SPIDER INPUT',.TRUE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

        CALL OPAUXFILE(.TRUE.,PDBFILE,NULL,LUN3,LENREC,'N',
     &                 'PDB OUTPUT',.TRUE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

	NCOUNT=0
	DO
55	READ(LUN2,30,END=701) TRFILE
	IF (TRFILE(1:2) .EQ.' ;') GOTO 55
	
	READ(TRFILE,33)NATOM,MMM,XO,YO,ZO,CATOM,TEMPERATURE,XSEN,XNUL

cc==decoding

	INUL=INT(XNUL)
	IRES=INT(CATOM/10000.)
	NR2=INT(XSEN/1000.)
	ISEN=INT(REAL(MOD(INT(XSEN),1000))/10.)
	OCCUPANCY=REAL(MOD(INT(XSEN),2))	
	    IF (IRES.EQ.1) THEN
		RESIDUE='LEU'
	    ElSEIF (IRES.EQ.2) THEN
		RESIDUE='THR'
	    ELSEIF (IRES.EQ.3) THEN
		RESIDUE='GLY'
	    ELSEIF (IRES.EQ.4) THEN
		RESIDUE='SER'
	    ELSEIF (IRES.EQ.5) THEN
		RESIDUE='VAL'
	    ELSEIF (IRES.EQ.6) THEN
		RESIDUE='PHE'
	    ELSEIF (IRES.EQ.7) THEN
		RESIDUE='ALA'
	    ELSEIF (IRES.EQ.8) THEN
		RESIDUE='LYS'
       	    ELSEIF (IRES.EQ.9) THEN
		RESIDUE='ARG'
	    ELSEIF (IRES.EQ.10) THEN
		RESIDUE='PRO'
            ELSEIF (IRES.EQ.11) THEN
		RESIDUE='TAU'
	    ELSEIF (IRES.EQ.12) THEN
		RESIDUE='HIS'
    	    ELSEIF (IRES.EQ.13) THEN
		RESIDUE='GLU'
	    ELSEIF (IRES.EQ.14) THEN
		RESIDUE='GLY'
	    ELSEIF (IRES.EQ.15) THEN
		RESIDUE='TYR'
	    ELSEIF (IRES.EQ.16) THEN
		RESIDUE='CYS'
  	    ELSEIF (IRES.EQ.17) THEN
		RESIDUE='TRY'
	    ELSEIF (IRES.EQ.18) THEN
		RESIDUE='ISO'
	    ELSEIF (IRES.EQ.19) THEN
		RESIDUE='MET'
	    ELSEIF (IRES.EQ.20) THEN
		RESIDUE='ILE'
	    ENDIF
cc
 	    IF (isen.EQ.1) THEN
		SENQ='A'
	    ELSEIF (isen.EQ.2) THEN
		SENQ='B'
	    ELSEIF (isen.EQ.3) THEN
		SENQ='C'
	    ELSEIF (isen.EQ.4) THEN
		SENQ='D'
	    ELSEIF (isen.EQ.5) THEN
		SENQ='E'
	    ELSEIF (isen.EQ.6) THEN
		SENQ='F'
	    ELSEIF (isen.EQ.7) THEN
		SENQ='G'
	    ELSEIF (isen.EQ.8) THEN
		SENQ='H'
     	    ELSEIF (isen.EQ.9) THEN
		SENQ='I'
	    ELSEIF (isen.EQ.10) THEN
		SENQ='J'
	    ELSEIF (isen.EQ.11) THEN
		SENQ='K'
	    ELSEIF (isen.EQ.12) THEN
		SENQ='L'
	    ELSEIF (isen.EQ.13) THEN
		SENQ='M'
	   ELSEIF (isen.EQ.14) THEN
		SENQ='N'
	    ELSEIF (isen.EQ.15) THEN
		SENQ='O'
	    ELSEIF (isen.EQ.16) THEN
		SENQ='P'
	    ELSEIF (isen.EQ.17) THEN
		SENQ='Q'
	   ENDIF
C
	   IF (INUL.EQ.0) THEN
		NULL=' '
	    ELSEIF (INUL.EQ.1) THEN
		NULL='A'
	    ELSEIF (INUL.EQ.2) THEN
		NULL='B'
	    ELSEIF (INUL.EQ.3) THEN
		NULL='C'
	    ELSEIF (INUL.EQ.4) THEN
		NULL='D'
	    ELSEIF (INUL.EQ.5) THEN
		NULL='E'
	    ELSEIF (INUL.EQ.6) THEN
		NULL='F'
	    ELSEIF (INUL.EQ.7) THEN
		NULL='G'
	    ELSEIF (INUL.EQ.8) THEN
		NULL='H'
     	    ELSEIF (INUL.EQ.9) THEN
		NULL='I'
	    ELSEIF (INUL.EQ.10) THEN
		NULL='J'
	    ELSEIF (INUL.EQ.11) THEN
		NULL='K'
	    ELSEIF (INUL.EQ.12) THEN
		NULL='L'
	    ELSEIF (INUL.EQ.13) THEN
		NULL='M'
	    ELSEIF (INUL.EQ.14) THEN
		NULL='N'
	    ELSEIF (INUL.EQ.15) THEN
		NULL='O'
	    ELSEIF (INUL.EQ.16) THEN
		NULL='P'
	    ELSEIF (INUL.EQ.17) THEN
		NULL='Q'
	    ELSEIF (INUL.EQ.18) THEN
		NULL='R'
	    ELSEIF (INUL.EQ.19) THEN
		NULL='S'
	    ELSEIF (INUL.EQ.20) THEN
		NULL='T'
	    ELSEIF (INUL.EQ.21) THEN
		NULL='U'
	    ELSEIF (INUL.EQ.22) THEN
		NULL='V'
	   ELSE
	        NULL=' '
	   ENDIF
	   NCATOM = INT(REAL(MOD(INT(CATOM),10000))/100.)	
	   NPOS   = MOD(INT(CATOM),100)
	
	   IF (NCATOM.EQ.1) THEN
		ANAM=' H  '
	   ELSEIF (NCATOM.EQ.6) THEN
		
		IF (NPOS.EQ.0) THEN
		    ANAM=' C  '
		ELSEIF (NPOS.EQ.1) THEN
		    ANAM=' CA '
		ELSEIF (NPOS.EQ.2) THEN
		    ANAM=' CB '
		ELSEIF (NPOS.EQ.3) THEN
		    ANAM=' CE1'
		ELSEIF (NPOS.EQ.4) THEN
		    ANAM=' CE2'
		ELSEIF (NPOS.EQ.5) THEN
		    ANAM=' CD '
		ELSEIF (NPOS.EQ.6) THEN
		    ANAM=' CD1'
		ELSEIF (NPOS.EQ.7) THEN
		    ANAM=' CD2'
		ELSEIF (NPOS.EQ.8) THEN
		    ANAM=' CG '
		ELSEIF (NPOS.EQ.9) THEN
		    ANAM=' CG1'
		ELSEIF (NPOS.EQ.10) THEN
		    ANAM=' CG2'
		ENDIF
	   ELSEIF (NCATOM.EQ.7) THEN
		IF (NPOS.EQ.0) THEN
		    ANAM=' N  '
		ELSEIF (NPOS.EQ.1) THEN
		    ANAM=' N1'
		ELSEIF (NPOS.EQ.2) THEN
		    ANAM=' NH2'
		ELSEIF (NPOS.EQ.3) THEN
		    ANAM=' NE2'
		ELSEIF (NPOS.EQ.4) THEN
		    ANAM=' NZ '
		ENDIF
	   ELSEIF (NCATOM.EQ.8) THEN
		IF (NPOS.EQ.0) THEN
		    ANAM=' O  '
		ELSEIF (NPOS.EQ.1) THEN
		    ANAM=' OH'
		ELSEIF (NPOS.EQ.2) THEN
		    ANAM=' OE1'
		ELSEIF (NPOS.EQ.3) THEN
		    ANAM=' OE2'
		ELSEIF (NPOS.EQ.4) THEN
		    ANAM=' OG1'
		ELSEIF (NPOS.EQ.5) THEN
		    ANAM=' OD1'
		ELSEIF (NPOS.EQ.6) THEN
		    ANAM=' OXT'
		ENDIF
	   ELSEIF (NCATOM.EQ.16) THEN
		IF (NPOS.EQ.0) THEN
		    ANAM=' S  '
		ELSEIF (NPOS.EQ.1) THEN
		    ANAM=' SG '
		ELSEIF (NPOS.EQ.2) THEN
		    ANAM=' SD '
		ENDIF
	   ENDIF
	   NCOUNT=NCOUNT+1

	   WRITE(98,*)XNUL,INUL
	   WRITE(LUN3,76)PHEAD,NCOUNT,ANAM,RESIDUE,SENQ,NR2,
     &                   NULL,XO,YO,ZO,
     &                   OCCUPANCY, TEMPERATURE,ANAM(2:2)
76	   FORMAT(A5,1X,I5,1X,A4,1X,A3,1X,A1,I4,1A,3X,3F8.3,2F6.2,
     &        '           ',1A)
	   IF (NCATOM.EQ.8.AND.NPOS.EQ.6) THEN
	      NCOUNT=NCOUNT+1
	      WRITE(LUN3,77)NCOUNT,RESIDUE,SENQ,NR2
	   ENDIF
77	   FORMAT('TER   ',I5,6X,A3,1X,A1,I4)	
	ENDDO	

33	FORMAT(I5,1X,I1,7G12.6)
701     WRITE(LUN3,62)
62	FORMAT('END')

	CLOSE(LUN2)
	CLOSE(LUN3)
	CLOSE(99)
	END

C **********************************************************************
C*
C* Copy PDB file to a SPIDER document file
C*
C********************************************

	SUBROUTINE PDBTOSP

	INCLUDE 'CMBLOCK.INC'
	INCLUDE 'CMLIMIT.INC'

	PARAMETER  (NLIST=8)
      	REAL	DLIST(NLIST)
        CHARACTER (LEN=MAXNAM) :: PDBFILE,TRFILE   
        CHARACTER *95 RECLIN,HEAD*10, NULL*1,ATOM*4,RESIDUE*3
	CHARACTER *95 NULL1*1,NULL2*1
	LOGICAL EX,FLAGCELLDIM,FLAGATOM,ZNUM,COMBIND

	DATA  LUN2,LUN3/25,26/

30	FORMAT(A95)
	LENREC = 0
	CALL OPAUXFILE(.TRUE.,PDBFILE,NULL,LUN2,LENREC,'O',
     &                 'PDB INPUT',.TRUE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

31	FORMAT(' ;  ',A60)
	NATOM=0

70	READ(LUN2,30) RECLIN
	IF (RECLIN(1:4) .EQ. 'ATOM') THEN
	   READ(RECLIN,75)HEAD,N,ATOM,NULL,RESIDUE,NULL1,NR2,
     &          NULL2,XO,YO,ZO,
     &         OCCUPANCY,TEMPERATURE,N
75	   FORMAT(A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,1X,I3)	

C I AM NOT SURE ABOUT THE NEXT LINE....
C	     IF (OCCUPANCY.GE.1) THEN

		IF (ATOM(2:2) .EQ. 'H') THEN
		   VATOM=100.
CCarbon
		ELSEIF (ATOM(2:2) .EQ. 'C') THEN
		   VATOM=600.
			IF (ATOM(3:3).EQ.'A') THEN
				 VATOM=601
			ELSEIF (ATOM(3:3).EQ.'B') THEN
				 VATOM=602
			ELSEIF (ATOM(3:4).EQ.'E1')THEN
			  	VATOM=603
			ELSEIF (ATOM(3:4).EQ.'E2') THEN
				VATOM=604
			ELSEIF (ATOM(3:3).EQ.'D') THEN	
				VATOM=605
			ELSEIF (ATOM(3:4).EQ.'D1') THEN
				VATOM=606
			ELSEIF (ATOM(3:4).EQ.'D2') THEN
				VATOM=607				
			ELSEIF (ATOM(3:3).EQ.'G') THEN
				VATOM=608
			ELSEIF (ATOM(3:4).EQ.'G1') THEN
				VATOM=609	
			ELSEIF (ATOM(3:4).EQ.'G2') THEN
				VATOM=610
			ELSEIF (ATOM(3:3).EQ.'Z') THEN
				VATOM=611
			ENDIF
CNitrogen
		ELSEIF (ATOM(2:2) .EQ. 'N') THEN
		   VATOM=700.
			IF (ATOM(3:4).EQ.'H1') THEN
			  	VATOM=701
			    ELSEIF (ATOM(3:4).EQ.'H2') THEN
				VATOM=702.
			    ELSEIF (ATOM(3:3).EQ.'E')THEN
				VATOM=703.
			    ELSEIF (ATOM(3:4).EQ.'E2') THEN
				VATOM=704.
			    ELSEIF (ATOM(3:3).EQ.'Z') THEN
				VATOM=705.
			ENDIF
CS	
	        ELSEIF (ATOM(2:2) .EQ. 'S') THEN
		   VATOM=1600.
			IF (ATOM(3:3).EQ.'G') THEN
				VATOM=1601.
			ELSEIF (ATOM(3:3).EQ.'D') THEN
				VATOM=1602.
			ENDIF
CP	
		ELSEIF (ATOM(2:2) .EQ. 'P') THEN
		      VATOM=1500.
COxygen
		ELSEIF (ATOM(2:2) .EQ. 'O') THEN
		   VATOM=800.
		   IF (ATOM(3:3) .EQ. 'H') THEN
			 VATOM=801.
		   ELSEIF (ATOM(3:4) .EQ. 'E1') THEN
			 VATOM=802.
		   ELSEIF (ATOM(3:4) .EQ. 'E2') THEN
			VATOM=803.
		   ELSEIF (ATOM(3:4) .EQ. 'G1') THEN
			VATOM=804.
		   ELSEIF (ATOM(3:4) .EQ. 'D1') THEN
			VATOM=805.
		   ELSEIF (ATOM(3:4) .EQ. 'XT') THEN
			VATOM=806.
		   ENDIF
		ENDIF
CResidue
	    IF (RESIDUE(1:3).EQ.'LEU') THEN
		IRES=1
	    ElSEIF (RESIDUE(1:3).EQ.'THR') THEN
		IRES=2
	    ELSEIF (RESIDUE(1:3).EQ.'GLY') THEN
		IRES=3
	    ELSEIF (RESIDUE(1:3).EQ.'SER') THEN
		IRES=4
	    ELSEIF (RESIDUE(1:3).EQ.'VAL') THEN
		IRES=5
	    ELSEIF (RESIDUE(1:3).EQ.'PHE') THEN
		IRES=6
	    ELSEIF (RESIDUE(1:3).EQ.'ALA') THEN
		IRES=7
	    ELSEIF (RESIDUE(1:3).EQ.'LYS') THEN
		IRES=8
       	    ELSEIF (RESIDUE(1:3).EQ.'ARG') THEN
		IRES=9
	    ELSEIF (RESIDUE(1:3).EQ.'PRO') THEN
		IRES=10
            ELSEIF (RESIDUE(1:3).EQ.'TAU') THEN
		IRES=11
	    ELSEIF (RESIDUE(1:3).EQ.'HIS') THEN
		IRES=12
    	    ELSEIF (RESIDUE(1:3).EQ.'GLU') THEN
		IRES=13
	    ELSEIF (RESIDUE(1:3).EQ.'GLY') THEN
		IRES=14
	    ELSEIF (RESIDUE(1:3).EQ.'TYR') THEN
		IRES=15
	    ELSEIF (RESIDUE(1:3).EQ.'CYS') THEN
		IRES=16
  	    ELSEIF (RESIDUE(1:3).EQ.'TRY') THEN
		IRES=17
	    ELSEIF (RESIDUE(1:3).EQ.'ISO') THEN
		IRES=18
	    ELSEIF (RESIDUE(1:3).EQ.'MET') THEN
		IRES=19
	    ELSEIF (RESIDUE(1:3).EQ.'ILE') THEN
		IRES=20
	    ENDIF
CCSequence
	    IF (NULL1(1:1).EQ.'A') THEN
		ISEN=1
	    ELSEIF (NULL1(1:1).EQ.'B') THEN
		ISEN=2
	    ELSEIF (NULL1(1:1).EQ.'C') THEN
		ISEN=3
	    ELSEIF (NULL1(1:1).EQ.'D') THEN
		ISEN=4
	    ELSEIF (NULL1(1:1).EQ.'E') THEN
		ISEN=5
	    ELSEIF (NULL1(1:1).EQ.'F') THEN
		ISEN=6
	    ELSEIF (NULL1(1:1).EQ.'G') THEN
		ISEN=7
	    ELSEIF (NULL1(1:1).EQ.'H') THEN
		ISEN=8
     	    ELSEIF (NULL1(1:1).EQ.'I') THEN
		ISEN=9
	    ELSEIF (NULL1(1:1).EQ.'J') THEN
		ISEN=10
	    ELSEIF (NULL1(1:1).EQ.'K') THEN
		ISEN=11
	    ELSEIF (NULL1(1:1).EQ.'L') THEN
		ISEN=12
	    ELSEIF (NULL1(1:1).EQ.'M') THEN
		ISEN=13
	   ELSEIF (NULL1(1:1).EQ.'N') THEN
		ISEN=14
	    ELSEIF (NULL1(1:1).EQ.'O') THEN
		ISEN=15
	    ELSEIF (NULL1(1:1).EQ.'P') THEN
		ISEN=16
	    ELSEIF (NULL1(1:1).EQ.'Q') THEN
		ISEN=17
	   ENDIF
	   IF (NULL2.EQ.' ') THEN
		INUL=0
	   ELSEIF (NULL2.EQ.'A') THEN
	        INUL=1	 	
	   ELSEIF (NULL2.EQ.'B') THEN
		INUL=2
	   ELSEIF (NULL2.EQ.'C') THEN
	        INUL=3	 	
	   ELSEIF (NULL2.EQ.'D') THEN
	        INUL=4	
	   ELSEIF (NULL2.EQ.'E') THEN
	        INUL=5	 	
	   ELSEIF (NULL2.EQ.'F') THEN
		INUL=6
	   ELSEIF (NULL2.EQ.'G') THEN
	        INUL=7	 	
	   ELSEIF (NULL2.EQ.'H') THEN
	        INUL=8	
	   ELSEIF (NULL2.EQ.'I') THEN
	        INUL=9	
	   ELSEIF (NULL2.EQ.'J') THEN
	        INUL=10	 	
	   ELSEIF (NULL2.EQ.'K') THEN
		INUL=11
	   ELSEIF (NULL2.EQ.'L') THEN
	        INUL=12	 	
	   ELSEIF (NULL2.EQ.'M') THEN
	        INUL=13
	   ELSEIF (NULL2.EQ.'N') THEN
	        INUL=14	 	
	   ELSEIF (NULL2.EQ.'O') THEN
		INUL=15
	   ELSEIF (NULL2.EQ.'P') THEN
	        INUL=16	 	
	   ELSEIF (NULL2.EQ.'Q') THEN
	        INUL=17	
	   ELSEIF (NULL2.EQ.'R') THEN
	        INUL=18	
	   ELSEIF (NULL2.EQ.'S') THEN
	        INUL=19	 	
	   ELSEIF (NULL2.EQ.'T') THEN
		INUL=20
	   ELSEIF (NULL2.EQ.'U') THEN
	        INUL=21	 	
	   ELSEIF (NULL2.EQ.'V') THEN
		INUL=22
	   ENDIF  
	    

           NATOM=NATOM+1

	   DLIST(1)=NATOM
	   DLIST(2)=XO
	   DLIST(3)=YO
           DLIST(4)=ZO
           DLIST(5)=IRES*10000.+VATOM
	   DLIST(6)=TEMPERATURE
           DLIST(7)=NR2*1000.+ISEN*10.+occupancy
	   DLIST(8)=INUL
	   CALL SAVD(LUN3,DLIST,NLIST,IRTFLG)
	   GOTO 70
        ELSEIF (RECLIN(1:3) .EQ. 'END') THEN
	   CLOSE(LUN2)
	   CALL SAVDC
	   CLOSE(LUN3)
	   RETURN
	ELSE
	   GOTO 70
	ENDIF
	END

C **********************************************************************
C*
C* Calculate center of gravity of a PDB file
C*
C********************************************

	SUBROUTINE PDBCG3

	INCLUDE 'CMBLOCK.INC'
        INCLUDE 'CMLIMIT.INC' 
	REAL, ALLOCATABLE, DIMENSION (:) :: R_TMP
        CHARACTER (LEN=MAXNAM) :: PDBFILE,TRFILE        
        CHARACTER *80 RECLIN,HEAD*10, NULL*1,ATOM*4,RESIDUE*3
	LOGICAL       EX,FLAGCELLDIM,FLAGATOM,ZNUM,COMBIND

	DATA  LUN2,LUN3/25,26/

	LENREC = 0
	CALL OPAUXFILE(.TRUE.,PDBFILE,NULL,LUN2,LENREC,'O',
     &                 'PDB INPUT',.TRUE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN
 
	NATOM=0
	UX=0.0
	UY=0.0
	UZ=0.0
	TMASS=0.0

70	READ(LUN2,30) RECLIN
30	FORMAT(A80)

	IF (RECLIN(1:6) .EQ. 'HEADER') THEN
           GOTO 70

	ELSEIF (RECLIN(1:4) .EQ. 'ATOM') THEN
	   READ(RECLIN,75) HEAD,N,ATOM,NULL,RESIDUE,NULL,NR2,
     &                     NULL,XO,YO,ZO,
     &                     OCCUPANCY,TEMPERATURE,N
75	   FORMAT(A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,1X,I3)	
C	   WRITE(*,*) XO,YO,ZO

C            I AM NOT SURE ABOUT THE NEXT LINE....
	     IF (OCCUPANCY .GE. 1) THEN
		IF (ATOM(2:2) .EQ. 'H') THEN
		   VATOM=1.0
		ELSEIF (ATOM(2:2) .EQ. 'C') THEN
		   VATOM=6.0
		ELSEIF (ATOM(2:2) .EQ. 'N') THEN
		   VATOM=7.0
		ELSEIF (ATOM(2:2) .EQ. 'O') THEN
		   VATOM=8.0
	        ELSEIF (ATOM(2:2) .EQ. 'S') THEN
		   VATOM=16.0
		ELSEIF (ATOM(2:2) .EQ. 'P') THEN
		   VATOM=15.0
		ELSEIF (ATOM(2:2) .EQ. 'X' .AND. ATOM(1:3) .EQ. 'OXT')
     &             THEN
C                  TRNA THREE PRIME TERMINAL
		   VATOM=8.0
		ELSE
		   WRITE(*, *) 'SPECIAL ATOM ENCOUNTERED IN:'
		   WRITE(*,*) RECLIN 
	           STOP
		ENDIF
C	        VATOM = VATOM*OCCUPANCY
	     ELSE 
	        VATOM=OCCUPANCY
	     ENDIF
	   TMASS=TMASS+VATOM
           UX=UX+XO*VATOM
           UY=UY+YO*VATOM
           UZ=UZ+ZO*VATOM
           NATOM=NATOM+1
	   GOTO 70

        ELSEIF (RECLIN(1:3) .EQ. 'END') THEN
	   UX=UX/TMASS
	   UY=UY/TMASS
	   UZ=UZ/TMASS

	   WRITE(*,780)  NATOM
780	   FORMAT(' Number of atoms encountered:',I6)

	   WRITE(*,781)  UX,UY,UZ
781	   FORMAT('  Center of gravity of PDB file',/,3(2x,g12.4))

C          CHANGE COORDINATE SYSTEM
	   SX = UY
	   SY = UX
	   SZ = -UZ

	   WRITE(*,782)  SX,SY,SZ
782	   FORMAT('  Center of gravity of PDB file in SPIDER coords'
     &		  ,/,3(2x,g12.4))

	   ALLOCATE (R_TMP(6))
	   NREG     = 7
	   R_TMP(1) = UX
	   R_TMP(2) = UY
	   R_TMP(3) = UZ
	   R_TMP(4) = SX
	   R_TMP(5) = SY
	   R_TMP(6) = SZ
	   R_TMP(7) = NATOM
	   CALL REG_SET_NSELA(NREG,R_TMP,.FALSE.,IRTFLG)
	   DEALLOCATE(R_TMP)
	   CLOSE(LUN2) 
	   RETURN	
	ELSE
	   GOTO 70
	ENDIF

	CLOSE(LUN2)
	END

C **********************************************************************
C*
C* Copy PDB file to a SPIDER document file
C*
C********************************************

	SUBROUTINE PDBIF

	INCLUDE 'CMBLOCK.INC'
	INCLUDE 'CMLIMIT.INC' 
	PARAMETER  (NLIST=4,NRG=6)
      	REAL	DLIST(NLIST)
	REAL	DLIST1(NLIST)
	REAL    R_TMP(NRG)   
        CHARACTER (LEN=MAXNAM) :: PDBFILE,TRFILE
        CHARACTER *80 RECLIN,HEAD*10
	CHARACTER *80 NULL*6,ATOM*4,RESIDUE*3
	LOGICAL       EX,FLAGCELLDIM,FLAGATOM,ZNUM,COMBIND

	DATA  LUN2,LUN3/25,26/

	LENREC=0
	CALL OPAUXFILE(.TRUE.,PDBFILE,NULL,LUN2,LENREC,'O',
     &                 'PDB INPUT',.TRUE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

30	FORMAT(A80)
 
31	FORMAT(' ;  ',A60)
33	FORMAT(A30)
35	FORMAT(A6,3(2X,F7.4),3F7.4)
36	FORMAT(A6,4X,3F10.4,5X,F10.4)
37	FORMAT(' ;',A21)
	NCOUNT=0
	NATOM=0
	NTERM=0
	NHET=0
	NPR=0

70	READ(LUN2,30) RECLIN
	IF ( RECLIN(1:4) .EQ. 'JRNL') THEN
	   WRITE(*,30)RECLIN         
	   GOTO 70

	ELSEIF (RECLIN(13:23) .EQ. 'FIT TO DATA') THEN
            WRITE(*,*) RECLIN(12:60)
            DO NPR=1,6
               READ(LUN2,30) RECLIN
               WRITE(*,*)RECLIN(12:60)
            ENDDO
            GOTO 70

	ELSEIF (RECLIN(13:21) .EQ. 'DATA USED') THEN
            WRITE(*,*) RECLIN(12:60)
            DO NPR=1,6
               READ(LUN2,30) RECLIN
               WRITE(*,*)RECLIN(12:60)
            ENDDO
            GOTO 70

	ELSEIF (RECLIN(14:20) .EQ. 'PROGRAM') THEN
            WRITE(*,*) 'USED PROGRAM',RECLIN(21:60)         
            GOTO 70

	ELSEIF (RECLIN(1:6) .EQ. 'CRYST1') THEN
	    WRITE(*,*) RECLIN	
	    READ(RECLIN(1:54),35)NULL,DLIST(2:4),DLIST1(2:4)
	    NCOUNT=NCOUNT+1
	    DLIST(1)=NCOUNT
	    R_TMP(1:3)=DLIST(2:4)	
	    CALL SAVD(LUN3,DLIST,NLIST,IRTFLG)
	    NCOUNT=NCOUNT+1
	    DLIST1(1)=NCOUNT
	    DLIST(:)=DLIST1(:)		
	    CALL SAVD(LUN3,DLIST,NLIST,IRTFLG)
	    GOTO 70

	ELSEIF (RECLIN(1:3) .EQ. 'ORI') THEN	
	    NCOUNT=NCOUNT+1
	    DLIST(1)=NCOUNT	
	    READ(RECLIN,36)NULL,DLIST(2:4),XX
	    CALL SAVD(LUN3,DLIST,NLIST,IRTFLG)
	    GOTO 70

	ELSEIF (RECLIN(1:3) .EQ. 'SCA') THEN
	    NCOUNT=NCOUNT+1
	    DLIST(1)=NCOUNT	
	    READ(RECLIN,36)NULL,DLIST(2:4),XX
	    CALL SAVD(LUN3,DLIST,NLIST,IRTFLG)
	    GOTO 70

	ELSEIF (RECLIN(1:4) .EQ. 'ATOM') THEN
	    NATOM=NATOM+1
	    GOTO 70

	ELSEIF (RECLIN(1:3) .EQ. 'TER') THEN
	    NTERM=NTERM+1
	    GOTO 70

	ELSEIF (RECLIN(1:3) .EQ. 'HET') THEN
	    NHET=NHET+1
	    GOTO 70

        ELSEIF (RECLIN(1:3) .EQ. 'END') THEN
            CLOSE(LUN2)
            NCOUNT=NCOUNT+1
            DLIST(1)=NCOUNT
            DLIST(2)=NATOM
            DLIST(3)=NTERM
            DLIST(4)=NHET
            CALL SAVD(LUN3,DLIST,NLIST,IRTFLG)
            CALL SAVDC
            CLOSE(LUN3)
            XNATOM=REAL(NATOM)
            XNTERM=REAL(NTERM)
            XNHET=REAL(NHET)
            R_TMP(4)=XNATOM
            R_TMP(5)=xNTERM
            R_TMP(6)=XNHET
            CALL REG_SET_NSELA(NRG,R_TMP,.FALSE.,IRTFLG)
            RETURN
	ELSE
            GOTO 70 
	ENDIF

	END

C **********************************************************************

	SUBROUTINE PDBRT3

	INCLUDE 'CMBLOCK.INC'
	INCLUDE 'CMLIMIT.INC'
        CHARACTER (LEN=MAXNAM) :: PDBFILE,TRFILE
        CHARACTER *80  RECLIN,HEAD*10, NULL*1,ATOM*4,RESIDUE*3
	LOGICAL        EX,FLAGCELLDIM,FLAGATOM,ZNUM,COMBIND
	DOUBLE PRECISION RM(3,3)

	DATA  LUN2,LUN3/25,26/

	NATOM  = 1
	LENREC = 0
	CALL OPAUXFILE(.TRUE.,PDBFILE,NULL,LUN2,LENREC,'O',
     &                 'PDB INPUT',.TRUE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

	CALL OPAUXFILE(.TRUE.,PDBFILE,NULL,LUN3,LENREC,'N',
     &                 'PDB OUTPUT',.TRUE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

   	CALL RDPRM3S(PHI, THETA, PSI,NOT_USED,
     &              'ROTATION ANGLES PHI, THETA, PSI: ',IRTFLG)

	CALL BLDR(RM,PSI,THETA,PHI)

70	READ(LUN2,30) RECLIN
30	FORMAT(A80)

	IF (RECLIN(1:4) .NE. 'ATOM' .AND. 
     &          RECLIN(1:4) .NE. 'END' .AND.
     & 	        RECLIN(1:4) .NE. 'TER' ) THEN
	   WRITE(LUN3,30) RECLIN
	   GOTO 70 

	ELSEIF (RECLIN(1:4) .EQ. 'ATOM') THEN
	   READ(RECLIN,75) HEAD,N,ATOM,NULL,RESIDUE,NULL,NR2,
     &                     NULL,XO,YO,ZO,
     &                     OCCUPANCY,TEMPERATURE,N
75	   FORMAT(A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,1X,I3)	
C	   WRITE(*,*) XO,YO,ZO

c          CHANGE COORDINATE SYSTEM
	   W  = XO
	   XO = YO
	   YO = W
	   ZO = -ZO

C          AFTER ROTATION CHANGE IT BACK, THAT'S WHY THE ORDER OF X,Y 
C          IS CHANGED AND Z HAS INVERTED SIGN
	   Y=RM(1,1)*XO+RM(1,2)*YO+RM(1,3)*ZO
	   X=RM(2,1)*XO+RM(2,2)*YO+RM(2,3)*ZO
	   Z=-(RM(3,1)*XO+RM(3,2)*YO+RM(3,3)*ZO)
	   WRITE(RECLIN(7:11),80) NATOM
80	   FORMAT(I5)
	   NATOM=NATOM+1
	   WRITE(RECLIN(31:54),85)  X,Y,Z
85	   FORMAT(3F8.3)
	   WRITE(LUN3,30)RECLIN
	   GOTO 70

	ELSEIF (RECLIN(1:3) .EQ. 'TER') THEN
	   WRITE(RECLIN(7:11),80)NATOM
	   WRITE(LUN3,30)RECLIN
	   NATOM=NATOM+1
	   GOTO 70

        ELSEIF (RECLIN(1:3) .EQ. 'END') THEN
	   WRITE(LUN3,30)RECLIN
	   GOTO 999

	ELSE
	   GOTO 70
	ENDIF

999	CLOSE(LUN2)
	CLOSE(LUN3)
	END