      subroutine import_pdb(nmax,filename,nres,residues)
      use protein_information
      use common_sense
      use constant
      implicit none
      type(residue_type) residues(nmax)
      integer nmax,natom,nres
      integer atomindex
      integer i,j,k,eof,filter,nindex,preindex
      real*8 x0,y0,z0
      character filename*255,head*13,type*3
      character gap*2,res*3,gap2*1,chid*1,other*4,last*26
      character precode*1,alt_now*1,alt*1,icode*1

!      filter=index(filename,".pdb")
      open(unit=15,file=trim(filename),status='old')
      i=0
      j=1
      eof=0
      do
         read(15,2001,IOSTAT=eof) head,type,gap,res,
     c        gap2,chid,nindex,other,x0,y0,z0,last
         if(eof.lt.0)then
            exit
         endif
         if(head(1:6).eq."ENDMDL")exit
         if((head(1:4).eq."ATOM").and.(gap.ne."B"))then
            if((i.eq.0).and.(res(1:3).eq.'ACE'))cycle
            if(type(1:1).eq.'H')cycle

            if((j.eq.1).or.
     c           ((preindex.ne.nindex).or.(icode.ne.precode).or.
     c           (residues(i)%res_type.ne.index(res_list,res)/4+1)))then
!            if((j.eq.1).or.(preindex.ne.nindex))then
               i=i+1
               residues(i)%res_type=index(res_list,res)/4+1
               preindex=nindex
               precode=icode
               alt_now=" "
               do k=1,14
                  residues(i)%r(1,k)=large
               enddo
            endif

            if((alt_now.eq." ").and.(alt.ne." "))then
               alt_now=alt
            endif

            if((i.eq.1).and.(type.eq."NT"))type="N"

            select case(type)
            case ('AE1 ')
               type='OE1 '
            case ('AE2 ')
               type='NE2 '
            case ('AD1 ')
               type='OD1 '
            case ('AD2 ')
               type='ND2 '
            end select

            atomindex=index(atom_list(residues(i)%res_type),type)
            if(atomindex.ne.0)then
               residues(i)%r(1,atomindex/3+1)=x0
               residues(i)%r(2,atomindex/3+1)=y0
               residues(i)%r(3,atomindex/3+1)=z0
               j=j+1
            elseif((type.ne."OXT").and.(type(1:1).ne."H").and.
     c              (type(1:1).ne.'Q').and.(type(1:1).ne.'D').and.
     c              (type.ne."OT"))then
               print *,"unknown atom type in residue ",res,' : ',type
            endif
         endif         
      enddo
      natom=j-1
      nres=i
!     print *,natom,nres
 2001 format(A13,A3,A1,A3,A1,A1,I4,A4,3F8.3,A26)
      close(15)
      return
      end subroutine import_pdb

!---------------------------------------      
      subroutine export_pdb(nres,residues,filename)
      use protein_information
      use common_sense
      use constant
      implicit none
      integer nres,filter,i,j
      type(residue_type) residues(nres)
      character filename*100
      real*8,parameter :: ratio=(1.00),bfactor=(10.00)
      character*7,parameter :: head=('ATOM   ')
      integer atomno,resno,restypeno
      character res*3

      filter=index(filename,".pdb")+3

      open(unit=10,file=filename(1:filter),status='unknown')
      atomno=0
      do i=1,nres
         restypeno=residues(i)%res_type
         res=res_list(restypeno*4-3:restypeno*4-1)
         do j=1,14
            if(residues(i)%r(1,j).ne.large)then
               atomno=atomno+1
               write(10,3001)head,atomno,atom_list(restypeno)(j*3-2:j*3)
     &              ,res,i,residues(i)%r(1,j),residues(i)%r(2,j),
     &              residues(i)%r(3,j),ratio,bfactor
            endif
         enddo
      enddo

 3001 format(A7,I4,2X,A3,X,A3,2X,I4,4X,3F8.3,2F6.2)
      close(10)
!     note: no "OXT" atom at the end, need to build additionally
      return
      end subroutine export_pdb

c---------------------------------------------
      SUBROUTINE crossnew(r1,r2,r3)
C
C     Cross product of two vectors.
C     
C     Aaron R. Dinner
C
      implicit none
      REAL*8,dimension(3)::r1,r2,r3

      r3(1)=r1(2)*r2(3)-r1(3)*r2(2)
      r3(2)=r1(3)*r2(1)-r1(1)*r2(3)
      r3(3)=r1(1)*r2(2)-r1(2)*r2(1)

      RETURN
      END

!-----------------------------------
      real*8 function dihedral(a,b,c,d)
!      use common_sense
      implicit none
      real*8,dimension(3) :: a,b,c,d
      real*8,dimension(3) :: d1,d2,d3,d12,d23
      real*8 sum,sign
c
      d1=b-a
      d2=c-b
      d3=d-c
      call cross_product(d1,d2,d12)
      call cross_product(d2,d3,d23)
      sum=dot_product(d12,d23)
      sign=dot_product(d12,d3)
      if(sign.ge.0)then
         dihedral=dacos(sum)
      else
         dihedral=-dacos(sum)
      endif
      return
      end
c
!-----------------------------------
      subroutine cross_product(a,b,c)
      implicit none
      real*8,dimension(3) :: a,b,c
      real*8 sum
      c(1)=a(2)*b(3)-a(3)*b(2)
      c(2)=a(3)*b(1)-a(1)*b(3)
      c(3)=a(1)*b(2)-a(2)*b(1)
      sum=sqrt(c(1)**2+c(2)**2+c(3)**2)
      c=c/sum
      return
      end
c
!--------------------------------------
      subroutine inv_dihedral(r1,r2,bond,dihedral,r3)
      implicit none
      integer i
      real*8,dimension(3) :: r1,r2,r3,ry,rz,weight
      real*8 bond,dihedral
      r2=r2/sqrt(dot_product(r2,r2))
      call cross_product(r1,r2,rz)
      call cross_product(rz,r2,ry)
      weight(1)=dcos(bond)
      weight(2)=dcos(dihedral)*dsin(bond)
      weight(3)=dsin(dihedral)*dsin(bond)
      do i=1,3
         r3(i)=r2(i)*weight(1)+ry(i)*weight(2)+rz(i)*weight(3)
      enddo
      return
      end subroutine inv_dihedral
c
!-----------------------------------------
      subroutine protein_rmsd(nres,res0,res1,rdiff,method,rot)
!     method: 1: all built atoms;  2: all backbone atoms; 3: CA atoms
      use protein_information
      use constant
      implicit none
      integer nres,i,j,atompair(2,14*nres),natom,count,method,rot,j1,j2
      real*8,dimension(3,14*nres) :: r,r0
      real*8 rdiff,vectij(3)
      type(residue_type) res1(nres),res0(nres)
c
      natom=0
      count=0
      do i=1,nres
         do j=1,14
            if(res1(i)%r(1,j).ne.large)then
               natom=natom+1
               r(:,natom)=res1(i)%r(:,j)
               r0(:,natom)=res0(i)%r(:,j)
               select case(method)
               case(1)
                  count=count+1
                  atompair(1,natom)=natom
                  atompair(2,natom)=natom
               case(2)
                  if(j.le.5)then
                     if((i.eq.nres).and.(j.eq.4))exit
                     count=count+1
                     atompair(1,count)=natom
                     atompair(2,count)=natom
                  endif
               case(3)
                  if(j.eq.2)then
                     count=count+1
                     atompair(1,count)=natom
                     atompair(2,count)=natom
                  endif
               case default
                  print *,'error select rmsd method!'
                  stop
               end select
            endif
         enddo
      enddo
c
      if(rot.eq.1)then
         call ROTLSQ(r0(1,:),r0(2,:),r0(3,:),natom,r(1,:),r(2,:),r(3,:),
     &        natom,atompair,count)
      endif
      rdiff=0.
      do i=1,count
         j1=atompair(1,i)
         j2=atompair(2,i)
         vectij=r(:,j2)-r0(:,j1)
         rdiff=rdiff+dot_product(vectij,vectij)
      enddo
      rdiff=sqrt(rdiff/count)
      if(rot.eq.1)then
         natom=0
         do i=1,nres
            do j=1,14
               if(res1(i)%r(1,j).ne.large)then
                  natom=natom+1
                  res1(i)%r(:,j)=r(:,natom)
               endif
            enddo
         enddo
      endif
      return
      end subroutine protein_rmsd

