c Read in coordinates and rotation angles around x and z axes.
c Peform rotation first around x, and then around z.
c Michael R. Sawaya August 19, 2022
       character pdbinfile1*100
c      Output PDB file after rotation
       character pdboutfile1*100 
c      lineinput is a temporary storage of pdb input lines to discriminate atoms from header
       character lineinput*80
c      Chain id of contacting molecules
       character ch*1
c      Resid of contacting molecules
       integer resid
c      Number of atoms read from PDB file.
       integer n
c      Counting how many chains there are in input PDB.
       character prevchain*1,chainlist(62)*1
       integer nchains
c      PDB items to read and parse.
       real x0,y0,z0
       real x1,y1,z1
       real x2,y2,z2
       character resn*3,atnam*5
c      rotation angles around x and z axes.
       real rotxdeg, rotzdeg
       real rotxrad, rotzrad
c      Matrix elements of operation
       real r1c1(2),r1c2(2),r1c3(2)
       real r2c1(2),r2c2(2),r2c3(2)
       real r3c1(2),r3c2(2),r3c3(2)
c      Solvation Free Energy for an atom
       real sfe(9000)
c
        write(*,*) 'Please type name of input .pdb file'
        read(*,'(a100)')pdbinfile1
        open (1,file=pdbinfile1, status='old', form='formatted')
        print*, pdbinfile1
        print*, 'Type value for rotation around x axis in degrees'
        read (*,"(F6.1)")rotxdeg
        write(6,101)rotxdeg
101     format('You said rotate',f6.1,' degrees around x')
        print*, 'Type value for rotation around z axis in degrees'
        read (*,"(F6.1)")rotzdeg
        write(6,102)rotzdeg
102     format('You said rotate',f6.1,' degrees around z')
        write(*,*) 'Please type name of output .pdb file'
        read(*,'(a100)')pdboutfile1
 	open (2,file=pdboutfile1,status='new',form='formatted')
c       Formulate matrix for rotation around Z axis.
        pi=3.141592653589793238
        rotxrad=rotxdeg*pi/180.0
        r1c1(1)=1
        r1c2(1)=0
        r1c3(1)=0
        r2c1(1)= 0
        r2c2(1)= cos(rotxrad)
        r2c3(1)=-sin(rotxrad)
        r3c1(1)=0
        r3c2(1)= sin(rotxrad)
        r3c3(1)= cos(rotxrad)
c
        rotzrad=rotzdeg*pi/180.0
        r1c1(2)= cos(rotzrad)
        r1c2(2)=-sin(rotzrad)
        r1c3(2)=0
        r2c1(2)= sin(rotzrad)
        r2c2(2)= cos(rotzrad)
        r2c3(2)=0
        r3c1(2)=0
        r3c2(2)=0
        r3c3(2)=1
c       Now I am reading in the coordinates from pdbinfile1.
5       n=1
2       format(a80)
4       read (1,2,end=6)lineinput
c       Header gets written to pdboutfile1 (unit 2).
        if((lineinput(1:6).ne.'HETATM').and.
     1     (lineinput(1:4).ne.'ATOM'))then
         write(2,2)lineinput
         goto 4
        endif
c       All atoms get transformed by rotxdeg and rotzdeg  and then written to pdb file.
        if((lineinput(1:6).eq.'HETATM').or.
     1     (lineinput(1:4).eq.'ATOM'))then
         atnam=lineinput(13:17)
         resn=lineinput(18:20)
         ch=lineinput(22:22)
         read(lineinput(23:26),'(i4)')resid
         read(lineinput(31:38),*)x0
         read(lineinput(39:46),*)y0
         read(lineinput(47:54),*)z0
c        Rotate atom around X axis.
         x1=(x0*r1c1(1))+(y0*r1c2(1))+(z0*r1c3(1))
         y1=(x0*r2c1(1))+(y0*r2c2(1))+(z0*r2c3(1))
         z1=(x0*r3c1(1))+(y0*r3c2(1))+(z0*r3c3(1))
c        Rotate atom around Z axis.
         x2=(x1*r1c1(2))+(y1*r1c2(2))+(z1*r1c3(2))
         y2=(x1*r2c1(2))+(y1*r2c2(2))+(z1*r2c3(2))
         z2=(x1*r3c1(2))+(y1*r3c2(2))+(z1*r3c3(2))
         write(2,103)lineinput(1:30),x2,y2,z2,lineinput(55:80)
103      format(a30,3f8.3,a26)
         n=n+1
        endif
        goto 4
6       write(6,104)n
104     format('Rotated', i6,' atoms.')
        print*,'Adioooooos amigoooooos'
c
 999	stop
	end
