CACmb/src/io.f90

module io

    use elements
    use parameters
    use atoms
    use box

    implicit none

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

   public 
   contains

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Subroutines for writing out data files !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    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
                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', 'mb', 'restart')
                outfilenum=outfilenum+1
                outfiles(outfilenum) =  temp_outfile
                exit
            case('cac')
                lmpcac = .true. 
                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 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('mb')
                call write_mb(outfiles(i))
            case('restart')
                call write_pycac(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 :: i, inod, ibasis

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

        !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

    subroutine write_pycac(file)
        !This subroutine writes restart files meant to be used with the McDowell Group CAC code.
        !NOTE: This code doesn't work for arbitrary number of basis atoms per node. It assumes that the 
        !each element has only 1 atom type at the node.
        character(len=100), intent(in) :: file
        integer :: interp_max, i, j,  lat_size, inod, ibasis, ip
        real(kind=dp) :: box_vec(3)

1  format('time' / i16, f23.15)
2  format('number of elements' / i16)
3  format('number of nodes' / i16)
4  format('element types' / i16)
5  format('number of atoms' / i16)
6  format('number of grains' / i16)
7  format('boundary ' / 3a1)
8  format('box bound' / 6f23.15)
9  format('box length' / 3f23.15)
10 format('box matrix')
11 format(3f23.15)
12 format('coarse-grained domain')
13 format('ie ele_type grain_ele lat_type_ele'/ 'ip ibasis type x y z')
14 format('atomistic domain' / 'ia grain_atom type_atom x y z')
15 format('maximum lattice periodicity length' / 3f23.15)
16 format('Number of lattice types and atom types '/ 2i16)
17 format('lattice type IDs')
18 format('lattice types for grains')
19 format('max nodes per element' / i16)
20 format('max interpo per element' / i16)
21 format('atom types to elements')

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

        write(11,1) timestep, total_time
        write(11,2) ele_num

        !Below writes the header information for the restart file

        !Calculate the max number of atoms per element
        select case(max_ng_node)
        case(8)
            interp_max = (max_esize)**3
        end select
        write(11,20) interp_max
        write(11,3) node_num
        write(11,19) max_ng_node
        write(11,4) lattice_types
        write(11,5) atom_num
        write(11,6) 1 !Grain_num is ignored 
        write(11,16) lattice_types, atom_types
        write(11,21)
        do i = 1, atom_types
            write(11,*) i, type_to_name(i)
        end do
        write(11,7) box_bc(1:1), box_bc(2:2), box_bc(3:3)
        write(11,18)
        write(11,'(2i16)') 1,1 !This is another throwaway line that is meaningless
        write(11,17)
        !This may have to be updated in the future but currently the only 8 node element is fcc
        do i = 1, lattice_types
            select case(ng_node(i))
            case(8)
                write(11, *) i, 'fcc'
            end select
        end do
        write(11,15) 1.0_dp, 1.0_dp, 1.0_dp !Another throwaway line that isn't needed
        write(11,8) box_bd
        write(11,9) box_bd(2)-box_bd(1), box_bd(4) - box_bd(3), box_bd(6)-box_bd(5)
        write(11,10)
        !Current boxes are limited to being rectangular 
        do i = 1,3 
            box_vec(:) = 0.0_dp
            box_vec(i) = box_bd(2*i) - box_bd(2*i-1)
            write(11,11) box_vec
        end do
        !We write this as box_mat ori and box_mat current
        do i = 1,3 
            box_vec(:) = 0.0_dp
            box_vec(i) = box_bd(2*i) - box_bd(2*i-1)
            write(11,11) box_vec
        end do

        !write the element information
        if(ele_num > 0) then
            write(11,12)
            do i = 1, lattice_types
                do j = 1, ele_num
                    if (lat_ele(j) == i) then 
                        lat_size = size_ele(j)-1
                        exit
                    end if
                end do 
                write(11,'(3i16)') i, lat_size, basis_type(1,i)
            end do
            ip = 0
            write(11,13)
            do i = 1, ele_num
                write(11, '(4i16)') i, lat_ele(i), 1, basis_type(1,lat_ele(i))
                do inod = 1, ng_node(lat_ele(i))
                    do ibasis = 1, basisnum(lat_ele(i))
                        ip = ip + 1
                        write(11, '(2i16, 3f23.15)') ip, ibasis, r_node(:, ibasis, inod, i)
                    end do
                end do
            end do 
        end if

        !Now write the atomic information
        if(atom_num /= 0) then 
            write(11,14)
            do i = 1, atom_num
                write(11, '(3i16, 3f23.15)') i, 1, type_atom(i), r_atom(:,i)
            end do
        end if

        close(11)
    end subroutine write_pycac

    subroutine write_mb(file)

        !This subroutine writes the cacmb formatted file which provides necessary information for building models
        character(len=100), intent(in) :: file

        integer :: i, j, k, inod, ibasis

        !Open the .mb file for writing
        open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')

        !First write the box boundary information
        !Write the global box boundaries
        write(11,*) box_bd(:)
        !Write the number of sub_boxes in the system
        write(11,*) sub_box_num
        !For every subbox write the orientation, sub box boundary, and sub_box_array_bds
        do i = 1, sub_box_num
            write(11,*) sub_box_ori(:,:,i)
            write(11,*) sub_box_bd(:,i)
            write(11,*) ((sub_box_array_bd(j,k,i), j = 1, 2), k = 1, 2)
        end do

        !Write the number of atom types in the current model and all of their names
        write(11,*) atom_types, (type_to_name(i)//' ', i=1, atom_types)
        !Write the number of lattice_types, basisnum and number of nodes for each lattice type
        write(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types)
        !Now for every lattice type write the basis atom types
        write(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = 1, lattice_types)

        !Now write the numbers of elements and atoms
        write(11,*) atom_num, ele_num

        !Write out atoms first
        do i = 1, atom_num
            write(11,*) i, type_atom(i), r_atom(:,i)
        end do

        !Write out the elements, this is written in two stages, one line for the element and then 1 line for 
        !every basis at every node
        do i = 1, ele_num
            write(11, *) i, lat_ele(i), size_ele(i), type_ele(i)
            do inod = 1, ng_node(lat_ele(i))
                do ibasis =1, basisnum(lat_ele(i))
                    write(11,*) inod, ibasis, r_node(:, ibasis, inod, i)
                end do
            end do
        end do

        close(11)
    end subroutine write_mb


    !!!!!!!!!!!!! Below are subroutines for reading files !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    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 or extension to output to:"
            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 does then ask user if they would like to overwrite the file
            inquire(file=trim(temp_infile), exist=file_exists)
            if (.not.file_exists) then 
                print *, "The file ", trim(adjustl(filename)), " does not exist. Please input a filename that exists"
                read(*,*) temp_infile
                cycle
            end if

            select case(temp_infile(scan(temp_infile,'.',.true.)+1:))
            case('xyz', 'lmp', 'vtk', 'mb')
                infilenum=infilenum+1
                infiles(infilenum) =  temp_infile
                exit
            case default    
                print *, "File type: ", trim(temp_infile(scan(temp_infile,'.',.true.):)), "not currently accepted. ", &
                         "please input a filename with extension from following list: mb."
                read(*,*) temp_infile

            end select
        end do

    end subroutine get_in_file 

    subroutine read_in(i, displace, temp_box_bd)
        !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine

        integer, intent(in) :: i 
        real(kind=dp), dimension(3), intent(in) :: displace
        real(kind=dp), dimension(6), intent(out) :: temp_box_bd

        !Pull out the extension of the file and call the correct write subroutine
        select case(trim(adjustl(infiles(i)(scan(infiles(i),'.',.true.)+1:))))
        case('mb')
            call read_mb(infiles(i), displace, temp_box_bd)
        case default
            print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
                     " is not accepted for writing. Please select from: mb and try again"
            stop 

        end select

    end subroutine read_in

    subroutine read_mb(file, displace, temp_box_bd)
        !This subroutine reads in an mb file for operation

        character(len=100), intent(in) :: file
        real(kind=dp), dimension(3), intent(in) :: displace
        real(kind = dp), dimension(6), intent(out) :: temp_box_bd

        integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles, new_atom_types, &
                   new_type_to_type(10), new_lattice_types
        character(len=100) :: etype
        real(kind=dp) ::  r(3), newdisplace(3)
        real(kind=dp), allocatable :: r_innode(:,:,:)
        character(len = 2) :: new_type_to_name(10)
        !First open the file
        open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')

        !Read in the box boundary and grow the current active box bd 
        read(11, *) temp_box_bd(:)

        do i = 1, 3
            newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
            temp_box_bd(2*i-1) = temp_box_bd(2*i-1) + newdisplace(i)
            temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i)
        end do 
        
        !Read in the number of sub_boxes and allocate the variables
        read(11, *) n

        if (sub_box_num == 0) then 
            call alloc_sub_box(n)
        else
            call grow_sub_box(n)
        end if

        !Read in subbox orientations and boundaries
        do i = 1,  n
            !Read in orientation with column major ordering
            read(11,*) ((sub_box_ori(j, k, sub_box_num+i), j = 1, 3), k = 1, 3)
            !Read in subbox boundaries
            read(11,*) sub_box_bd(:,sub_box_num+i)
            sub_box_bd(:,sub_box_num+i) = sub_box_bd(:, sub_box_num+i) + displace(:)
            !Read in sub_box_array_bd
            read(11,*) ((sub_box_array_bd(j, k, sub_box_num+i), j = 1, 2), k = 1, 2)
        end do

        !Add the existing element boundaries
        sub_box_array_bd(:,1,sub_box_num+1:) = sub_box_array_bd(:,1,sub_box_num+1:) + atom_num
        sub_box_array_bd(:,2,sub_box_num+1:) = sub_box_array_bd(:,2,sub_box_num+1:) + ele_num

        sub_box_num = sub_box_num + n

        !Read in the number of atom types and all their names
        read(11, *) new_atom_types, (new_type_to_name(i), i = 1, new_atom_types)
        !Now fit these into the global list of atom types, after this new_type_to_type is the actual global 
        !type of the atoms within this file
        do i = 1, new_atom_types
            call add_atom_type(new_type_to_name(i), new_type_to_type(i))
        end do
        !Read the number of lattice types, basisnum, and number of nodes for each lattice type
        read(11,*) new_lattice_types, (basisnum(i), i = lattice_types+1, lattice_types+new_lattice_types), &
                   (ng_node(i), i = lattice_types+1, lattice_types+new_lattice_types)
        !Define max_ng_node and max_basis_num
        max_basisnum = maxval(basisnum)
        max_ng_node = maxval(ng_node)
        !Read the basis atom types for every lattice
        read(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = lattice_types+1, lattice_types+new_lattice_types)
        !Convert the basis_atom types
        do j  = lattice_types+1, lattice_types+new_lattice_types
            do i = 1, basisnum(j)
                basis_type(i,j) = new_type_to_type(basis_type(i,j))
            end do 
        end do 
        !Read number of elements and atoms and allocate arrays
        read(11, *) in_atoms, in_eles
        call grow_ele_arrays(in_eles, in_atoms)
        allocate(r_innode(3,max_basisnum, max_ng_node))

        !Read the atoms 
        do i = 1, in_atoms
            read(11,*) j, type, r(:)
            call add_atom(new_type_to_type(type), r+newdisplace)
        end do

        !Read the elements
        do i = 1, in_eles
            read(11, *) n, type, size, etype
            do inod = 1, ng_node(type)
                do ibasis =1, basisnum(type)
                    read(11,*) j, k, r_innode(:, ibasis, inod)
                    r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace
                end do
            end do
            type = type + lattice_types
            call add_element(etype, size, type, r_innode)
        end do

        !Close the file being read
        close(11)

        !Only increment the lattice types if there are elements, if there are no elements then we 
        !just overwrite the arrays
        if(in_eles > 0) lattice_types = lattice_types + new_lattice_types
    end subroutine read_mb
end module io