!     subroutines for ca potential
!     init : initization procedure
!     ca_energy : calculation of energy
!     orginally coded by XWGG
!     modified by Lumy Dec,11,2005
!     latest modification (May,12,2006)

      subroutine init(dir)
      use energy_function
      use common_sense
      use orientation
      implicit none
      character*255 dir,tempstr
      integer i,j,k,filter2,ii,jj,kk,l,ll,mm,nn,num
      real*8 temp1
      integer eof
      character*255 wholeline

      filter2=index(dir," ")-1

      call block_para

      open(unit=1,file=trim(dir)//'cluster_index.dat',status='old')
      do i=1,num_type
         do j=1,num_type
            do k=1,u_class_num(i)
               do l=1,u_class_num(j)
                  read(1,"(6I6)")ii,jj,kk,ll,u_clu_ind(k,l,i,j) 
     c                 ,u_clu_orient(k,l,i,j)
               enddo
            enddo
         enddo
      enddo
  
      do i=1,num_type
         do k=1,u_class_num(i)
            do l=k,u_class_num(i)
               read(1,"(5I6)")ii,kk,ll,u_clu_ind_sym(k,l,i)
     c          ,u_clu_orient_sym(k,l,i)
            enddo
         enddo
      enddo

      read(1,*)wholeline
      
      do i=1,num_type
         do j=1,num_type
            read(1,"(3I6)")ii,jj,u_clu_num(i,j)
         enddo
      enddo
  
      do i=1,num_type
         read(1,"(2I6)")ii,u_clu_num_sym(i)
      enddo
      
      close(1)
  
!      open(unit=1,file=trim(dir)//
!     c     'out_clu_num.dat',status='unknown')
!      do i=1,num_block
!         do j=i,num_block
!            read(1,"(3I6)")ii,jj,uclu(i,j)
!         enddo
!      enddo
!      close(1)

!      open(unit=1,file=trim(dir)//
!     c     'reduced_index.dat',status='unknown')
!      eof=0
!      do
!         read(1,"(6I6)",IOSTAT=eof)ii,jj,kk,ll,mm,nn
!         if(eof.lt.0)exit
!         u_map(mm,kk,ll,ii,jj)=nn
!         if(ii.eq.jj)then
!            u_map(mm,ll,kk,ii,jj)=nn
!         endif
!      enddo
!      close(1)

      open(unit=1,file=trim(dir)//
     c     'energy_subunit.dat',status='unknown')
      do i=1,num_block
         do j=i,num_block
            if(i.eq.j)then
               num=u_clu_num_sym(block_2_type(i))
            else
               num=u_clu_num(block_2_type(i),block_2_type(j))
            endif 
!            do k=1,uclu(i,j)
            do k=1,num
               read(1,"(3I6,2F15.5)")ii,jj,kk,temp1
     c              ,energy_subunit(k,i,j)
            enddo
         enddo
      enddo
      close(1)

!      open(unit=1,file=trim(dir)//
!     c     'stat_subunit_pair.dat',status='unknown')
!      do i=1,num_block
!         do j=i,num_block
!            read(1,"(2I6,2F10.5)")ii,jj,temp1,energy_u_pair(i,j)
!         enddo
!      enddo
!      close(1)

      open(unit=1,file=trim(dir)//'wei.dat',status='old')
      do i=1,num_block
         do j=i,num_block
            read(1,*)ii,jj,wei(i,j)
         enddo
      enddo
      close(17)

      open(unit=17,file=trim(dir)//'vdw.dat',status='old')
      do i=1,18
         read(17,*)ii,radii(i),welldepth(i)
      enddo
      close(17)

      return
      end subroutine init

      subroutine cal_energy(nres,residues,energy)
      use protein_information
      use energy_function
      use common_sense
      use constant
      implicit none
      integer nres
      type(residue_type) residues(nres)
      real*8 energy

      call energy_subunit_packing(nres,residues)

      call van_de_walls_energy(nres,residues)

      energy=subunit_energy+r_lj*lj_energy
      return
      end subroutine cal_energy
c
      subroutine energy_subunit_packing(nres,residues)
      use protein_information
      use common_sense
      use constant
      use energy_function
      use orientation
      implicit none
      integer nres,indexi,i,j,k,l,indexj,count,ii,jj,iii,jjj
      type(residue_type)residues(nres)
      real*8 dis,ri,rj
      real*8,dimension(3)::vx,vy,vz,origin,rn,rc,ro,r1,r2,r3
      real*8,dimension(3)::vectij,v1

      real*8,dimension(3,max_num_single,nres)::u_o
      real*8,dimension(3,3,max_num_single,nres)::u_r
      integer,dimension(max_num_single,nres)::u_type,u_block
!      integer,dimension(max_num_single,nres):: u_natom
!      integer,dimension(9,max_num_single,nres)::u_atom
      integer,dimension(nres) :: u_num
      integer u_class
      logical,dimension(max_num_single,nres) :: u_exist

!      real*8,dimension(3,max_num_single)::u_o1,u_o2
!      real*8,dimension(3,3,max_num_single)::u_r1,u_r2
!      integer,dimension(max_num_single)::u_type1,u_block1,
!     c     u_type2,u_block2
!      integer u_natom1,u_natom2,u_num1,u_num2
      integer,dimension(9)::u_atom1,u_atom2
      integer u_natom1,u_natom2
!      integer u_class
      real*8,parameter::rcut=(5.)
      logical ncon,match(nres)
!      logical,dimension(max_num_single) :: u_exist1,u_exist2
      integer n_count      

      subunit_energy=0.

      do i=1,nres
         indexi=residues(i)%res_type
         match(i)=.true.
         if(indexi.le.0)then
            match(i)=.false.
         else
            if((i.ne.1).and.(indexi.ne.20).and.(match(i-1)))then
               call add_h(residues(i),residues(i-1))
            endif
            if(indexi.eq.1)call add_cb(residues(i))

            call blocks(residues(i),u_num(i),u_type(:,i),u_block(:,i)
     c           ,u_o(:,:,i),u_r(:,:,:,i),u_exist(:,i))

         endif
      enddo

      do i=1,nres-2
         if(.not.match(i))cycle
         indexi=residues(i)%res_type

         do k=1,u_num(i)
            if(.not.u_exist(k,i))cycle
            call u_atomlist(indexi,u_block(k,i),u_natom1,u_atom1)

            do j=i+2,nres
               if(.not.match(j))cycle
               indexj=residues(j)%res_type

               do l=1,u_num(j)
                  if(.not.u_exist(l,j))cycle
                  call u_atomlist(indexj,u_block(l,j),u_natom2,u_atom2)

!                  if(((u_block(k,i).eq.1).or.(u_block(k,i).eq.7)).and.
!     c                 ((u_block(l,j).eq.1).or.(u_block(l,j).eq.7)))
!     c                 cycle

                  if(((u_block(k,i).eq.1).or.
     c                 ((u_block(k,i).ge.5).and.(u_block(k,i).le.7)))
     c                 .and.((u_block(l,j).eq.1).or.
     c                 ((u_block(l,j).ge.5).and.(u_block(l,j).le.7))))
     c                 cycle
                 
!     check contact start
                  ncon=.false.
                  n_count=0
                  do ii=1,u_natom1
                     iii=u_atom1(ii)
!                     ri=radii(atom_ind(iii,indexi))
                     do jj=1,u_natom2
                        jjj=u_atom2(jj)

!                        rj=radii(atom_ind(jjj,indexj))
!                        rcut=ri+rj+1.
!                        rcut=5.
                        vectij=residues(i)%r(:,iii)-residues(j)%r(:,jjj)
                        dis=sqrt(dot_product(vectij,vectij))
                        if(dis.le.rcut)then
                           ncon=.true.
                           exit
                        endif
                     enddo
                     if(ncon)exit
                  enddo

!     check contact end

!     make statistics
                  if(.not.ncon)cycle

!                  goto 333
!                  if(u_block(k,i).le.u_block(l,j))then
!                     subunit_energy=subunit_energy+
!!     c                    wei(u_block(k,i),u_block(l,j))*
!     c                    energy_u_pair(u_block(k,i),u_block(l,j))
!                  else
!                     subunit_energy=subunit_energy+
!!     c                    wei(u_block(l,j),u_block(k,i))*
!     c                    energy_u_pair(u_block(l,j),u_block(k,i))
!                  endif
!                  cycle

 333              if(u_block(k,i).le.u_block(l,j))then

                     call u_classify(u_type(k,i),u_block(k,i),
     c                    u_o(:,k,i),u_r(:,:,k,i),u_type(l,j),
     c                    u_block(l,j),u_o(:,l,j),u_r(:,:,l,j),u_class)

                     subunit_energy=subunit_energy+
     c                    wei(u_block(k,i),u_block(l,j))*
     c                    energy_subunit(
     c                    u_class,u_block(k,i),u_block(l,j))
!                     subunit_energy=subunit_energy+
!     c                    energy_subunit(
!     c                    u_class,u_block(k,i),u_block(l,j))
!                     write (37,"(5I6,3F10.5)"),i,j,u_block(k,i),
!     c                    u_block(l,j),u_class,
!     c                    wei(u_block(k,i),u_block(l,j)),energy_subunit
!     c                    (u_class,u_block(k,i),u_block(l,j)),
!     c                    wei(u_block(k,i),u_block(l,j))*energy_subunit(
!     c                    u_class,u_block(k,i),u_block(l,j))
                  else
                     
                     call u_classify(u_type(l,j),u_block(l,j),
     c                    u_o(:,l,j),u_r(:,:,l,j),u_type(k,i),
     c                    u_block(k,i),u_o(:,k,i),u_r(:,:,k,i),u_class)

                     subunit_energy=subunit_energy+
     c                    wei(u_block(l,j),u_block(k,i))*
     c                    energy_subunit(
     c                    u_class,u_block(l,j),u_block(k,i))
!                     subunit_energy=subunit_energy+
!     c                    energy_subunit(
!     c                    u_class,u_block(l,j),u_block(k,i))
!                     write (37,"(5I6,3F10.5)"),i,j,
!     c                    u_block(l,j),u_block(k,i),u_class,
!     c                    wei(u_block(l,j),u_block(k,i)),energy_subunit
!     c                    (u_class,u_block(l,j),u_block(k,i)),
!     c                    wei(u_block(l,j),u_block(k,i))*energy_subunit(
!     c                    u_class,u_block(l,j),u_block(k,i))
                  endif

               enddo
            enddo
         enddo
      enddo

      return
      end subroutine energy_subunit_packing

      subroutine van_de_walls_energy(nres,residues)
      use protein_information
      use energy_function
      use common_sense
      implicit none
      integer nres
      type(residue_type)residues(nres)
      integer i,j,k,index_i,index_j,t
      real*8 r2,vectij(3)
      real*8 lj_part
      integer ii,typei,max_atom,max_atom2,iii,jj,jjj,typej,ai,aj
      integer atom_index(14,nres)

      call build_atom_index(nres,residues,atom_index)

      lj_energy=0.
      do ii=1,nres-1
         typei=residues(ii)%res_type
         max_atom=atoms_per_res(typei)
         do iii=1,max_atom

            do jj=ii+1,nres
               typej=residues(jj)%res_type
               max_atom2=atoms_per_res(typej)
               do jjj=1,max_atom2

                  if((iii.le.5).and.(jjj.le.5))cycle
                  
                  if((typej.eq.20).and.(jjj.eq.7).and.(jj.eq.(ii+1))
     c                 .and.(iii.eq.3))cycle

                  vectij=residues(ii)%r(:,iii)-residues(jj)%r(:,jjj)
                  r2=dot_product(vectij,vectij)

                  ai=atom_index(iii,ii)
                  aj=atom_index(jjj,jj)

                  call cal_lj(ai,aj,r2,lj_part)
                  lj_energy=lj_energy+lj_part

               enddo
            enddo
         enddo
      enddo

      return
      end subroutine van_de_walls_energy

!------------ L-J potential-----------
      subroutine cal_lj(ai,aj,r2,lj_part)
      use energy_function
      implicit none
      real*8,parameter :: rc_lj=(6.25D0),r0=(0.5653D0) !1/1.33/1.33
      integer ai,aj
      real*8 dis_lj,r2,rc2,lj_part
      real*8 ir4,ir6,ir12,r1
      real*8 ratio

      real*8,parameter :: rc_atr=(0.7972 ) !1/1.12/1.12

      rc2=dis_lj(ai,aj)

      rc2=rc2*rc2
      r2=r2/rc2
      
!      if(r2.le.rc_lj)then
      if(r2.le.rc_atr)then

         if(r2.gt.r0)then
            ir4=1/r2/r2
            ir6=ir4/r2
            ir12=ir6*ir6
            lj_part=sqrt(welldepth(ai)*welldepth(aj))*
     c           (ir12-2*ir6)
         else
            r1=sqrt(r2)
            lj_part=sqrt(welldepth(ai)*welldepth(aj))*
     c           (49.69D0-40.06D0*r1)
         endif
      else
         lj_part=0.D0
      endif

      return
      end subroutine cal_lj
!
      subroutine build_atom_index(nres,residues,atom_index)
      use protein_information
      use common_sense
      implicit none
      integer nres
      type(residue_type)residues(nres)
      integer atom_index(14,nres),i,j,rtype

!     build atom_index

      atom_index=0
!      scale_long=1.D0
      do i=1,nres
!     backbone atoms
         atom_index(1,i)= 9 !N
         atom_index(2,i)= 1 !CA
         atom_index(3,i)= 2 !C
         atom_index(4,i)=13 !O
      enddo

      do i=1,nres
         rtype=residues(i)%res_type
         select case(rtype)
         case (1)               !GLY
         case (2)               !ALA
            atom_index(5,i)=5 !CB, 3H
         case (3)               !VAL
            atom_index(5,i)=3 !CB, 1H
            atom_index(6,i)=5 !CG1,CG2, 3H
            atom_index(7,i)=5
         case (4)               !ILE
            atom_index(5,i)=3 !CB, 1H
            atom_index(6,i)=4 !CG1, 2H
            atom_index(7,i)=5 !CG2, 3H
            atom_index(8,i)=5 !CD1, 3H
         case (5)               !LEU
            atom_index(5,i)=4 !CB, 2H
            atom_index(6,i)=3 !CG, 1H
            atom_index(7,i)=5 !CD1, 3H
            atom_index(8,i)=5 !CD2, 3H
         case (6)               !SER
            atom_index(5,i)=4 !CB, 2H
            atom_index(6,i)=16!OG, hydroxyl
         case (7)               !THR
            atom_index(5,i)=3 !CB, 1H
            atom_index(6,i)=16!OG1, hydroxyl
            atom_index(7,i)=5 !CG2, 3H
         case (8)               !ASP
            atom_index(5,i)=4 !CB, 2H
            atom_index(6,i)=7 !CG, carboxyl
            atom_index(7,i)=15!OD1, carboxyl
            atom_index(8,i)=15!OD2, carboxyl
         case (9)               !ASN
            atom_index(5,i)=4 !CB, 2H
            atom_index(6,i)=7 !CG, carbonyl
            atom_index(7,i)=14!OD1, carbonyl
            atom_index(8,i)=10!ND2
         case (10)              !GLU
            atom_index(5,i)=4 !CB, 2H
            atom_index(6,i)=4 !CG, 2H
            atom_index(7,i)=7 !CD, carboxyl
            atom_index(8,i)=15!OE1, carboxyl
            atom_index(9,i)=15!OE2, carboxyl
         case (11)              !GLN
            atom_index(5,i)=4 !CB, 2H
            atom_index(6,i)=4 !CG, 2H
            atom_index(7,i)=7 !CD, carbonyl
            atom_index(8,i)=14!OE1, carbonyl
            atom_index(9,i)=10!NE2
         case (12)              !LYS
            atom_index(5,i)=4 !CB, 2H
            atom_index(6,i)=4 !CG, 2H
            atom_index(7,i)=4 !CD, 2H
            atom_index(8,i)=4 !CE, 2H
            atom_index(9,i)=10!NZ
         case (13)              !ARG
            atom_index(5,i)=4 !CB, 2H
            atom_index(6,i)=4 !CG, 2H
            atom_index(7,i)=4 !CD, 2H
            atom_index(8,i)=10!NE
            atom_index(9,i)=7 !CZ, guanidyl
            atom_index(10,i)=10!NH1
            atom_index(11,i)=10!NH2
         case (14)              !CYS
            atom_index(5,i)=8 !CB, 2H
            atom_index(6,i)=17!SG
         case (15)              !MET
            atom_index(5,i)=4 !CB, 2H
            atom_index(6,i)=4 !CG, 2H
            atom_index(7,i)=18!SD
            atom_index(8,i)=5 !CE, 3H
         case (16)              !PHE
            atom_index(5,i)=4 !CB, 2H
            atom_index(6,i)=6 !CG, aromatic ring
            atom_index(7,i)=6 !CD1, aromatic ring
            atom_index(8,i)=6 !CD2, aromatic ring
            atom_index(9,i)=6 !CE1, aromatic ring
            atom_index(10,i)=6 !CE2, aromatic ring
            atom_index(11,i)=6 !CZ, aromatic ring
         case (17)              !TYR
            atom_index(5,i)=4  !CB, 2H
            atom_index(6,i)=6  !CG, aromatic ring
            atom_index(7,i)=6  !CD1, aromatic ring
            atom_index(8,i)=6  !CD2, aromatic ring
            atom_index(9,i)=6  !CE1, aromatic ring
            atom_index(10,i)=6 !CE2, aromatic ring
            atom_index(11,i)=6 !CZ, aromatic ring
            atom_index(12,i)=16!OH, hydroxyl
         case (18)              !TRP
            atom_index(5,i)=4 !CB, 2H
            atom_index(6,i)=6 !CG, aromatic ring
            atom_index(7,i)=6 !CD1, aromatic ring
            atom_index(8,i)=6 !CD2, aromatic ring
            atom_index(9,i)=10!NE1
            atom_index(10,i)=6 !CE2, aromatic ring
            atom_index(11,i)=6 !CE3, aromatic ring
            atom_index(12,i)=6 !CZ2, aromatic ring
            atom_index(13,i)=6 !CZ3, aromatic ring
            atom_index(14,i)=6   !CH2, aromatic ring
         case (19)              !HIS
            atom_index(5,i)=4 !CB, 2H
            atom_index(6,i)=6 !CG, aromatic ring
            atom_index(7,i)=10!ND1
            atom_index(8,i)=6 !CD2, aromatic ring
            atom_index(9,i)=6 !CE1, aromatic ring
            atom_index(10,i)=10!NE2
         case (20)              !PRO
            atom_index(5,i)=4 !CB, 2H
            atom_index(6,i)=4 !CG, 2H
            atom_index(7,i)=4 !CD, 2H
            atom_index(1,i)=12 !N, PRO
         end select
         
      enddo

      return
      end subroutine build_atom_index
!
      real*8 function dis_lj(atomi,atomj)
      use energy_function
      use constant
      implicit none
      integer atomi,atomj
      logical acci,accj,doni,donj
      integer ai,aj

      if((atomi.eq.17).and.(atomj.eq.17))then
         dis_lj=2.0D0
      elseif((atomi.eq.8).and.(atomj.eq.17))then
         dis_lj=3.0D0
      elseif((atomi.eq.17).and.(atomj.eq.8))then
         dis_lj=3.0D0
      else
         if(atomi.gt.atomj)then
            ai=atomj
            aj=atomi
         else
            ai=atomi
            aj=atomj
         endif
         call hbond_type(ai,aj,dis_lj)
         if(dis_lj.lt.small)then
            dis_lj=radii(atomi)+radii(atomj)
         endif
      endif

      return
      end function dis_lj

      subroutine hbond_type(atomi,atomj,dis_hbond)
      implicit none
      integer atomi,atomj
      real*8 dis_hbond

      dis_hbond=0.D0

      select case(atomi)
      case(9:11)
         select case(atomj)
         case(11,12)
            dis_hbond=3.1D0
         case(13,14)
            dis_hbond=2.9D0
         case(15,16)
            dis_hbond=2.9D0
         end select
      case(12)
         if((atomj.eq.15).or.(atomj.eq.16))then
            dis_hbond=2.9D0
         endif
      case(13)
         if((atomj.eq.15).or.(atomj.eq.16))then
            dis_hbond=2.8D0
         endif
      case(15,16)
         if(atomj.eq.16)then
            dis_hbond=2.8D0
         endif
      end select

      return
      end subroutine hbond_type

