module orientation
  implicit none
  integer,private,parameter :: dbl= selected_real_kind(p=13) !double precision

  !constant data:
  !num_type     :   total # of geometric types
  !num_block    :   total # of block types
  !num_dom      :   total # of dominant directions
  !max_class    :   max # of classes based on direction
  !num_rotate   :   max # of possible rotation
  !max_atom     :   max # of atoms in a block type
  !block_2_type :   geometric type of each block type
  !block_2_atom :   # of atoms in each block type
  !u_plane      :   plane or non-plane of each geometric type

  integer,parameter :: num_type = 9 ,num_dom = 27 ,num_block = 19 ,&
       max_class = 17 ,max_rotate = 8 ,max_atom = 9 , max_num_single = 5
  integer,parameter,dimension(num_block) :: block_2_type = (/ &
       1,1,1,1,1,2,3,3,3,3,4,4,4,5,5,6,7,8,9/) ,block_2_natom = (/ &
       1,1,1,1,1,2,2,2,2,2,3,3,4,3,3,6,5,9,5/)
  logical,parameter,dimension(num_type) :: u_plane=(/&
       .false.,.false.,.false.,.true.,.true.,&
       .false.,.true.,.true.,.true./)

  !sharable data :
  !u_class_num   :  # of classes of each block type
  !u_class_sub   :  class indeces
  !u_sym         :  symmetry of each class
  !u_axis_dom    :  unit vector of dominant directions
  !u_axis_1,2    :  first and second axices for each dominant direction of each block type

  integer,save,dimension(num_type) :: u_class_num = 0
  integer,save,dimension(num_dom,num_type):: u_class_sub = 0 ,u_sym = 0
  real(kind=dbl),save,dimension(3,num_dom,num_type) :: &
       u_axis_dom = 0.D0 ,u_axis_1 = 0.D0 ,u_axis_2 = 0.D0 
!
  !u_pos         :  relative atomic positions for each block type
  !u_clu_ind,u_clu_ind_sym : cluster indeces
  !u_clu_orient,u_clu_orient_sym : # of rotation
  !u_clu_num,u_clu_num_sym : # of clusters (classes)
  !wei : weight
  !u_dis_max : max atomic distance between two blocks
  !u_map : map towards reduced clusters
  !u_lj (omit)
  !uclu :# of reduced clusters
  integer,save,dimension(max_class,max_class,num_type,num_type) :: &
       u_clu_ind,u_clu_orient
  integer,save,dimension(max_class,max_class,num_type) :: &
       u_clu_ind_sym,u_clu_orient_sym
  integer,save,dimension(num_type) :: u_clu_num_sym
  integer,save,dimension(num_type,num_type) :: u_clu_num
  integer,save,dimension(max_rotate,max_class,max_class,num_block,num_block) &
       :: u_map
  real(kind=dbl),save,dimension(num_block,num_block) :: wei
  real(kind=dbl),save,dimension(3,max_atom,num_block) :: u_pos
  real(kind=dbl),save,dimension(num_block) :: u_dis_max
!  integer,save,dimension(num_type,num_block) :: u_lj
  integer,save,dimension(num_block,num_block) :: uclu

!contained procedures

  private cross_product,ch_ref,ch_ref_inv,ref_pos,istr,inv_dihedral

contains
  subroutine block_para !initialize parameters
    implicit none
    !u_type : type index
    !ux,uy,uz : unit vectors
    !nx,ny,nz : direction indeces , -1,0,+1
    !dom_dir  : dominant direction value
    !vi : direction vector
    integer :: u_type,nx,ny,nz,dom_dir
    real(kind=dbl),dimension(3):: ux,uy,uz,vi

    do u_type=1,9
       select case(u_type)
       case(1:3)
          u_class_num(u_type)=5
          u_class_sub(:,u_type)=(/1,2,3,4,5,0,0,0,0,0,&
               0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
          u_sym(:,u_type)=(/1,3,3,3,1,0,0,0,0,0,&
               0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
       case(6)
          u_class_num(u_type)=3
          u_class_sub(:,u_type)=(/1,2,3,2,1,0,0,0,0,0,&
               0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
          u_sym(:,u_type)=(/1,3,2,3,1,0,0,0,0,0,&
               0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
       case(4)
          u_class_num(u_type)=11
          u_class_sub(:,u_type)=(/9,10,9,6,11,6,7,8,7,4,&
               5,4,1,0,1,2,3,2,9,10,9,6,11,6,7,8,7/)
          u_sym(:,u_type)=(/4,-1,4,1,-1,1,1,-1,1,-3,&
               2,-3,1,0,1,1,2,1,4,-1,4,1,-1,1,1,-1,1/)!modified
       case(5)
          u_class_num(u_type)=17
          u_class_sub(:,u_type)=(/15,16,14,10,17,9,12,13,11,7,&
               8,6,2,0,1,4,5,3,15,16,14,10,17,9,12,13,11/)
          u_sym(:,u_type)=(/4,4,4,1,4,1,1,4,1,-3,&
               -3,-3,1,0,1,1,-3,1,4,4,4,1,4,1,1,4,1/)
       case(7)
          u_class_num(u_type)=17
          u_class_sub(:,u_type)=(/15,16,14,10,17,9,12,13,11,7,&
               8,6,2,0,1,4,5,3,15,16,14,10,17,9,12,13,11/)
          u_sym(:,u_type)=(/4,4,4,4,4,4,4,4,4,-3,&
               -3,-3,-3,0,-3,-3,-3,-3,4,4,4,4,4,4,4,4,4/)
       case(8)
          u_class_num(u_type)=15
          u_class_sub(:,u_type)=(/13,14,12,9,15,8,10,11,10,6,&
               7,5,2,0,1,3,4,3,13,14,12,9,15,8,10,11,10/)
          u_sym(:,u_type)=(/4,4,4,4,4,4,4,-1,4,-3,&
               -3,-3,-3,0,-3,-3,2,-3,4,4,4,4,4,4,4,-1,4/)
       case(9)
          u_class_num(u_type)=11
          u_class_sub(:,u_type)=(/9,10,9,6,11,6,7,8,7,4,&
               5,4,1,0,1,2,3,2,9,10,9,6,11,6,7,8,7/)
          u_sym(:,u_type)=(/4,-1,4,4,-1,4,4,-1,4,-3,&
               2,-3,-3,0,-3,-3,2,-3,4,-1,4,4,-1,4,4,-1,4/)
       end select

       ux=(/1.,0.,0./)
       uy=(/0.,1.,0./)
       uz=(/0.,0.,1./)

       if(u_plane(u_type))then !for plane
          do nx=-1,1
             do ny=-1,1
                do nz=-1,1
                   if((nx.eq.0).and.(ny.eq.0).and.(nz.eq.0))cycle
                   dom_dir=nx+ny*3+nz*9+14
                   vi=nx*ux+ny*uy+nz*uz
                   u_axis_dom(:,dom_dir,u_type)=vi/sqrt(dot_product(vi,vi))

                   select case(dom_dir)
                   case(13,15)   !x
                      if(u_sym(dom_dir,u_type).eq.1)then
                         u_axis_1(:,dom_dir,u_type)=0.
                         u_axis_2(:,dom_dir,u_type)=0.
                      elseif(u_sym(dom_dir,u_type).eq.-3)then
                         u_axis_1(:,dom_dir,u_type)=uz
                         u_axis_2(:,dom_dir,u_type)=uy
                      else
                         u_axis_1(:,dom_dir,u_type)=uy
                         u_axis_2(:,dom_dir,u_type)=uz
                      endif
                   case(11,17)   !y
                      if(u_sym(dom_dir,u_type).eq.1)then
                         u_axis_1(:,dom_dir,u_type)=0.
                         u_axis_2(:,dom_dir,u_type)=0.
                      elseif(u_sym(dom_dir,u_type).eq.-1)then
                         u_axis_1(:,dom_dir,u_type)=ux
                         u_axis_2(:,dom_dir,u_type)=uz
                      else
                         u_axis_1(:,dom_dir,u_type)=uz
                         u_axis_2(:,dom_dir,u_type)=ux
                      endif
                   case(5,23)    !z
                      if(u_sym(dom_dir,u_type).eq.1)then
                         u_axis_1(:,dom_dir,u_type)=0.
                         u_axis_2(:,dom_dir,u_type)=0.
                      elseif(u_sym(dom_dir,u_type).eq.-2)then
                         u_axis_1(:,dom_dir,u_type)=uy
                         u_axis_2(:,dom_dir,u_type)=ux
                      else
                         u_axis_1(:,dom_dir,u_type)=ux
                         u_axis_2(:,dom_dir,u_type)=uy
                      endif
                   case(1,3,7,9,19,21,25,27) !xyz
                      if(u_sym(dom_dir,u_type).eq.1)then
                         u_axis_1(:,dom_dir,u_type)=0.
                         u_axis_2(:,dom_dir,u_type)=0.
                      else
                         u_axis_1(:,dom_dir,u_type) &
                              =(2*nz*uz-nx*ux-ny*uy)/sqrt(6.)
                         u_axis_2(:,dom_dir,u_type)= &
                             (ny*uy-nx*ux)/sqrt(2.)
                      endif
                   case default !xy ,yz or zx
                      if(u_sym(dom_dir,u_type).eq.1)then
                         u_axis_1(:,dom_dir,u_type)=0.
                         u_axis_2(:,dom_dir,u_type)=0.
                         cycle
                      endif
                      if(nx.eq.0)then                        
                         u_axis_1(:,dom_dir,u_type)=ux
                         u_axis_2(:,dom_dir,u_type)=(ny*uy-nz*uz)/sqrt(2.)
                      elseif(ny.eq.0)then
                         u_axis_1(:,dom_dir,u_type)=uy
                         u_axis_2(:,dom_dir,u_type)=(nz*uz-nx*ux)/sqrt(2.)
                      elseif(nz.eq.0)then
                         u_axis_1(:,dom_dir,u_type)=uz
                         u_axis_2(:,dom_dir,u_type)=(nx*ux-ny*uy)/sqrt(2.)
                      endif                      
                   end select
                enddo
             enddo
          enddo
       endif
    enddo
    
  end subroutine block_para
!
  subroutine test_para !check parameters in subroutine block_para
    implicit none
    integer ::u_type,dom_dir,class,count,u_block
    real(kind=dbl),dimension(3):: d1,d2,u_o
    real(kind=dbl),dimension(3,3)::u_r
    character(len=255) :: letter
    integer :: i,j,ii,jj

    !read block position data
    open(unit=1,file='block_pos.dat',status='old')
    do i=1,num_block
       do j=1,block_2_natom(i)
          read(1,"(2I4,3F10.5,I5)")ii,jj,u_pos(:,j,i)
       enddo
    enddo
    close(1)

    u_o=0.
    do u_block=1,num_block
       letter=istr(u_block)
       open(unit=1,file="test_"//trim(letter)//".pdb",status='replace')
       u_type=block_2_type(u_block)
       do dom_dir=1,num_dom
          if(u_class_sub(dom_dir,u_type).eq.0)cycle
          do i=1,block_2_natom(u_block)
             write(1,"(A6,I5,2X,A3,A1,A3,1X,I5,4X,"// &
                  "3F8.3,A12)")"ATOM   ",i,"C  "," ","BLK" &
                  ,dom_dir,u_pos(:,i,u_block),"  1.00 20.00"
          enddo
          u_r(:,1)=u_axis_1(:,dom_dir,u_type)
          u_r(:,2)=u_axis_2(:,dom_dir,u_type)
          u_r(:,3)=u_axis_dom(:,dom_dir,u_type)
          count=dom_dir-1
          call export_block(count,1,u_type,u_o,u_r)
          write(1,"(A3)")"END"
       enddo
       close(1)
    enddo
  end subroutine test_para
!
  subroutine export_block(count,res_type,u_type,u_o,u_r) !export blocks
    use common_sense      
    implicit none
    integer,intent(inout) :: count
    integer,intent(in) :: res_type,u_type
    real(kind=dbl),intent(in),dimension(3) :: u_o
    real(kind=dbl),intent(in),dimension(3,3) :: u_r

    count=count+1
    write(1,"(A6,I5,2X,A3,A1,A3,1X,I5,4X,3F8.3,A12)") &
         "ATOM   ",1,"C  "," ",res_list(res_type*4-3:res_type*4-1) &
         ,count,u_o,"  1.00 20.00"
    if(u_type.ne.1)then
       write(1,"(A6,I5,2X,A3,A1,A3,1X,I5,4X,3F8.3,A12)") &
            "ATOM   ",2,"N  "," ",res_list(res_type*4-3:res_type*4-1) &
            ,count,u_o+1.5*u_r(:,1),"  1.00 20.00"
       if(u_plane(u_type))then
          write(1,"(A6,I5,2X,A3,A1,A3,1X,I5,4X,3F8.3,A12)") &
               "ATOM   ",3,"O  "," ",res_list(res_type*4-3:res_type*4-1) &
               ,count,u_o+1.5*u_r(:,2),"  1.00 20.00"
          write(1,"(A6,I5,2X,A3,A1,A3,1X,I5,4X, 3F8.3,A12)") &
               "ATOM   ",4,"S  "," ",res_list(res_type*4-3:res_type*4-1) &
               ,count,u_o+1.5*u_r(:,3),"  1.00 20.00"
       endif
    endif
  end subroutine export_block
!
  subroutine blocks(res,u_num,u_type,u_block,u_o,u_r,u_exist) !find blocks
    use protein_information
    !max_num_single  :  max # of possible blocks for one residue
    !res             :  residue information
    !u_num           :  num of blocks for 'res'
    !u_type          :  geometric type list 'res'
    !u_block         :  block type list for 'res'
    !u_o             :  origin for each block in 'res'
    !u_r             :  axes for each block in 'res'
    !u_exist         :  existance of each block in 'res' (may not be there owning to missing atom)
    type(residue_type),intent(in) :: res
    integer,intent(out) :: u_num
    integer,intent(out),dimension(max_num_single) :: u_type,u_block
    real(kind=dbl),intent(out),dimension(3,max_num_single) :: u_o
    real(kind=dbl),intent(out),dimension(3,3,max_num_single) :: u_r
    logical,intent(out),dimension(max_num_single) ::  u_exist

    real(kind=dbl),dimension(3) :: r1,r2,r3,r4,r5,r6
    real(kind=dbl),dimension(3) :: hatom,cb

    select case(res%res_type)
    case(1)                   !GLY:1,7,4
       u_num=3
       u_type=(/1,3,1,0,0/)
       u_block=(/1,7,4,0,0/)

       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)          
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,2)
       r2=res%r(:,2)
       r3=res%r(:,5)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
    case(2)                   !ALA:1,7,3
       u_num=3
       u_type=(/1,3,1,0,0/)
       u_block=(/1,7,3,0,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,2)
       r3=res%r(:,5)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
    case(3)                   !VAL:1,7,12
       u_num=3
       u_type=(/1,3,4,0,0/)
       u_block=(/1,7,12,0,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,6)
       r3=res%r(:,7)
       call u_ref(4,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
    case(4)                   !ILE: 1,7,12,3
       u_num=4
       u_type=(/1,3,4,1,0/)
       u_block=(/1,7,12,3,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,6)
       r3=res%r(:,7)
       call u_ref(4,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,8)
       r2=res%r(:,6)
       r3=res%r(:,8)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
    case(5)                   !LEU: 1,7,5,12
       u_num=4
       u_type=(/1,3,1,4,0/)
       u_block=(/1,7,5,12,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,2)
       r3=res%r(:,6)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,6)
       r2=res%r(:,7)
       r3=res%r(:,8)
       call u_ref(4,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
    case(6)                   !SER: 1,7,5,2
       u_num=4
       u_type=(/1,3,1,1,0/)
       u_block=(/1,7,5,2,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,2)
       r3=res%r(:,6)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,6)
       r2=res%r(:,5)
       r3=res%r(:,6)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
    case(7)                   !THR: 1,7,15
       u_num=3
       u_type=(/1,3,5,0,0/)
       u_block=(/1,7,15,0,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,6)
       r3=res%r(:,7)
       call u_ref(5,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
    case(8)                   !ASP: 1,7,5,11
       u_num=4
       u_type=(/1,3,1,4,0/)        
       u_block=(/1,7,5,11,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,2)
       r3=res%r(:,6)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,6)
       r2=res%r(:,7)
       r3=res%r(:,8)
       call u_ref(4,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
    case(9)                   !ASN: 1,7,5,14
       u_num=4
       u_type=(/1,3,1,5,0/)        
       u_block=(/1,7,5,14,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,2)
       r3=res%r(:,6)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,6)
       r2=res%r(:,7)
       r3=res%r(:,8)
       call u_ref(5,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
    case(10)                   !GLU: 1,7,6,11
       u_num=4
       u_type=(/1,3,2,4,0/)        
       u_block=(/1,7,6,11,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,6)
       call u_ref(2,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,7)
       r2=res%r(:,8)
       r3=res%r(:,9)
       call u_ref(4,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
    case(11)                   !GLN: 1,7,6,14
       u_num=4
       u_type=(/1,3,2,5,0/)
       u_block=(/1,7,6,14,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,6)
       call u_ref(2,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,7)
       r2=res%r(:,8)
       r3=res%r(:,9)
       call u_ref(5,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
    case(12)                  !LYS: 1,7,5,6,10
       u_num=5
       u_type=(/1,3,1,2,3/)        
       u_block=(/1,7,5,6,10/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,2)
       r3=res%r(:,6)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,6)
       r2=res%r(:,7)
       call u_ref(2,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
       r1=res%r(:,8)
       r2=res%r(:,9)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,5),u_r(:,:,5),u_exist(5))
    case(13)                  !ARG: 1,7,5,6,13
       u_num=5
       u_type=(/1,3,1,2,4/)
       u_block=(/1,7,5,6,13/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,2)
       r3=res%r(:,6)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,6)
       r2=res%r(:,7)
       call u_ref(2,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
       r1=res%r(:,8)
       r2=res%r(:,10)
       r3=res%r(:,11)
       call u_ref(4,r1,r2,r3,r4,r5,r6,u_o(:,5),u_r(:,:,5),u_exist(5))
    case(14)                   !CYS: 1,7,8
       u_num=3
       u_type=(/1,3,3,0,0/)
       u_block=(/1,7,8,0,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,6)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
    case(15)                  !MET: 1,7,6,9
       u_num=4
       u_type=(/1,3,2,3,0/)   
       u_block=(/1,7,6,9,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,6)
       call u_ref(2,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,7)
       r2=res%r(:,8)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
    case(16)                  !PHE: 1,7,5,16
       u_num=4
       u_type=(/1,3,1,6,0/)   
       u_block=(/1,7,5,16,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,2)
       r3=res%r(:,6)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,6)
       r2=res%r(:,7)
       r3=res%r(:,8)
       r4=res%r(:,9)
       r5=res%r(:,10)
       r6=res%r(:,11)
       call u_ref(6,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
    case(17)                  !TYR: 1,7,5,16,2
       u_num=5
       u_type=(/1,3,1,6,1/)  
       u_block=(/1,7,5,16,2/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,2)
       r3=res%r(:,6)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,6)
       r2=res%r(:,7)
       r3=res%r(:,8)
       r4=res%r(:,9)
       r5=res%r(:,10)
       r6=res%r(:,11)
       call u_ref(6,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
       r1=res%r(:,12)
       r2=res%r(:,11)
       r3=res%r(:,12)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,5),u_r(:,:,5),u_exist(5))
    case(18)                  !TRP: 1,7,5,18
       u_num=4
       u_type=(/1,3,1,8,0/)   
       u_block=(/1,7,5,18,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,2)
       r3=res%r(:,6)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,6)
       r2=res%r(:,9)
       r3=res%r(:,13)
       r4=res%r(:,14)
       call u_ref(8,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
    case(19)                  !HIS: 1,7,5,17
       u_num=4
       u_type=(/1,3,1,7,0/)   
       u_block=(/1,7,5,17,0/)
       r1=res%r(:,1)
       r2=res%r(:,1)
       r3=res%r(:,15)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
       r1=res%r(:,5)
       r2=res%r(:,2)
       r3=res%r(:,6)
       call u_ref(1,r1,r2,r3,r4,r5,r6,u_o(:,3),u_r(:,:,3),u_exist(3))
       r1=res%r(:,6)
       r2=res%r(:,9)
       r3=res%r(:,10)
       call u_ref(7,r1,r2,r3,r4,r5,r6,u_o(:,4),u_r(:,:,4),u_exist(4))
    case(20)                  !PRO: 7,19
       u_num=2
       u_type=(/3,9,0,0,0/)   
       u_block=(/7,19,0,0,0/)
       r1=res%r(:,3)
       r2=res%r(:,4)
       call u_ref(3,r1,r2,r3,r4,r5,r6,u_o(:,1),u_r(:,:,1),u_exist(1))
       r1=res%r(:,1)
       r2=res%r(:,2)
       r3=res%r(:,5)
       r4=res%r(:,6)
       r5=res%r(:,7)
       call u_ref(9,r1,r2,r3,r4,r5,r6,u_o(:,2),u_r(:,:,2),u_exist(2))
    case default
       u_num=0
    end select
  end subroutine blocks
!
  subroutine u_ref(u_type,r1,r2,r3,r4,r5,r6,u_o,u_r,u_exist)  
    !calculate reference frame
    use constant
    implicit none
    !u_type  :  geometric type of current block
    !r1~r6   :  atomic positions (maxinum 6 atoms required to cal. ref. frame)
    !u_o     :  origin of current block
    !u_r     :  axices of current block
    !u_exist    existance of current block
    integer,intent(in) ::u_type
    real(kind=dbl),intent(in),dimension(3) :: r1,r2,r3,r4,r5,r6
    real(kind=dbl),intent(out),dimension(3) :: u_o
    real(kind=dbl),intent(out),dimension(3,3) :: u_r    
    logical,intent(out) ::  u_exist

    integer i,j,k
    real(kind=dbl),dimension(3) :: u1,u2

    select case(u_type)
    case(1)                   !N(backbone),OH,-C,-CA-,-C-
       if((abs(r1(1)-large).lt.small).or.(abs(r2(1)-large).lt.small).or. &
            (abs(r3(1)-large).lt.small))then
          u_exist=.false.
       else
          u_o=r1
          u_r(:,1)=r3-r2
          u_r(:,1)=u_r(:,1)/sqrt(dot_product(u_r(:,1),u_r(:,1)))
          u_exist=.true.
       endif
    case(2,3)                 !-C-C- ; C=0(backbone),C-S,S-C,C-N
       if((abs(r1(1)-large).lt.small).or.(abs(r2(1)-large).lt.small))then
          u_exist=.false.
       else
          u_o=(r1+r2)/2
          u_r(:,1)=r2-r1
          u_r(:,1)=u_r(:,1)/sqrt(dot_product(u_r(:,1),u_r(:,1)))
          u_exist=.true.
       endif
    case(4,5,7)               !0-C-O,C-C-C,N-CN-N;0-C-N,O-C-C;CCNCN(5-ring)
       if((abs(r1(1)-large).lt.small).or.(abs(r2(1)-large).lt.small).or. &
            (abs(r3(1)-large).lt.small))then
          u_exist=.false.
       else
          u_o=(r1+r2+r3)/3
          u_r(:,1)=r3-r2
          u_r(:,1)=u_r(:,1)/sqrt(dot_product(u_r(:,1),u_r(:,1)))
          u_r(:,2)=r3-r1
          u_r(:,2)=u_r(:,2)-dot_product(u_r(:,1),u_r(:,2))*u_r(:,1)
          u_r(:,2)=u_r(:,2)/sqrt(dot_product(u_r(:,2),u_r(:,2)))
          call cross_product(u_r(:,1),u_r(:,2),u_r(:,3))
          u_exist=.true.
       endif
    case(6)                   !C6 (6-ring)
       if((abs(r1(1)-large).lt.small).or.(abs(r2(1)-large).lt.small).or. &
            (abs(r3(1)-large).lt.small).or.(abs(r4(1)-large).lt.small).or. &
            (abs(r5(1)-large).lt.small).or.(abs(r6(1)-large).lt.small))then
          u_exist=.false.
       else
          u_o=(r1+r6)/2
          u1=r6-r1
          u2=r2-r5
          call cross_product(u1,u2,u_r(:,1))
          !     can omit
!!          u_r(:,2)=u1-dot_product(u_r(:,1),u1)*u_r(:,1)
!!          u_r(:,2)=u_r(:,2)/sqrt(dot_product(u_r(:,2),u_r(:,2)))
!!          call cross_product(u_r(:,1),u_r(:,2),u_r(:,3))
          !     this part end
          u_exist=.true.
       endif
    case(8)                   !C8N (6-ring+5-ring)
       if((abs(r1(1)-large).lt.small).or.(abs(r2(1)-large).lt.small).or. &
            (abs(r3(1)-large).lt.small).or.(abs(r4(1)-large).lt.small))then
          u_exist=.false.
       else
          u_o=(r1+r2+r3+r4)/4.
          u_r(:,1)=r1+r3-r2-r4
          u_r(:,1)=u_r(:,1)/sqrt(dot_product(u_r(:,1),u_r(:,1)))
          u_r(:,2)=r3-r1
          u_r(:,2)=u_r(:,2)-dot_product(u_r(:,1),u_r(:,2))*u_r(:,1)
          u_r(:,2)=u_r(:,2)/sqrt(dot_product(u_r(:,2),u_r(:,2)))
          call cross_product(u_r(:,1),u_r(:,2),u_r(:,3))
          u_exist=.true.
       endif
    case(9)                   !NC4 (5-ring)
       if((abs(r1(1)-large).lt.small).or.(abs(r2(1)-large).lt.small).or. &
            (abs(r3(1)-large).lt.small).or.(abs(r4(1)-large).lt.small).or. &
            (abs(r5(1)-large).lt.small))then
          u_exist=.false.
       else
          u_o=(r1+r2+r3+r4+r5)/5
          u1=r2-r4
          u2=r3-r5
          call cross_product(u1,u2,u_r(:,3))
          u_r(:,1)=u1+u2
          u_r(:,1)=u_r(:,1)-dot_product(u_r(:,3),u_r(:,1))*u_r(:,3)
          u_r(:,1)=u_r(:,1)/sqrt(dot_product(u_r(:,1),u_r(:,1)))
          call cross_product(u_r(:,3),u_r(:,1),u_r(:,2))
          u_exist=.true.
       endif
    end select
    
  end subroutine u_ref
!
  subroutine u_atomlist(res_type,u_block,u_natom,u_atom) 
    !atom indeces for a block
    implicit none
    !res_type  :  residue type
    !u_block   :  block type
    !u_natom   :  # of atoms
    !u_atom    :  atom indeces
    integer,intent(in)::res_type,u_block
    integer,intent(out)::u_natom
    integer,intent(out),dimension(9)::u_atom    

    select case(res_type)
    case(1)                   !GLY
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(4)
          u_natom=1
          u_atom(1)=2
       case default
          u_natom=0
       end select

    case(2)                   !ALA
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(3)
          u_natom=1
          u_atom(1)=5
       case default
          u_natom=0
       end select

    case(3)                   !VAL
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(12)
          u_natom=3
          u_atom(1)=5
          u_atom(2)=6
          u_atom(3)=7
       case default
          u_natom=0
       end select

    case(4)                   !ILE
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(12)
          u_natom=3
          u_atom(1)=5
          u_atom(2)=6
          u_atom(3)=7
       case(3)
          u_natom=1
          u_atom(1)=8
       case default
          u_natom=0
       end select

    case(5)                   !LEU
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(5)
          u_natom=1
          u_atom(1)=5
       case(12)
          u_natom=3
          u_atom(1)=6
          u_atom(2)=7
          u_atom(3)=8
       case default
          u_natom=0
       end select

    case(6)                   !SER
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(5)
          u_natom=1
          u_atom(1)=5
       case(2)
          u_natom=1
          u_atom(1)=6
       case default
          u_natom=0
       end select

    case(7)                   !THR
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(15)
          u_natom=3
          u_atom(1)=5
          u_atom(2)=6
          u_atom(3)=7
       case default
          u_natom=0
       end select

    case(8)                   !ASP
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(5)
          u_natom=1
          u_atom(1)=5
       case(11)
          u_natom=3
          u_atom(1)=6
          u_atom(2)=7
          u_atom(3)=8
       case default
          u_natom=0
       end select

    case(9)                   !ASN
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(5)
          u_natom=1
          u_atom(1)=5
       case(14)
          u_natom=3
          u_atom(1)=6
          u_atom(2)=7
          u_atom(3)=8
       case default
          u_natom=0
       end select

    case(10)                   !GLU
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(6)
          u_natom=2
          u_atom(1)=5
          u_atom(2)=6
       case(11)
          u_natom=3
          u_atom(1)=7
          u_atom(2)=8
          u_atom(3)=9
       case default
          u_natom=0
       end select

    case(11)                   !GLN
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(6)
          u_natom=2
          u_atom(1)=5
          u_atom(2)=6
       case(14)
          u_natom=3
          u_atom(1)=7
          u_atom(2)=8
          u_atom(3)=9
       case default
          u_natom=0
       end select

    case(12)                  !LYS
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(5)
          u_natom=1
          u_atom(1)=5
       case(6)
          u_natom=2
          u_atom(1)=6
          u_atom(2)=7
       case(10)
          u_natom=2
          u_atom(1)=8
          u_atom(2)=9
       case default
          u_natom=0
       end select

    case(13)                  !ARG
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(5)
          u_natom=1
          u_atom(1)=5
       case(6)
          u_natom=2
          u_atom(1)=6
          u_atom(2)=7
       case(13)
          u_natom=4
          u_atom(1)=8
          u_atom(2)=9
          u_atom(3)=10
          u_atom(4)=11
       case default
          u_natom=0
       end select

    case(14)                   !CYS
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(8)
          u_natom=2
          u_atom(1)=5
          u_atom(2)=6
       case default
          u_natom=0
       end select

    case(15)                  !MET
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(6)
          u_natom=2
          u_atom(1)=5
          u_atom(2)=6
       case(9)
          u_natom=2
          u_atom(1)=7
          u_atom(2)=8
       case default
       end select

    case(16)                  !PHE
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(5)
          u_natom=1
          u_atom(1)=5
       case(16)
          u_natom=6
          u_atom(1)=6
          u_atom(2)=7
          u_atom(3)=8
          u_atom(4)=9
          u_atom(5)=10
          u_atom(6)=11
       case default
          u_natom=0
       end select

    case(17)                  !TYR
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(5)
          u_natom=1
          u_atom(1)=5
       case(16)
          u_natom=6
          u_atom(1)=6
          u_atom(2)=7
          u_atom(3)=8
          u_atom(4)=9
          u_atom(5)=10
          u_atom(6)=11
       case(2)
          u_natom=1
          u_atom(1)=12
       case default
          u_natom=0
       end select

    case(18)                  !TRP
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(5)
          u_natom=1
          u_atom(1)=5
       case(18)
          u_natom=9
          u_atom(1)=6
          u_atom(2)=7
          u_atom(3)=8
          u_atom(4)=9
          u_atom(5)=10
          u_atom(6)=11
          u_atom(7)=12
          u_atom(8)=13
          u_atom(9)=14
       case default
          u_natom=0
       end select

    case(19)                  !HIS
       select case(u_block)
       case(1)
          u_natom=1
          u_atom(1)=1
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(5)
          u_natom=1
          u_atom(1)=5
       case(17)
          u_natom=5
          u_atom(1)=6
          u_atom(2)=7
          u_atom(3)=8
          u_atom(4)=9
          u_atom(5)=10
       case default
          u_natom=0
       end select

    case(20)                  !PRO
       select case(u_block)
       case(7)
          u_natom=2
          u_atom(1)=3
          u_atom(2)=4
       case(19)
          u_natom=5
          u_atom(1)=1
          u_atom(2)=2
          u_atom(3)=5
          u_atom(4)=6
          u_atom(5)=7
       case default
          u_natom=0
       end select
    end select
  end subroutine u_atomlist
!
  subroutine u_classify(u_type1,u_block1,u_o1,u_r1 &  !classify block pairs
       ,u_type2,u_block2,u_o2,u_r2,u_class)
!,u_dirall)
    implicit none
    !u_type1,2   :   geometic type
    !u_block1,2  :   block type
    !u_o1,2      :   origins
    !u_r1,2      :   reference frames
    !u_class     :   classified type
    integer,intent(in) :: u_type1,u_block1,u_type2,u_block2
    real(kind=dbl),intent(in),dimension(3) :: u_o1,u_o2
    real(kind=dbl),intent(in),dimension(3,3) :: u_r1,u_r2
    integer,intent(out) :: u_class
!    integer,intent(out),dimension(3) :: u_dirall

    integer :: dom_dir1,dom_dir2,rotate,class1,class2

    !     classification of orientaion of two contacting blocks

    call dom_cube(u_type1,u_block1,u_o1,u_r1,u_o2,dom_dir1)
    call dom_cube(u_type2,u_block2,u_o2,u_r2,u_o1,dom_dir2)
!      print *,u_type1,u_type2,dom_dir1,dom_dir2
    call cube_rotation(u_block1,u_type1,u_o1,u_r1,dom_dir1, &
         u_block2,u_type2,u_o2,u_r2,dom_dir2,class1,class2,rotate)

!      u_class=u_map(rotate,class1,class2,u_block1,u_block2)

    call get_class(u_block1,u_type1,class1,u_block2,u_type2,class2 &
         ,rotate,u_class)
!    u_dirall(1)=dom_dir1
!    u_dirall(2)=dom_dir2
!    u_dirall(3)=rotate
  end subroutine u_classify
!
  subroutine u_classify_reduced(u_type1,u_block1,u_o1,u_r1 &  
       !classify block pairs
       ,u_type2,u_block2,u_o2,u_r2,u_class)
    implicit none
    !u_type1,2   :   geometic type
    !u_block1,2  :   block type
    !u_o1,2      :   origins
    !u_r1,2      :   reference frames
    !u_class     :   classified type
    integer,intent(in) :: u_type1,u_block1,u_type2,u_block2
    real(kind=dbl),intent(in),dimension(3) :: u_o1,u_o2
    real(kind=dbl),intent(in),dimension(3,3) :: u_r1,u_r2
    integer,intent(out) :: u_class

    integer :: dom_dir1,dom_dir2,rotate,class1,class2

    !     classification of orientaion of two contacting blocks

    call dom_cube(u_type1,u_block1,u_o1,u_r1,u_o2,dom_dir1)
    call dom_cube(u_type2,u_block2,u_o2,u_r2,u_o1,dom_dir2)
!      print *,u_type1,u_type2,dom_dir1,dom_dir2
    call cube_rotation(u_block1,u_type1,u_o1,u_r1,dom_dir1, &
         u_block2,u_type2,u_o2,u_r2,dom_dir2,class1,class2,rotate)
    u_class=u_map(rotate,class1,class2,u_block1,u_block2)

!    call get_class(u_block1,u_type1,class1,u_block2,u_type2,class2 &
!         ,rotate,u_class)
  end subroutine u_classify_reduced
!
  subroutine dom_cube(u_type,u_block,u_o,u_r,u_r2,dom_dir) 
    !derive dominant direction
    use common_sense
    use constant
    implicit none
    !u_type, u_block,u_o,u_r   : as above
    !u_r2                      : the origin of the partener
    !dom_dir                   : dominant direction
    integer,intent(in) :: u_type,u_block
    real(kind=dbl),intent(in),dimension(3) :: u_o,u_r2
    real(kind=dbl),intent(in),dimension(3,3) :: u_r
    integer,intent(out) :: dom_dir
    !rij : pair distance
    !vij : pair vector
    !cos_theta : cos(theta)
    !tan_fi,fi,rx,ry,bin_fi : related to fi
    !ind_xyz : indeces for calculating 'dom_dir'
    real(kind=dbl),parameter :: sin45 = 0.70710678D0, sin60 = 0.86602540D0,&
         sin30 = 0.500000D0
    real(kind=dbl) :: rij,cos_theta,tan_fi,fi,rx,ry
    real(kind=dbl),dimension(3) :: vij
    integer :: bin_fi
    integer,dimension(3) :: ind_xyz

    select case(u_type)
    case default
       dom_dir=1
!    case(2)
!       !     1:x+ ; 2: y ; 3: x-
!       vij=u_r2-u_o
!       rij=sqrt(dot_product(vij,vij))
!       cos_theta=dot_product(vij,u_r(:,1))/rij
!       if(cos_theta.ge.sin45)then
!          dom_dir=1
!       elseif(cos_theta.ge.-sin45)then            
!          dom_dir=2
!       else
!          dom_dir=3
!       endif
    case(1:3,6)                   !line
       !     1: x+ ; 2: x+y ; 3: y ; 4: x-y ; 5: x-
       vij=u_r2-u_o
       rij=sqrt(dot_product(vij,vij))
       cos_theta=dot_product(vij,u_r(:,1))/rij
       if(cos_theta.ge.sin60)then
          dom_dir=1
       elseif(cos_theta.ge.sin30)then
          dom_dir=2
       elseif(cos_theta.ge.-sin30)then
          dom_dir=3
       elseif(cos_theta.ge.-sin60)then
          dom_dir=4
       else
          dom_dir=5
       endif
    case(4,5,7:9)             !plane
       vij=u_r2-u_o
       rij=sqrt(dot_product(vij,vij))
       vij=vij/rij
       cos_theta=dot_product(vij,u_r(:,3))

       if(cos_theta.ge.sin60)then
          dom_dir=23
          return
       elseif(cos_theta.ge.sin30)then
          ind_xyz(3)=1
       elseif(cos_theta.ge.-sin30)then
          ind_xyz(3)=0
       elseif(cos_theta.ge.-sin60)then
          ind_xyz(3)=-1
       else
          dom_dir=5
          return
       endif

       rx=dot_product(vij,u_r(:,1))
       ry=dot_product(vij,u_r(:,2))

       if(abs(rx).lt.small)then
          if(ry.ge.zero)then
             fi=0.5
          else
             fi=1.5
          endif
       else
          tan_fi=ry/rx
          fi=atan(tan_fi)/pi
          if(rx.lt.0.)fi=fi+1.
          if(fi.lt.0.)fi=fi+2.
       endif

       bin_fi=int(fi*8)+1

       if(bin_fi.ge.17)then
          bin_fi=16
       endif

       select case(bin_fi)
       case(1,16)
          ind_xyz(1)=1
          ind_xyz(2)=0
       case(2,3)
          ind_xyz(1)=1
          ind_xyz(2)=1
       case(4,5)
          ind_xyz(1)=0
          ind_xyz(2)=1
       case(6,7)
          ind_xyz(1)=-1
          ind_xyz(2)=1
       case(8,9)
          ind_xyz(1)=-1
          ind_xyz(2)=0
       case(10,11)
          ind_xyz(1)=-1
          ind_xyz(2)=-1
       case(12,13)
          ind_xyz(1)=0
          ind_xyz(2)=-1
       case(14,15)
          ind_xyz(1)=1
          ind_xyz(2)=-1
       end select

       dom_dir=ind_xyz(3)*9+ind_xyz(2)*3+ind_xyz(1)+14
    end select
  end subroutine dom_cube
!
  subroutine cube_rotation(u_block1,u_type1,u_o1,u_r1,dom_dir1, &
       u_block2,u_type2,u_o2,u_r2,dom_dir2,class1,class2,rotate)
    !calculate rotation degree
    implicit none
    !variables as above
    integer,intent(in) :: u_block1,u_type1,u_block2,u_type2, &
         dom_dir1,dom_dir2
    real(kind=dbl),dimension(3) ::  u_o1,u_o2
    real(kind=dbl),dimension(3,3) :: u_r1,u_r2
    integer,intent(out) :: class1,class2,rotate
    !nplane1,2 : plane or not
    !sym1,2 : symmetry status
    !ui_1,2,uj_1,2 : reference axices
    logical :: nplane1,nplane2
    integer :: sym1,sym2
    real(kind=dbl),dimension(3) :: ui_1,ui_2,uj_1,uj_2
    !     rotate:
    !     dot involving case : 1
    !     line, line : 1, parallel; -1, antiparallel; 2, orthogonal
    !     line, plane : 1,-1,2,-2,3,-3: the dir that the head of line points to
    !     plane, plane :1,2,3,4; (-3:,symmetic direction)

    class1=u_class_sub(dom_dir1,u_type1)
    class2=u_class_sub(dom_dir2,u_type2)
    sym1=u_sym(dom_dir1,u_type1)
    sym2=u_sym(dom_dir2,u_type2)

    !      print *,dom_dir1,dom_dir2,class1,class2,sym1,sym2

    if((sym1.eq.1).or.(sym2.eq.1))then
       rotate=1
       return
    endif

    nplane1=u_plane(u_type1)
    nplane2=u_plane(u_type2)

    if(.not.nplane1)then
       if(.not.nplane2)then
          !     line,line
          select case(sym1+sym2)
          case(4,5)           !2,2   2,3
             call u_2_line(u_o1,u_r1(:,1),u_o2,u_r2(:,1),rotate)
          case(6)             !3,3
             call u_3_line(u_o1,u_r1(:,1),u_o2,u_r2(:,1),rotate)
          case default
             print *,'error line line'
             stop
          end select
       else
          !     line,plane
          ui_1=u_r1(:,1)
          call ch_ref(u_r2(:,1),u_r2(:,2),u_r2(:,3), &
               u_axis_1(:,dom_dir2,u_type2),uj_1)
          call ch_ref(u_r2(:,1),u_r2(:,2),u_r2(:,3), &
               u_axis_2(:,dom_dir2,u_type2),uj_2)

          select case(sym1)
          case(2)
             call u_2_plane(u_o1,ui_1,u_o2,uj_1,uj_2,rotate)
          case(3)
             select case(sym2)
             case(2)
                call u_2_plane(u_o1,ui_1,u_o2,uj_1,uj_2,rotate)
             case(-1,-2,-3)
                call u_3_plane_line(u_o1,ui_1,u_o2,uj_1,uj_2,rotate)
             case(4)
                call u_4_plane_line3(u_o1,ui_1,u_o2,uj_1,uj_2,rotate)
             case default
                print *,'error line plane 3'
                stop
             end select
          case default
             print *,'error line plane'
             stop
          end select
       endif
    else
       if(.not.nplane2)then
          !     plane,line
          uj_1=u_r2(:,1)
          call ch_ref(u_r1(:,1),u_r1(:,2),u_r1(:,3), &
               u_axis_1(:,dom_dir1,u_type1),ui_1)
          call ch_ref(u_r1(:,1),u_r1(:,2),u_r1(:,3), &
               u_axis_2(:,dom_dir1,u_type1),ui_2)
          uj_1=u_r2(:,1)

          select case(sym2)
          case(2)
             call u_2_plane(u_o2,uj_1,u_o1,ui_1,ui_2,rotate)
          case(3)
             select case(sym1)
             case(2)
                call u_2_plane(u_o2,uj_1,u_o1,ui_1,ui_2,rotate)
             case(-1,-2,-3)
                call u_3_plane_line(u_o2,uj_1,u_o1,ui_1,ui_2,rotate)
             case(4)
                call u_4_plane_line3(u_o2,uj_1,u_o1,ui_1,ui_2,rotate)
             case default
                print *,'error line plane 3'
                stop
             end select
          case default
             print *,'error line plane'
             stop
          end select
       else
          !     plane,plane
          call ch_ref(u_r1(:,1),u_r1(:,2),u_r1(:,3), &
               u_axis_1(:,dom_dir1,u_type1),ui_1)
          call ch_ref(u_r1(:,1),u_r1(:,2),u_r1(:,3), &
               u_axis_2(:,dom_dir1,u_type1),ui_2)
          call ch_ref(u_r2(:,1),u_r2(:,2),u_r2(:,3), &
               u_axis_1(:,dom_dir2,u_type2),uj_1)
          call ch_ref(u_r2(:,1),u_r2(:,2),u_r2(:,3), &
               u_axis_2(:,dom_dir2,u_type2),uj_2)

          select case(sym1)
          case(2)               
             select case(sym2)
             case(2,4)
                call u_2_plane(u_o1,ui_1,u_o2,uj_1,uj_2,rotate)
             case(-1,-2,-3)
                call u_2_plane(u_o2,uj_1,u_o1,ui_1,ui_2,rotate)
             case default
                print *,'error plane plane 2'
                stop
             end select
          case(-1,-2,-3)
             select case(sym2)
             case(2)
                call u_2_plane(u_o1,ui_1,u_o2,uj_1,uj_2,rotate)
             case(-1,-2,-3)
                call u_3_plane(u_o1,ui_1,ui_2,u_o2,uj_1,uj_2,rotate)
             case(4)
                call u_4_plane_sym3(u_o1,ui_1,ui_2, &
                     u_o2,uj_1,uj_2,rotate)
             case default
                print *,'error plane plane 3'
                stop
             end select
          case(4)
             select case(sym2)
             case(2)
                call u_2_plane(u_o2,uj_1,u_o1,ui_1,ui_2,rotate)
             case(-1,-2,-3)
                call u_4_plane_sym3(u_o2,uj_1,uj_2, &
                     u_o1,ui_1,ui_2,rotate)
             case(4)
                if(u_block1.eq.u_block2)then
                   if(class1.eq.class2)then
                      call u_4_plane_degen(u_o1,ui_1,ui_2,u_o2, &
                           uj_1,uj_2,rotate)
                   elseif(class1.lt.class2)then
                      call u_4_plane(u_o1,ui_1,ui_2,u_o2, &
                           uj_1,uj_2,rotate)
                   else
                      call u_4_plane(u_o2,uj_1,uj_2,u_o1, &
                           ui_1,ui_2,rotate)
                   endif
                else
                   call u_4_plane(u_o1,ui_1,ui_2,u_o2, &
                        uj_1,uj_2,rotate)
                endif
             case default
                print *,'error plane plane 4'
                stop
             end select
          case default
             print *,'error'
             stop
          end select
       endif

    endif
  end subroutine cube_rotation
!
  subroutine u_2_line(u_o1,u_x1,u_o2,u_x2,rotate)
    implicit none
    real(kind=dbl),intent(in),dimension(3):: u_x1,u_x2,u_o1,u_o2
    integer,intent(out) :: rotate
    real(kind=dbl),dimension(3) :: u_12,d12,d23
    real(kind=dbl) :: inner

    u_12=u_o2-u_o1
    call cross_product(u_x1,u_12,d12)
    call cross_product(u_x2,u_12,d23)
    inner=dot_product(d12,d23)
!      sign=dot_product(d12,u_x2)
    if(abs(inner).ge.0.7071)then
       rotate=1
    else
       rotate=2
    endif
  end subroutine u_2_line
!
  subroutine u_3_line(u_o1,u_x1,u_o2,u_x2,rotate)
    implicit none
    real(kind=dbl),intent(in),dimension(3):: u_x1,u_x2,u_o1,u_o2
    integer,intent(out) :: rotate
    real(kind=dbl),dimension(3) :: u_12,d12,d23
    real(kind=dbl) :: inner

    u_12=u_o2-u_o1
    call cross_product(u_x1,u_12,d12)
    call cross_product(u_x2,u_12,d23)
    inner=dot_product(d12,d23)
!     sign=dot_product(d12,u_x2)
    if(abs(inner).le.0.7071)then
       rotate=1
    elseif(inner.gt.0.7071)then
       rotate=2
    else
       rotate=3
    endif
  end subroutine u_3_line
!
  subroutine u_2_plane(u_o1,u_x1,u_o2,u_x2,u_y2,rotate)
    implicit none
    real(kind=dbl),intent(in),dimension(3):: u_x1,u_x2,u_y2,u_o1,u_o2
    integer,intent(out) :: rotate
    real(kind=dbl),dimension(3) :: u_12,d12,d23,d24
    real(kind=dbl) :: inner,inner2,dom_dir2

    u_12=u_o2-u_o1
    call cross_product(u_x1,u_12,d12)
    call cross_product(u_x2,u_12,d23)
    inner=dot_product(d12,d23)

    call cross_product(u_y2,u_12,d24)
    inner2=dot_product(d12,d24)
    if(abs(inner).ge.abs(inner2))then
       rotate=1
    else
       rotate=2
    endif
  end subroutine u_2_plane
!
  subroutine u_3_plane(u_o1,u_x1,u_y1,u_o2,u_x2,u_y2,rotate)
    implicit none
    real(kind=dbl),intent(in),dimension(3):: u_x1,u_y1,u_x2,u_y2,u_o1,u_o2
    integer,intent(out) :: rotate
    real(kind=dbl),dimension(3) :: u_12,d12,d23,d24,d25
    real(kind=dbl) :: inner,inner2,inner3,dom_dir2

    !     rotate = 1: two symmetic axes are orthogonal;
    !     rotate = 2: two non-symmetric axes are parrallel;
    !     rotate = 3: two non-symmetric axes are anti-parrallel

    u_12=u_o2-u_o1
    call cross_product(u_x1,u_12,d12)
    call cross_product(u_x2,u_12,d23)
    inner=dot_product(d12,d23)

    call cross_product(u_y2,u_12,d24)
    inner2=dot_product(d12,d24)

    if(abs(inner).le.abs(inner2))then
       rotate=1
    else
       call cross_product(u_y1,u_12,d25)
       inner3=dot_product(d25,d24)
       if(inner3.ge.0.)then
          rotate=2
       else
          rotate=3
       endif
    endif
  end subroutine u_3_plane
!
  subroutine u_3_plane_line(u_o1,u_x1,u_o2,u_x2,u_y2,rotate)
    implicit none
    real(kind=dbl),intent(in),dimension(3):: u_x1,u_x2,u_y2,u_o1,u_o2
    integer,intent(out) :: rotate
    real(kind=dbl),dimension(3) :: u_12,d12,d23,d24
    real(kind=dbl) :: inner,inner2,dom_dir2

    !     rotate = 1: two symmetic axes are orthogonal;
    !     rotate = 2: two non-symmetric axes are parrallel;
    !     rotate = 3: two non-symmetric axes are anti-parrallel

    u_12=u_o2-u_o1
    call cross_product(u_x1,u_12,d12)
    call cross_product(u_x2,u_12,d23)
    inner=dot_product(d12,d23)

    call cross_product(u_y2,u_12,d24)
    inner2=dot_product(d12,d24)

    if(abs(inner).ge.abs(inner2))then
       rotate=1
    else
       if(inner2.ge.0.)then
          rotate=2
       else
          rotate=3
       endif
    endif
  end subroutine u_3_plane_line
!
  subroutine u_4_plane(u_o1,u_x1,u_y1,u_o2,u_x2,u_y2,rotate)
    implicit none
    real(kind=dbl),intent(in),dimension(3):: u_x1,u_y1,u_x2,u_y2,u_o1,u_o2
    integer,intent(out) :: rotate
    real(kind=dbl),dimension(3) :: u_12,d12,d23,d24,d25
    real(kind=dbl) :: inner,inner2,inner3,inner4,dom_dir2
    !     rotate = 1: x1,x2  ,y1,y2
    !     rotate = 3: x1,-x2 ,y1,-y2
    !     rotate = 2: x1,y2  ,y1,-x2
    !     rotate = 4: x1,-y2 ,y1,x2

    !     rotate = 5: x1,-x2  ,y1,y2
    !     rotate = 7: x1,x2  ,y1,-y2
    !     rotate = 6: x1,y2 ,y1,x2
    !     rotate = 8: x1,-y2 ,y1,-x2

    u_12=u_o2-u_o1
    call cross_product(u_x1,u_12,d12)
    call cross_product(u_x2,u_12,d23)
    inner=dot_product(d12,d23) !x1,x2

    call cross_product(u_y2,u_12,d24)
    inner2=dot_product(d12,d24) !x1,y2

    call cross_product(u_y1,u_12,d25)

    if(abs(inner).ge.abs(inner2))then
       inner3=dot_product(d25,d24) !y1,y2
       if(inner.ge.0)then
          if(inner3.ge.0)then
             rotate=1
          else
             rotate=7
          endif
       else
          if(inner3.ge.0)then
             rotate=5
          else
             rotate=3
          endif
       endif
    else
       inner4=dot_product(d25,d23) !y1,x2
       if(inner2.ge.0)then
          if(inner4.ge.0)then
             rotate=6
          else
             rotate=2
          endif
       else
          if(inner4.ge.0)then
             rotate=4
          else
             rotate=8
          endif
       endif
    endif
  end subroutine u_4_plane
!
  subroutine u_4_plane_line3(u_o1,u_x1,u_o2,u_x2,u_y2,rotate)
    implicit none
    real(kind=dbl),intent(in),dimension(3):: u_x1,u_x2,u_y2,u_o1,u_o2
    integer,intent(out) :: rotate
    real(kind=dbl),dimension(3) :: u_12,d12,d23,d24
    real(kind=dbl) :: inner,inner2,dom_dir2

    !     rotate = 1: x1,x2
    !     rotate = 3: x1,-x2
    !     rotate = 2: x1,y2
    !     rotate = 4: x1,-y2

    u_12=u_o2-u_o1
    call cross_product(u_x1,u_12,d12)
    call cross_product(u_x2,u_12,d23)
    inner=dot_product(d12,d23)

    call cross_product(u_y2,u_12,d24)
    inner2=dot_product(d12,d24)

    if(abs(inner).ge.abs(inner2))then
       if(inner.ge.0)then
          rotate=1
       else
          rotate=3
       endif
    else
       if(inner2.ge.0.)then
          rotate=2
       else
          rotate=4
       endif
    endif
  end subroutine u_4_plane_line3
!
  subroutine u_4_plane_sym3(u_o1,u_x1,u_y1,u_o2,u_x2,u_y2,rotate)
    implicit none
    real(kind=dbl),intent(in),dimension(3):: u_x1,u_y1,u_x2,u_y2,u_o1,u_o2
    integer,intent(out) :: rotate
    real(kind=dbl),dimension(3) :: u_12,d12,d23,d24,d25
    real(kind=dbl) :: inner,inner2,inner3,inner4,dom_dir2

    !     rotate = 1: x1,+/-x2, y1,y2
    !     rotate = 3: x1,+/-x2, y1,-y2
    !     rotate = 2: x1,+/-y2, y1,-x2
    !     rotate = 4: x1,+/-y2, y1,x2

    u_12=u_o2-u_o1
    call cross_product(u_x1,u_12,d12)
    call cross_product(u_x2,u_12,d23)
    inner=dot_product(d12,d23)

    call cross_product(u_y2,u_12,d24)
    inner2=dot_product(d12,d24)

    call cross_product(u_y1,u_12,d25)

    if(abs(inner).ge.abs(inner2))then
       inner3=dot_product(d24,d25)
       if(inner3.ge.0)then
          rotate=1
       else
          rotate=3
       endif
    else
       inner4=dot_product(d23,d25)
       if(inner4.ge.0.)then
          rotate=4
       else
          rotate=2
       endif
    endif
  end subroutine u_4_plane_sym3
!
  subroutine u_4_plane_degen(u_o1,u_x1,u_y1,u_o2,u_x2,u_y2,rotate)
    implicit none
    real(kind=dbl),intent(in),dimension(3):: u_x1,u_y1,u_x2,u_y2,u_o1,u_o2
    integer,intent(out) :: rotate
    real(kind=dbl),dimension(3) :: u_12,d12,d23,d24,d25,d26
    real(kind=dbl) :: inner,inner2,inner3,inner4,dom_dir2
    !     rotate = 2: x1,x2  ,y1,y2
    !     rotate = 3: x1,-x2 ,y1,-y2
    !     rotate = 1: x1,y2  ,y1,-x2 ; or x1,-y2 ,y1,x2;

    !     rotate = 4: x1,-x2  ,y1,y2
    !     rotate = 6: x1,x2  ,y1,-y2
    !     rotate = 5: x1,y2 ,y1,x2
    !     rotate = 7: x1,-y2 ,y1,-x2

    u_12=u_o2-u_o1
    call cross_product(u_x1,u_12,d12)
    call cross_product(u_x2,u_12,d23)
    inner=dot_product(d12,d23)

    call cross_product(u_y2,u_12,d24)
    inner2=dot_product(d12,d24)

    call cross_product(u_y1,u_12,d25)

    if(abs(inner).ge.abs(inner2))then
       inner3=dot_product(d25,d24) !y1,y2
       if(inner.ge.0)then
          if(inner3.ge.0)then
             rotate=2
          else
             rotate=6
          endif
       else
          if(inner3.ge.0)then
             rotate=4
          else
             rotate=3
          endif
       endif
    else
       inner4=dot_product(d25,d23) !y1,x2
       if(inner2.ge.0)then
          if(inner4.ge.0)then
             rotate=5
          else
             rotate=1
          endif
       else
          if(inner4.ge.0)then
             rotate=1
          else
             rotate=7
          endif
       endif
    endif
  end subroutine u_4_plane_degen
!
  subroutine get_class(u_block1,u_type1,class1,u_block2, &  !derive class
       u_type2,class2,rotate,u_class)
    implicit none
    !variables as above
    integer,intent(in) :: u_block1,class1,u_type1, &
         u_block2,class2,u_type2,rotate
    integer,intent(out) :: u_class

    if(u_block1.eq.u_block2)then
       if(class1.le.class2)then
          u_class=u_clu_ind_sym(class1,class2,u_type1)+rotate
       else
          u_class=u_clu_ind_sym(class2,class1,u_type1)+rotate
       endif
    else
       u_class=u_clu_ind(class1,class2,u_type1,u_type2)+rotate
    endif
  end subroutine get_class
!
  subroutine gen_cluster_index  !generate general cluster indeces
    integer :: i,j,ii,jj,l,k,count,u1,u2
    integer,parameter,dimension(4,4) :: tot_rotate=reshape((/ &
         1,1,1,1, &
         1,2,2,2, &
         1,2,3,4, &
         1,2,4,8/),(/4,4/))
    open(unit=13,file='cluster_index.dat',status='replace')

    do i=1,num_type
       do j=1,num_type
          count=0
          do ii=1,u_class_num(i)
             do l=1,num_dom
                if(u_class_sub(l,i).ne.ii)cycle
                exit
             enddo
             do jj=1,u_class_num(j)
                do k=1,num_dom
                   if(u_class_sub(k,j).ne.jj)cycle
                   exit
                enddo
                if(u_sym(l,i).lt.0)then
                   u1=3
                else
                   u1=u_sym(l,i)
                endif
                if(u_sym(k,j).lt.0)then
                   u2=3
                else
                   u2=u_sym(k,j)
                endif

                u_clu_ind(ii,jj,i,j)=count
                count=count+tot_rotate(u1,u2)
                u_clu_orient(ii,jj,i,j)=tot_rotate(u1,u2)
             enddo
          enddo
          u_clu_num(i,j)=count
       enddo
    enddo

    do i=1,num_type
       count=0
       do ii=1,u_class_num(i)
          do l=1,num_dom
             if(u_class_sub(l,i).ne.ii)cycle
             exit
          enddo
          do jj=ii,u_class_num(i)
             do k=1,num_dom
                if(u_class_sub(k,i).ne.jj)cycle
                exit
             enddo
             if(u_sym(l,i).lt.0)then
                u1=3
             else
                u1=u_sym(l,i)
             endif
             if(u_sym(k,i).lt.0)then
                u2=3
             else
                u2=u_sym(k,i)
             endif

             u_clu_ind_sym(ii,jj,i)=count
             if((u1.eq.4).and.(u2.eq.4).and.(ii.eq.jj))then
                count=count+7
                u_clu_orient_sym(ii,jj,i)=7
             else
                count=count+tot_rotate(u1,u2)
                u_clu_orient_sym(ii,jj,i)=tot_rotate(u1,u2)
             endif

          enddo
       enddo
       u_clu_num_sym(i)=count
    enddo

    do i=1,num_type
       do j=1,num_type
          do k=1,u_class_num(i)
             do l=1,u_class_num(j)
                write(13,"(6I6)")i,j,k,l,u_clu_ind(k,l,i,j) &
                     ,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)
             write(13,"(5I6)")i,k,l,u_clu_ind_sym(k,l,i) &
                  ,u_clu_orient_sym(k,l,i)
          enddo
       enddo
    enddo

    write(13,"(A5)")"stat:"

    do i=1,num_type
       do j=1,num_type
          write(13,"(3I6)")i,j,u_clu_num(i,j)
       enddo
    enddo

    do i=1,num_type
       write(13,"(2I6)")i,u_clu_num_sym(i)
    enddo

    close(13)
  end subroutine gen_cluster_index
!
  subroutine export_unit_ref(u_block1,u_o1,u_r1,u_natom1,u_atom1, &
       indexi,u_block2,u_o2,u_r2,u_natom2,u_atom2,indexj)
    !output reference frame for two blocks
    use common_sense
    implicit none
    integer,intent(in) :: u_block1,u_block2,u_natom1,u_natom2,indexi,indexj
    integer,intent(in),dimension(:) :: u_atom1,u_atom2
    real(kind=dbl),intent(in),dimension(3) :: u_o1,u_o2
    real(kind=dbl),intent(in),dimension(3,3) :: u_r1,u_r2
    
    integer :: typei,typej
    integer :: count,i,ii,jj
    real(kind=dbl),dimension(3) :: r_ref,d_r,rij,u_o1_now,u_o2_now
    real(kind=dbl),dimension(3,3) :: u_r1_now,u_r2_now
    real(kind=dbl),dimension(3,max_atom) :: r1,r2
    
    typei=block_2_type(u_block1)
    typej=block_2_type(u_block2)

    do jj=1,u_natom2
       r2(:,jj)=u_pos(:,jj,u_block2)
    enddo

    rij=u_o1-u_o2

    call ch_ref_inv(u_r2,rij,u_o1_now)

    call ch_ref_inv(u_r2,u_r1(:,1),u_r1_now(:,1))
    call ch_ref_inv(u_r2,u_r1(:,2),u_r1_now(:,2))
    call ch_ref_inv(u_r2,u_r1(:,3),u_r1_now(:,3))

    do ii=1,u_natom1
       call ch_ref(u_r1_now(:,1),u_r1_now(:,2),u_r1_now(:,3) &
            ,u_pos(:,ii,u_block1),d_r)
       r1(:,ii)=u_o1_now+d_r
    enddo

    do i=1,u_natom2
       write(1,"(A6,I5,2X,A3,A1,A3,1X,I5,4X,3F8.3,A12)") &
            "ATOM   ",i,atom_list(indexj)(u_atom2(i)*3-2:u_atom2(i)*3) &
            ," ",res_list(indexj*4-3:indexj*4-1),1,r2(:,i),"  1.00 20.00"
    enddo

    count=0
    u_o2_now=0.
    u_r2_now=0.
    u_r2_now(1,1)=1.
    u_r2_now(2,2)=1.
    u_r2_now(3,3)=1.
    call export_block(count,1,typej,u_o2_now,u_r2_now)

    do i=1,u_natom1
       write(1,"(A6,I5,2X,A3,A1,A3,1X,I5,4X,3F8.3,A12)") &
            "ATOM   ",i,atom_list(indexi)(u_atom1(i)*3-2:u_atom1(i)*3) &
            ," ",res_list(indexi*4-3:indexi*4-1),2,r1(:,i),"  1.00 20.00"
    enddo

    call export_block(count,1,typei,u_o1_now,u_r1_now)

    write(1,"(A3)")"END"

  end subroutine export_unit_ref
!
  subroutine store_position(u_block,u_o,u_r,res,u_natom,u_atom,u_check)
    !used to generate position_pos.dat
    use protein_information
    implicit none
    type(residue_type),intent(in) :: res
    integer,intent(in) :: u_block,u_natom
    integer,intent(in),dimension(:) :: u_atom
    real(kind=dbl),intent(in),dimension(3) :: u_o
    real(kind=dbl),intent(in),dimension(3,3) :: u_r
    logical,intent(inout),dimension(num_block) :: u_check
    !nstop : stop collecting data
    !r_ref : position in reference fram
    logical :: nstop
    real(kind=dbl),dimension(3,3) :: u_r2
    real(kind=dbl),dimension(3) :: r_ref
    integer :: i

    nstop=.true.
    do i=1,19
       if(.not.u_check(i))then
          nstop=.false.
          exit
       endif
    enddo
    if(nstop)stop
    if(u_check(u_block))return
    u_check(u_block)=.true.

    open(unit=37,file='block_pos.dat',access='append')

    u_r2=u_r
    if(u_block.eq.16)then
       u_r2(:,2)=res%r(:,u_atom(1))-u_o
       u_r2(:,2)=u_r2(:,2)-dot_product(u_r2(:,1),u_r2(:,2))*u_r2(:,1)
       u_r2(:,2)=u_r2(:,2)/sqrt(dot_product(u_r2(:,2),u_r2(:,2)))
       call cross_product(u_r2(:,1),u_r2(:,2),u_r2(:,3))
    endif

    do i=1,u_natom
       call ref_pos(res%r(:,u_atom(i)),u_o,u_r2,r_ref)
       write(37,"(2I4,3F10.5)")u_block,i,r_ref
    enddo
    close(37)

  end subroutine store_position
!
  subroutine cross_product(a,b,c)  !perform cross product with normalization
    implicit none
    real(kind=dbl),intent(in),dimension(3) :: a,b
    real(kind=dbl),intent(out),dimension(3) :: c
    real(kind=dbl) :: norm

    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)
    norm=sqrt(c(1)**2+c(2)**2+c(3)**2)
    c=c/norm
    return
  end subroutine cross_product
!
  subroutine ch_ref(vx,vy,vz,r1,r2) !change reference
    !r1 in (vx,vy,vz)==> r2 in (nx,ny,nz)
    implicit none
    real(kind=dbl),intent(in),dimension(3) :: vx,vy,vz,r1
    real(kind=dbl),intent(out),dimension(3) :: r2
    
    r2=r1(1)*vx+r1(2)*vy+r1(3)*vz
  end subroutine ch_ref
!
  subroutine ch_ref_inv(r0,r1,r2) !change back reference
    !r0 in (nx,ny.nz) ==> r2 in r1 (3 components)
    implicit none
    real(kind=dbl),intent(in),dimension(3) :: r1
    real(kind=dbl),intent(in),dimension(3,3) :: r0
    real(kind=dbl),intent(out),dimension(3) :: r2
    
    r2(1)=dot_product(r0(:,1),r1)
    r2(2)=dot_product(r0(:,2),r1)
    r2(3)=dot_product(r0(:,3),r1)
  end subroutine ch_ref_inv
!
  subroutine ref_pos(r,origin,vr,r_ref) ! position in new reference frame
    ! r relative to reference frame of (origin,vr)
    implicit none
    real(kind=dbl),intent(in),dimension(3) :: r,origin
    real(kind=dbl),intent(out),dimension(3) :: r_ref

    real(kind=dbl),dimension(3) :: r1
    real(kind=dbl),dimension(3,3) :: vr
    
    r1=r-origin
    r_ref(1)=dot_product(r1,vr(:,1))
    r_ref(2)=dot_product(r1,vr(:,2))
    r_ref(3)=dot_product(r1,vr(:,3))
    
  end subroutine ref_pos
!
  character(len=255) function istr(number) !convert an integer to string
    implicit none
    integer,intent(in) :: number
    integer :: count,i,temp,dig,scale
    character(len=1),dimension(255) :: istring
    logical :: stopnow

    count=0
    istring=''
    stopnow=.false.
    scale=number
    do while (.not.stopnow)
       count=count+1
       temp=scale
       scale=scale/10
       dig=temp-scale*10
       istring(count)=char(dig+48)
       if(scale.eq.0)stopnow=.true.      
    enddo
    istr=istring(count)
    do i=1,count-1
       istr=istr(1:i)//istring(count-i)
    enddo
  end function istr

  subroutine add_h(res,res_pre)
    use protein_information
    use constant
    implicit none
    type(residue_type),intent(inout) :: res
    type(residue_type),intent(in) :: res_pre
    real(kind=dbl),dimension(3) :: vect

    if((abs(res_pre%r(1,3)-large).lt.small).or. &
         (abs(res_pre%r(1,4)-large).lt.small).or. &
         (abs(res%r(1,1)-large).lt.small))then
       res%r(1,15)=large
    else
       vect=res_pre%r(:,3)-res_pre%r(:,4)
       vect=vect/sqrt(dot_product(vect,vect))
       res%r(:,15)=res%r(:,1)+vect
    endif

  end subroutine add_h

  subroutine add_cb(res)
    use protein_information
    use common_sense
    use constant
    implicit none
    type(residue_type),intent(inout) :: res
    real(kind=dbl),dimension(3) :: va,vb,vect

    if((abs(res%r(1,1)-large).lt.small).or. &
         (abs(res%r(1,2)-large).lt.small).or. &
         (abs(res%r(1,3)-large).lt.small))then
       res%r(1,5)=large
    else
       va=res%r(:,3)-res%r(:,1)
       vb=res%r(:,2)-res%r(:,3)

       call inv_dihedral(va,vb,angle(1),angle(5),vect)
       res%r(:,5)=res%r(:,2)+vect*length(5)
    endif
  end subroutine add_cb

  subroutine inv_dihedral(r1,r2,bond,dihedral,r3)
    implicit none
    real(kind=dbl),dimension(3),intent(in) :: r1
    real(kind=dbl),intent(in) :: bond,dihedral
    real(kind=dbl),dimension(3),intent(out) :: r3
    real(kind=dbl),dimension(3),intent(inout) :: r2
    integer :: i
    real(kind=dbl),dimension(3) :: ry,rz,weight

    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

end module orientation

!program test
!  use orientation
!  use constant
!  implicit none

!  print *,large
!  call block_para
!  call test_para  
!  call gen_cluster_index

!end program test