!! Read a list of matchs from an alignment (may be non-sequential)
!! Get the ovrlay matrix based on those atoms and the superimpose.
!! Use CAs only.
!! Debugged and ready to go. Uses an alignment that is based on
!! pairs of resSeq numbers, target template. (No other data in alignment file)
!!  CB ---------------   Fri Sep 27 12:16:26 EST 2002
!! 24-APR-00
!! Canabalized FRAGMENT.f
!! This program overlays two fragments.
!!
!! C.Bystroff  8-APR-97
!!
!!
program ovrlayali
        implicit none
        integer*4 I,J,K,L,ndev,jres,kres,ires,nres,jarg,ilocal,lres
        integer*4 MAXRES,nlocal,ios
        parameter(MAXRES=5000)
        real*4 x,y,rms,rmsL,sumdev2,calpha_real(3,MAXRES)
        real*4 calpha_pred(3,MAXRES),lastrmsL,cf,maxrms
        real*4 match_tar(3,MAXRES),match_tpl(3,MAXRES)
        character*80 aline,pdbfile1,pdbfile2,arg1,outfile,alifile
        integer*4,parameter :: MAXLOCAL=40,MINLOCAL=6,LOCALSTEP=1
        real*4 :: confid(MAXRES), minconfid(MAXRES)
        integer*4 :: nseg,ntrue,msize,istart,iend,tstart
        integer*4, parameter :: MAXLIST=4000
        integer*4,dimension(MAXLIST) :: fragstart,fragend
        real*4,dimension(MAXLIST) :: rmslist
        character :: predchar*1,chainID*1,tchainID, alt
        character*3,dimension(MAXRES) :: resname
        integer :: resseq(MAXRES)
        integer*4 :: ilist,mlist,ounit,nfrag,nch,iunit,maxL,junit
        real*4 :: avgrms,btar(1000),btpl(1000)
        real*4,dimension(3) :: vec,Ivec=(/0.,0.,0./)
        real*4,dimension(3,3) :: mtx,Imtx=reshape((/1.,0.,0.,0.,1.,0.,0.,0.,1./),(/3,3/))
        logical :: writing=.true.
        integer(kind=2) :: nmatch, match(2,MAXRES)
!! funtions
        integer iargc
        real getmda
!! initialization
        predchar = 'X'
        alt = '?'
        mlist = MAXLIST
        calpha_pred(1,:) = 999.
        calpha_real(1,:) = 999.
!!
        jarg = iargc()
        !do i=1, jarg
        !  call getarg(i,alifile)
        !  write(*,*)"arg[",i,"]",alifile
        !enddo
        if (jarg.lt.6) then
          write(*,*)   'Usage: xovrlayali alignment targetpdb chain I templatepdb chain I J ', & 
          '[outputfile] '
          write(*,*) 'Outputfile is optional'
          stop 'ovrlayali.f90 v. 17-NOV-04'
        endif
        write(*,*) 'ovrlayali.f90 v. 17-NOV-04 by Xin'
        call getarg(1,alifile)
        call getarg(2,pdbfile1)
        call getarg(3,tchainID)
        if (tchainID == '_' ) tchainID = ' '
        call getarg(4,arg1)
        read(arg1,*) tstart
        call getarg(5,pdbfile2)
        call getarg(6,chainID)
        if (chainID == '_' ) chainID = ' '
        call getarg(7,arg1)
        read(arg1,*) istart
        call getarg(8,arg1)
        read(arg1,*) iend
        if (jarg.ge.9) then
          call getarg(9,outfile)
          writing = .true.
        else
          writing = .false.
        endif
        write(*,'("============= OVRLAYALI ==============")')
        write(*,'(" Superimpose fragments: "," Target  : ",a30," chain=",a1,2x,i5)'),& 
          trim(pdbfile1),tchainID,tstart
        write(*,'("                          Template: ",a30," chain=",a1,2x,i5,i5)'),& 
          trim(pdbfile2),chainID,istart,iend
        write(*,'('' Output file ========> '',a)') outfile(1:index(outfile,'      ')-1)
!! Read alignment file
        open(22,file=alifile,status='old',form='formatted',iostat=ios)
        if (ios/=0) stop 'Error opening alifile'
        i = 0
        do
          i = i + 1
          read(22,*,iostat=ios) match(1:2,i)
          if (ios/=0) exit
          if (match(1,i)==0.or.match(2,i)==0) i = i - 1
        enddo
        nmatch = i -1
        write(*,'(i9," matches read from ",a)') nmatch, trim(alifile)
        close(22)
!! read target structure. PDBfile1 can be .casp file
        open(10,file=pdbfile1,form='formatted',status='old')
        ires = 0
        jres = 0
        lres = -99
        alt = '?'
!! read target ATOM coords
        do  
          read(10,'(a)',iostat=ios) aline
          if(ios/=0) exit
          if (aline(1:5).ne.'ATOM ') cycle 
          if (aline(13:16).ne.' CA ') cycle 
          if (tchainID.eq.'.') tchainID = aline(22:22)
          if (aline(22:22).ne.tchainID)  cycle
          if ((alt.eq.'?').and.(aline(17:17).ne.' ')) alt = aline(17:17)  !! if altLoc is not a space, use it
          if ((aline(17:17).ne.' ').and.(aline(17:17).ne.alt)) cycle 
          read(aline(23:26),*) ires
          if(ires > MAXRES) stop "Error in reading the ires"
          if (ires < tstart) cycle 
          if (ires > tstart+iend-istart)  cycle
          if (ires==lres) cycle 
          lres = ires
          !write(*,*) "in reading targ.pdb: res= ", lres
          jres = jres + 1
          !! THIS VERSION USES ires AS THE INDEX (numbering same as PDB file)
          read(aline(31:54),'(3f8.3)') calpha_pred(1:3,ires)
        enddo
        close(10) 
        nres = jres
        write(*,'(i5,a,a,a,i5)') jres,'residues read from ',&
           trim(pdbfile1),' last res: ',lres
!! read template file, output chain " "
        kres = 0;lres = -99
        open(10,file=pdbfile2,form='formatted',status='old')
        alt = '?'
        do  
          read(10,'(a)',iostat=ios) aline
          if(ios/=0) exit
          if (aline(1:5).ne.'ATOM ') cycle 
          if (aline(13:16).ne.' CA ') cycle 
          if (chainID.eq.'.') chainID = aline(22:22)
          if (aline(22:22).ne.chainID) cycle 
          if ((alt.eq.'?').and.(aline(17:17).ne.' ')) alt = aline(17:17)  !! if altLoc is not a space, use it
          if ((aline(17:17).ne.' ').and.(aline(17:17).ne.alt)) cycle 
          read(aline(23:26),*) ires
          !write(*,*) "in read tmpl, ires=",ires
          if (ires < istart) cycle 
          if (ires > iend) cycle 
          if (ires==lres) cycle 
          lres = ires
          kres = kres + 1
          !write(*,*) "in reading tmpl.pdb: kres= ", kres
          resname(kres) = aline(18:20)
          resseq(kres) = ires
          !!read(aline(31:54),'(3f8.3)') calpha_real(1:3,kres)
          !! THIS VERSION USES ires AS THE INDEX (numbering same as PDB file)
          read(aline(31:54),'(3f8.3)') calpha_real(1:3,ires)
        enddo
        close(10)
        write(*,'(i5,a,a,a,i5)') kres,'residues read from ',&
           trim(pdbfile2),' last res: ',lres
        nres = min(jres,kres)
        !! iend = nres
        write(*,'('' Using residues labels '',i3,'' to '',i4)') resseq(1),&
          resseq(nres)
        write(*,'(i9,'' C-alpha atoms total. '')') nres
!! save only aligned residues (match)
        btpl(:) = 10.0
        btar(:) = 1.0
        do i=1,nmatch
          match_tar(1:3,i) = calpha_pred(1:3,match(1,i))
          match_tpl(1:3,i) = calpha_real(1:3,match(2,i))
          !! diagnostic    
          !write(*,'(i7,3f8.3)') match(1,i),match_tar(1:3,i),match(2,i),match_tpl(1:3,i)
          btpl(match(2,i)) = 42.
          btar(match(1,i)) = 50.
        enddo
!! Superimpose each fragment on target, write out coords
        call superimp(match_tpl,match_tar,rms,nmatch,mtx,vec)
        if (writing) then
          junit = 9
          open(junit,file=pdbfile2,form='formatted',status='old')
          ounit = 20
          open(ounit,file=outfile,form='formatted',status='replace')
          rewind(ounit)  ! overwrite
!! write template fragment, no rotation
          I = istart
          !! J = istart + nres - 1
          J = iend
          call writepdb(junit,I,J,Imtx,Ivec,btpl,rms,chainID,chainID,ounit)
          write(ounit,'(a)') 'TER   '
          close(junit)
!! write target fragment
          iunit = 10
          open(iunit,file=pdbfile1,form='formatted',status='old')
          I = tstart
          !! J = tstart + nres - 1
          J = tstart + (iend - istart)
          chainID = ' '
          call writepdb(iunit,I,J,mtx,vec,btar,rms,tchainID,predchar,ounit)
          write(ounit,'(a)') 'TER   '
!! Done.
          close(ounit)
          close(iunit)
        endif
        write(*,'('' RMSD(CA) ='',f8.3)') rms
!!
        open(11,file="matvec.dat",status='replace',form='formatted')
        write(11,'(3f12.5)') mtx, vec
        close(11)
        if (writing) stop 'ovrlayali.f v.27-SEP-02'
end program ovrlayali
!!
!!
        subroutine writepdb(iunit,istart,iend,mtx,vec,bfac,x,chainID,chain,ounit)
        implicit none
        integer*4 :: N,iunit,ounit,I,J,K,ires,istart,iend
        real*4,dimension(3) :: coord
        real*4 :: bfac(1000),x
        real*4,dimension(3) :: vec
        real*4,dimension(3,3) :: mtx
        character*1 :: chainID,chain,alt
        character*80 :: aline
!!
        rewind(iunit)
        alt = '?'
 5      read(iunit,'(a)',end=1000) aline
        if (aline(1:6) .ne. 'ATOM  ') goto 5
        if (aline(22:22) .ne. chainID) goto 5
        if ((alt.eq.'?').and.(aline(17:17).ne.' ')) alt = aline(17:17)  !! if altLoc is not a space, use it
        if ((aline(17:17).ne.' ').and.(aline(17:17).ne.alt)) goto 5
        read(aline(23:26),*) ires
        if (ires <= iend .and. ires >= istart) then
          read(aline(31:54),'(3f8.3)') coord
          call move(coord,mtx,vec)
          write(aline(31:66),'(3f8.3,2f6.2)') coord,x,bfac(ires)
          aline(22:22) = chain
          write(ounit,'(a80)') aline
        endif
        goto 5
 1000   return
        end
!!
        subroutine move(x,mtx,vec)
        real*4,dimension(3) :: x,vec,y
        real*4,dimension(3,3) :: mtx
        integer*4 :: I,J
        do J=1,3
          y(J) = 0.
          do I=1,3
            y(J) = y(J) + x(I)*mtx(J,I)
          enddo
        enddo
        do I=1,3
          x(I) = y(I) + vec(I)
        enddo
        return
        end
!!************************************************************
!************************************************************************
        subroutine superimp(coord0,coord1,err,natm,mtx,vec)
        !!  06-SEP-02
!** coord0 is template.
!** coord1 is target, returned
!** err is rmsd of natm atoms
!** get the rms deviation between two arbitrarilly oriented
!** sets of 'natm' atoms.
        implicit none
        integer(kind=2),intent(in) ::  natm
        real,intent(out) :: err
        real,dimension(3,natm),intent(in):: coord0
        real,dimension(3,natm),intent(inout):: coord1
        real,intent(out) ::  mtx(3,3),vec(3)
        real(kind=4) :: aa(3,3),rot(3,3),t(3),sig,sg,gam,bet,cc,bb
        real,dimension(:,:),allocatable :: x0,x1
        real ::  mat(3,3),xc0(3),xc1(3)
        integer :: ix,iy,iz,i,j,k,iflag,ict,nnn
        real :: y
        real,parameter :: tol=0.0001
!**
        allocate(x0(3,natm),x1(3,natm),stat=i)
        !! diagnostic
        !! write(*,'("Number of atoms to allocate=",i9)') natm
        if (i/=0) then
          stop 'Failed to allocate space in SUPERIMP'
        endif
!* note: if any atom is missing, omit it from the calculation
!*       if less than 5 atoms left, return err=-1 as a flag
        nnn = 0
        xc0 = 0.;  xc1 = 0.
        do I=1,natm
          if (coord0(1,I) > 998.0) cycle
          if (coord1(1,I) > 998.0) cycle
          nnn = nnn + 1
          x0(1:3,nnn) = coord0(1:3,I)
          x1(1:3,nnn) = coord1(1:3,I)
          xc0 = xc0 + x0(1:3,nnn)
          xc1 = xc1 + x1(1:3,nnn)
        enddo
        if (nnn < 5) then
          write(*,'("Number of template atoms   =",i9,"  <===too small to superimpose.")') nnn
          err = -1
          return
        endif
        xc0 = xc0/real(nnn)   !! center of mass, target
        xc1 = xc1/real(nnn)   !! center of mass, template
        do I=1,nnn
          x0(:,I) = x0(:,I) - xc0
          x1(:,I) = x1(:,I) - xc1
        enddo
!**
        aa = 0.
        do i=1,3
          do j=1,3
            aa(i,j) = dot_product(x1(i,:),x0(j,:))
          enddo
        enddo
        !  do k=1,nnn
        !    do i=1,3
        !      do j=1,3
        !        aa(i,j)=aa(i,j)+x1(i,k)*x0(j,k)
        !      enddo
        !    enddo
        !  enddo
        rot = 0.
        do i=1,3
          rot(i,i) = 1.
        enddo
!! new f90 version
        ict=0
        iflag=0
        ix=1
        do while (ict < 1000)
          ict=ict+1
          do ix=1,3
            iy=mod(ix,3)+1
            iz=6-ix-iy
            sig=aa(iz,iy)-aa(iy,iz)
            gam=aa(iy,iy)+aa(iz,iz)
            sg=sqrt(sig*sig+gam*gam)
            if(sg == 0.) cycle
            sg=1./sg
            if(abs(sig) <= tol*abs(gam)) cycle
            do k=1,3
              bb=gam*aa(iy,k)+sig*aa(iz,k)
              cc=gam*aa(iz,k)-sig*aa(iy,k)
              aa(iy,k)=bb*sg
              aa(iz,k)=cc*sg
              bb=gam*rot(iy,k)+sig*rot(iz,k)
              cc=gam*rot(iz,k)-sig*rot(iy,k)
              rot(iy,k)=bb*sg
              rot(iz,k)=cc*sg
            enddo
            iflag=1
          enddo
          if (iflag==0) exit
        enddo
        err=0.
        do I=1,natm
          if (coord0(1,I) > 998.0) cycle
          if (coord1(1,I) > 998.0) cycle
          x1(:,I) = coord1(:,I) - xc1
          do j=1,3
            t(j)=dot_product(rot(j,:),x1(:,I))
          enddo
          !! coord1(:,I)= t + xc0  !! coord1 is now superimposed on coord0
          t(:)= t + xc0  !! t is now superimposed on coord0
          t = coord0(:,I) - t(:)
          err = err + dot_product(t,t)
        enddo
        err=sqrt(err/nnn)
        deallocate(x0,x1,stat=i)
        do i=1,3
          t(i)=dot_product(rot(i,:),xc1(:))
          do j=1,3
            mtx(i,j)=rot(i,j)
          enddo
        enddo
        vec = xc0 - t
        !! diagnostic
        do I=1,natm
          if (coord0(1,I) > 998.0) cycle
          if (coord1(1,I) > 998.0) cycle
          call move(coord1(1,I),mtx,vec)
          !write(*,'("Superimposed ATOMS ",i5,6f8.3)') I,coord0(1:3,I),coord1(1:3,I)
        enddo
        end subroutine superimp
!************************************************************************
