program wc2emo
 implicit none
 character*256 :: title1, title2
 integer :: i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i
 integer :: glength
 integer :: nMO
 real :: a3, b3, c3
 real :: a4, b4, c4
 real :: a5, b5, c5
 real :: a6, b6, c6
 real :: a7, b7, c7, d7
 real :: a8, b8, c8, d8
 real :: a9, b9, c9, d9
 real :: a10, b10, c10, d10
 real :: a11, b11, c11, d11
 real :: a12, b12, c12, d12
 real, dimension(:,:), allocatable :: Wf
 real, dimension(:,:), allocatable :: TM
 real, dimension(6) :: p
 integer :: step
 character(len=16), dimension(4) :: Wfindex
 integer :: iWf
 integer :: EMOindex
 character(len=10) :: fileroot
 character(len=128) :: filename

 ! Define an array containing the result of the unitary transformation
 real, dimension(:), allocatable :: EMO     

 ! Indices of Wannier functions to read (4 and no more or less than four)
 ! These are read from Wfindex.dat
 open(unit=97, file='Wfindex.dat', status='old')
 do iWf = 1, 4
   read(97,*) Wfindex(iWf)
 end do
 close(97)

 ! Indices of EMOs to read
 EMOindex = 5

 ! Number of grid points (I assume that the grid is 60x60x60) 
 i4 = 60
 i5 = 60
 i6 = 60
 glength = i4 * i5 * i6

 ! Number of EMOs (= number  of KS orbitals)                                                             
 nMO = 8

 ! Wannier function array, Wf(4, gridpoints)
 allocate(Wf(4,glength))
 Wf = 0.0

 ! Read Wannier functions 
 fileroot = "WANNIER_1."
 do iWf = 1, 4
   filename = trim(fileroot) // trim(Wfindex(iWf)) // ".cube"
   print *, 'Reading ', filename
   open(unit=10, file=trim(filename), status='old')
   ! Read cube titles and lattice/grid sizes
   read(10,*) title1
   read(10,*) title2
   read(10,*) i3, a3, b3, c3
   read(10,*) i4, a4, b4, c4
   read(10,*) i5, a5, b5, c5
   read(10,*) i6, a6, b6, c6
   if (i4.ne.60 .and. i5.ne.60 .and. i6.ne.60) stop 'Change i4, i5, i6'
   read(10,*) i7, a7, b7, c7, d7
   read(10,*) i8, a8, b8, c8, d8
   read(10,*) i9, a9, b9, c9, d9
   read(10,*) i10, a10, b10, c10, d10
   read(10,*) i11, a11, b11, c11, d11
   read(10,*) i12, a12, b12, c12, d12
   read(10,*) Wf(iWf,:)
   close(10) 
 end do

 ! Transformation matrix array, TM(nMO,nMO)
 allocate (TM(nMO,nMO))
 TM=0.0

 ! Read Transformation matrix
 open(unit=98,file='TM.dat', status='old')
 read(98,*) TM
 close(98)
 write(99,*) TM
 close(99)
 print *,TM

 ! Allocate EMO array with the same size as gridpoints
 allocate (EMO(glength))

 ! Do unitary transformation, assuming transformation matrix is arranged in row
 ! if assuming transformation matrix is arranged in column, EMO(:)=EMO(:)+TM(EMOindex,iWf)*Wf(iWf,:) 
 EMO=0.0
 do iWf=1,4
 EMO(:)=EMO(:)+TM(iWf+4,EMOindex)*Wf(iWf,:) 
 print *, iWf+4, EMOindex, TM(iWf+4,EMOindex)
 end do

 ! Print out the EMO as a cube file
 open(unit=100, file='EMO5.cube', status='new')
 ! print out cube titles and lattice/grid sizes
 write(100,*) title1
 write(100,*) title2
 write(100,*) i3, a3, b3, c3
 write(100,*) i4, a4, b4, c4
 write(100,*) i5, a5, b5, c5
 write(100,*) i6, a6, b6, c6
 if (i4.ne.60 .and. i5.ne.60 .and. i6.ne.60) stop 'Change i4, i5, i6'
 write(100,*) i7, a7, b7, c7, d7
 write(100,*) i8, a8, b8, c8, d8
 write(100,*) i9, a9, b9, c9, d9
 write(100,*) i10, a10, b10, c10, d10
 write(100,*) i11, a11, b11, c11, d11
 write(100,*) i12, a12, b12, c12, d12
 ! print out EMO in lines each containing 6 numbers (before deallocating it)
 write(100,'(6e13.5)') EMO(:)
 close(100) 

 deallocate(EMO)

 deallocate(TM)

! Now Wf(i, :) contains the real-space coefficients of Wannier function i     

 deallocate(Wf)

end program wc2emo