program surfgen
! this program generates supercells of different plains of the 
! LJ bulk, for use in EPOCS
implicit none

integer, parameter :: dp=selected_real_kind(15,300)
integer :: i,j,k,num_cells,x_cell,y_cell,z_cell,gen_type,Ncell
integer, dimension(1:2) :: x_loc
real(kind=dp) :: sigma,length
real(kind=dp), dimension(1:4,1:3) :: prim_fcc
real(kind=dp), dimension(1:3,1:3) :: unit_cell
real(kind=dp), allocatable, dimension(:,:) :: primative
real(kind=dp), allocatable, dimension(:,:) :: output_cell
character(255) :: cells,file1,file2,gubbins
gen_type=0

! structure is FCC - space group Fm3-m

! just put in the positions of atoms:
!  (0,0,0)+; (0,0.5,0.5)+; (0.5,0,0.5)+; (0.5,0.5,0)+;

! specify the LJ parameter
print*, "please enter the Lennard-Jones parameter:"
read*, sigma
print*, "the lennard jones parameter is", sigma
length = sqrt(2.0_dp)*((2.0_dp)**(1.0_dp/6.0_dp))*sigma

! specify the number of cells
print*, "please enter the number of cells in AxBxC format"
read*, cells
j=1
do i=1,len(cells)
   if ((cells(i:i) == "x") .or. (cells(i:i) == "X")) then
      x_loc(j) = i
      j = j+1
   end if
   if (i == len(cells)) stop "delimiter not specified"
   if (j == 3) exit
end do

x_cell = 0
y_cell = 0
z_cell = 0

read (cells(1:x_loc(1)-1),'(i)') x_cell
read (cells(x_loc(1)+1:x_loc(2)-1),'(i)') y_cell
read (cells(x_loc(2)+1:len(cells)),'(i)') z_cell
!print*, trim(cells(1:x_loc(1)-1))
!print*, trim(cells(x_loc(1)+1:x_loc(2)-1))
!print*, trim(cells(x_loc(2)+1:len(cells)))

!print*, "you are outputting ", trim(cells(1:x_loc(1)-1)),"x", trim(cells(x_loc(1)+1:x_loc(2)-1)),"x",  trim(cells(x_loc(2)+1:len(cells)))," cells"
print '(a,3(i3,a))', "you are outputting ",x_cell," x",y_cell," x",z_cell," cells"

print*, "are you reading in from files(1), or using FCC(2)?"
read*, gen_type

if (gen_type==1) then
   !Ncell = 2
   print*, "please enter file name for cell:"
   read*, file1
   open(unit=70,file=trim(file1))
   do i=1,3
      read(70,*) unit_cell(i,:)
   end do
   close(70)
   print*, "please enter file name for primative:"
   read*, file2
   open(unit=71,file=trim(file2))
   read(71,*) Ncell
   allocate(primative(1:Ncell,1:3))
   read(71,*) gubbins
   do i=1,Ncell
      read(71,*) gubbins,primative(i,:)
   end do
   close(71)
else
print*, "generate FCC"
Ncell = 4
prim_fcc(1,:) = (/0.0_dp,0.0_dp,0.0_dp/)*length
prim_fcc(2,:) = (/0.0_dp,0.5_dp,0.5_dp/)*length
prim_fcc(3,:) = (/0.5_dp,0.0_dp,0.5_dp/)*length
prim_fcc(4,:) = (/0.5_dp,0.5_dp,0.0_dp/)*length

open(unit=99,file="primative.xyz")
write(99,*) 4
write(99,*) '*'
do i=1,4
   write(99,'(a,3f15.6)') 'Ne',prim_fcc(i,:)
end do
close(99)
allocate(primative(1:Ncell,1:3))
primative(:,:) = prim_fcc(:,:)
unit_cell(1,:) = (/length,0.0_dp,0.0_dp/)
unit_cell(2,:) = (/0.0_dp,length,0.0_dp/)
unit_cell(3,:) = (/0.0_dp,0.0_dp,length/)
end if



! generate output cell
allocate(output_cell(1:(Ncell*x_cell*y_cell*z_cell),1:3))
output_cell(:,:) = 0.0_dp
k=0
do i=1,Ncell
   k=k+1
   output_cell(i,:) = primative(i,:)
end do
do j=1,(x_cell-1)
   do i=1,Ncell
      k=k+1
      output_cell(Ncell*j+i,:) = output_cell(i,:) + j*unit_cell(1,:) 
   end do
end do
do j=1,(y_cell-1)
   do i=1,Ncell*x_cell
      k=k+1
      output_cell(Ncell*x_cell*j+i,:) = output_cell(i,:) + j*unit_cell(2,:)
   end do
end do
do j=1,(z_cell-1)
   do i=1,Ncell*x_cell*y_cell
      k=k+1
      output_cell(Ncell*x_cell*y_cell*j+i,:) = output_cell(i,:) + j*unit_cell(3,:)
   end do
end do


open(unit=98,file="outcell.xyz")
write(98,*) Ncell*x_cell*y_cell*z_cell
write(98,*) '*'
do i=1,Ncell*x_cell*y_cell*z_cell
   write(98,'(a,3f15.6)') 'Ne',output_cell(i,:)
end do
close(98)

open(unit=97,file="outcell.cell")
write(97,fmt='(3f9.4)') x_cell*unit_cell(1,:)
write(97,fmt='(3f9.4)') y_cell*unit_cell(2,:)
write(97,fmt='(3f9.4)') z_cell*unit_cell(3,:)
close(97)

print*, "the data can be found in outcell.cell and outcell.xyz"

print*, "the cell is:"
print '(f9.4,",",f9.4,","f9.4)', x_cell*unit_cell(1,:)
print '(f9.4,",",f9.4,","f9.4)', y_cell*unit_cell(2,:)
print '(f9.4,",",f9.4,","f9.4)', z_cell*unit_cell(3,:)

end program surfgen