cacmb/src/io.f90

module io

    use elements
    use parameters
    use atoms

    implicit none

    integer :: outfilenum = 0, infilenum = 0
    character(len=100) :: outfiles(10), infiles(10)

   public 
   contains

    subroutine get_out_file(filename)

        implicit none

        character(len=100), intent(in) :: filename
        character(len=100) :: temp_outfile
        character(len=1) :: overwrite
        logical :: file_exists

        !If no filename is provided then this function is called with none and prompts user input
        if (filename=='none') then 
            print *, "Please specify a filename or extension to output to:"
            read(*,*) temp_outfile
        else
            temp_outfile = filename
        end if

        !Infinite loop which only exists if user provides valid filetype
        overwrite = 'r'
        do while(.true.)

            !Check to see if file exists, if it does then ask user if they would like to overwrite the file
            inquire(file=trim(temp_outfile), exist=file_exists)
            if (file_exists) then 
                if (overwrite == 'r') print *, "File ", trim(temp_outfile), " already exists. Would you like to overwrite? (Y/N)"
                read(*,*) overwrite
                if((scan(overwrite, "n") > 0).or.(scan(overwrite, "N") > 0)) then 
                    print *, "Please specify a new filename with extension:"
                    read(*,*) temp_outfile
                    cycle
                else if((scan(overwrite, "y") > 0).or.(scan(overwrite, "Y") > 0)) then 
                    continue
                else
                    print *, "Please pick either y or n"
                    read(*,*) overwrite
                end if

            end if

            if (scan(temp_outfile,'.',.true.) == 0) then 
                print *, "No extension included on filename, please type a full filename that includes an extension."
                read(*,*) temp_outfile
                cycle
            end if
            select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:))
            case('xyz','lmp','vtk','cac')
                outfilenum=outfilenum+1
                outfiles(outfilenum) =  temp_outfile
                exit
            case default    
                print *, "File type: ", trim(temp_outfile(scan(temp_outfile,'.',.true.):)), " not currently accepted. ", &
                         "please input a filename with extension from following list: xyz, lmp, vtk, cac."
                read(*,*) temp_outfile

            end select
        end do

    end subroutine get_out_file 

    subroutine get_in_file(filename)

        implicit none

        character(len=100), intent(in) :: filename
        character(len=100) :: temp_infile
        logical :: file_exists

        !If no filename is provided then this function is called with none and prompts user input
        if (filename=='none') then 
            print *, "Please specify a filename with extension to read in:"
            read(*,*) temp_infile
        else
            temp_infile = filename
        end if

        !Infinite loop which only exists if user provides valid filetype
        do while(.true.)

            !Check to see if file exists, if it doesn't then ask the user for another input
            inquire(file=trim(temp_infile), exist=file_exists)
            if (.not.file_exists) then 
                print *, "The file ", temp_infile, "does not exist, please input an existing file to read in."
                read(*,*) temp_infile
                cycle
            end if

            if (scan(temp_outfile,'.',.true.) == 0) then 
                print *, "No extension included on filename, please type a full filename that includes an extension."
                read(*,*) temp_infile
                cycle
            end if
            select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:))
            case('cac')
                infilenum=infilenum+1
                infiles(infilenum) =  temp_infile
                exit
            case default    
                print *, "File type: ", trim(temp_infile(scan(temp_outfile,'.',.true.):)), " not currently accepted. ", &
                         "please input a filename with extension from following list: cac."
                read(*,*) temp_infile

            end select
        end do

    end subroutine get_in_file 

    subroutine write_out
        !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine

        integer :: i 

        !Find max esize which will be needed later
        call set_max_esize

        do i = 1, outfilenum
            !Pull out the extension of the file and call the correct write subroutine
            select case(trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))))
            case('xyz')
                call write_xyz(outfiles(i))
            case('lmp')
                call write_lmp(outfiles(i))
            case('vtk')
                call write_vtk(outfiles(i))
            case('cac')
                call write_lmpcac(outfiles(i))
            case default
                print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
                         " is not accepted for writing. Please select from: xyz and try again"
                stop 

            end select
        end do 
    end subroutine write_out



    subroutine write_xyz(file)
        !This is the simplest visualization subroutine, it writes out all nodal positions and atom positions to an xyz file
        character(len=100), intent(in) :: file 

        integer :: node_num, i, inod, ibasis

        open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')

        !Calculate total node number
        node_num=0
        do i = 1, ele_num
            node_num = node_num + basisnum(lat_ele(i))*ng_node(lat_ele(i))
        end do

        !Write total number of atoms + elements
        write(11, '(i16)') node_num+atom_num
         
        !Write comment line
        write(11, '(a)') "#Node + atom file created using cacmb"

        !Write nodal positions
        do i = 1, ele_num
            do inod = 1, ng_node(lat_ele(i))
                do ibasis = 1, basisnum(lat_ele(i))
                        write(11, '(a, 3f23.15)') basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i)
                end do 
            end do 
        end do 

        !Write atom positions
        do i = 1, atom_num
            write(11, '(a, 3f23.15)') type_atom(i), r_atom(:,i) 
        end do

        !Finish writing 
        close(11)
    end subroutine write_xyz

    subroutine write_lmp(file)
        !This subroutine writes out a .lmp style dump file
        character(len=100), intent(in) :: file
        integer :: write_num, i, iatom, type_interp(max_basisnum*max_esize**3)
        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), mass

        open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')

        !Comment line
        write(11, '(a)') '# lmp file made with cacmb'
        write(11, '(a)')
        !Calculate total atom number
        write_num = atom_num
        do i = 1,ele_num
            if(type_ele(i) == 'fcc') write_num = write_num + size_ele(i)**3
        end do
        !Write total number of atoms + elements
        write(11, '(i16, a)') write_num, ' atoms'
        !Write number of atom types 
        write(11, '(i16, a)') atom_types, ' atom types'

        write(11,'(a)') ' '
        !Write box bd
        write(11, '(2f23.15, a)') box_bd(1:2), ' xlo xhi'
        write(11, '(2f23.15, a)') box_bd(3:4), ' ylo yhi'
        write(11, '(2f23.15, a)') box_bd(5:6), ' zlo zhi'

        !Masses
        write(11, '(a)') 'Masses'

        write(11, '(a)') ' '
        do i =1, atom_types
            call atommass(type_to_name(i),mass)
            write(11, '(i16, f23.15)') i, mass
        end do
        write(11, '(a)') ' '

        !Write atom positions
        write(11, '(a)') 'Atoms'
        write(11, '(a)') ' '
        do i = 1, atom_num
            write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i) 
        end do

        !Write refined element atomic positions
        do i = 1, ele_num
            call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp)
            select case(trim(adjustl(type_ele(i))))
            case('fcc')
                do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
                    write(11, '(2i16, 3f23.15)') atom_num+iatom, type_interp(iatom), r_interp(:,iatom)
                end do
            end select
        end do
    end subroutine write_lmp

    subroutine write_lmpcac(file)
        !This subroutine writes out a .lmp style dump file
        character(len=100), intent(in) :: file
        integer :: write_num, i, inod, ibasis
        real(kind=dp) :: mass

1       format(i16, '   Eight_Node', 4i16)
2       format(i16, '   Atom', 4i16)
3       format(3i16,3f23.15)

        open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')

        !Comment line
        write(11, '(a)') '# CAC input file made with cacmb'
        write(11, '(a)')
        !Calculate total atom number
        write_num = atom_num + ele_num

        !Write total number of atoms + elements
        write(11, '(i16, a)') write_num, ' cac elements'
        !Write number of atom types 
        write(11, '(i16, a)') atom_types, ' atom types'

        write(11,'(a)') ' '
        !Write box bd
        write(11, '(2f23.15, a)') box_bd(1:2), ' xlo xhi'
        write(11, '(2f23.15, a)') box_bd(3:4), ' ylo yhi'
        write(11, '(2f23.15, a)') box_bd(5:6), ' zlo zhi'

        !Masses
        write(11, '(a)') 'Masses'

        write(11, '(a)') ' '
        do i =1, atom_types
            call atommass(type_to_name(i),mass)
            write(11, '(i16, f23.15, 2a)') i, mass, ' # ', type_to_name(i)
        end do
        write(11, '(a)') ' '

        write(11, '(a)') 'CAC Elements'
        write(11, '(a)') ' '

        !Write element nodal positions
        do i = 1, ele_num
            select case(trim(adjustl(type_ele(i))))
            case('fcc')
                !The first entry is the element specifier 
                write(11,1) i, basisnum(lat_ele(i)), size_ele(i), size_ele(i), size_ele(i)
                do ibasis = 1, basisnum(lat_ele(i))
                    do inod = 1, 8
                        !Nodal information for every node
                        write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i)
                    end do
                end do
            end select
        end do

        do i = 1, atom_num
            !Element specifier dictating that it is an atom
            write(11,2) ele_num+i, 1, 1, 1, 1
            !Write the atomic information
            write(11,3) 1, 1, type_atom(i), r_atom(:,i)
        end do
    end subroutine write_lmpcac

    subroutine write_vtk(file)
        !This subroutine writes out a vtk style dump file
        integer :: i, j, inod, ibasis
        character(len=100), intent(in) :: file

1       format('# vtk DataFile Version 4.0.1', / &
              'CAC output -- cg', / &
              'ASCII')
11      format('# vtk DataFile Version 4.0.1', / &
               'CACmb output -- atoms', / &
               'ASCII')
2       format('DATASET UNSTRUCTURED_GRID')
3       format('POINTS', i16, ' float')
4       format(/'CELLS', 2i16)
5       format(/'CELL_TYPES', i16)
12      format(/'CELL_DATA', i16)
16      format(/'POINT_DATA', i16)
17      format('SCALARS weight float', / &
               'LOOKUP_TABLE default')
18      format('SCALARS atom_type float', / &
               'LOOKUP_TABLE default')

20      format('SCALARS lattice_type float', /&
               'LOOKUP_TABLE default')

        !First we write the vtk file containing the atoms
        open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind')

        write(11, 11)
        write(11, 2)
        write(11, 3) atom_num
        do i = 1, atom_num
            write(11, '(3f23.15)') r_atom(:,i)
        end do 
        write(11,4) atom_num, atom_num*2
        do i = 1, atom_num
            write(11, '(2i16)') 1, i-1
        end do
        write(11, 5) atom_num
        do i = 1, atom_num
            write(11, '(i16)') 1 
        end do
        write(11, 16) atom_num
        write(11, 18)
        do i = 1, atom_num
            write(11, '(i16)') type_atom(i)
        end do
        close(11)

        open(unit=11, file='cg_'//trim(adjustl(file)), action='write', status='replace',position='rewind')
        write(11,1)
        write(11,2)
        write(11,3) node_num
        do i = 1, ele_num
            do inod=1, ng_node(lat_ele(i))  
                do ibasis = 1, basisnum(lat_ele(i))
                    write(11, '(3f23.15)') sum(r_node(:,:,inod,i),2)/basisnum(lat_ele(i))
                end do 
            end do 
        end do
        write(11, 4) ele_num, ele_num + node_num
        do i =1, ele_num
            write(11, '(9i16)') ng_node(lat_ele(i)), (j, j = (i-1)*ng_node(lat_ele(i)), i*ng_node(lat_ele(i))-1)
        end do
        write(11,5) ele_num
        do i = 1, ele_num
            if(trim(adjustl(type_ele(i))) == 'fcc') write(11, '(i16)') 12
        end do
        write(11,12) ele_num
        write(11,20)
        do i = 1, ele_num
            write(11, '(i16)') lat_ele(i)
        end do
        close(11)
    end subroutine


! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! READ SUBROUTINES !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     subroutine read_lmpcac(file, box_bd)
!         !This subroutine reads in a lmpcac file which can be used with different options and modes

!         !Arguments
!         character(len=100), intent(in) :: file
!         real(kind=wp), dimension(6), intent(out) :: box_bd

!         !Internal variables
!         character(len=1000) :: line
!         integer :: read_num, atom_lim, ele_lim

!         !Open the lmpcac file
!         open(unit=11, file=file, action='read', position='rewind')

!         !Skip header lines
!         read(11,*) line
!         read(11,*) line

!         !Read total number of elements
        

!     end subroutine read_lmpcac
end module io