diff --git a/src/box.f90 b/src/box.f90 index ad5fe55..3043851 100644 --- a/src/box.f90 +++ b/src/box.f90 @@ -5,7 +5,7 @@ module box implicit none real(kind=dp) :: box_bd(6) !Global box boundaries - + character(len=3) :: box_bc !Box boundary conditions (periodic or shrinkwrapped) !The subbox variables contain values for each subbox, being the boxes read in through some !command. Currently only mode_merge will require sub_boxes, for mode_create it will always !allocate to only 1 sub_box @@ -14,12 +14,19 @@ module box real(kind=dp), allocatable :: sub_box_ori(:,:,:)!Orientations for each of the subboxes real(kind=dp), allocatable :: sub_box_bd(:,:)!Boundaries for each of the sub_boxes + + !Below are some simulation parameters which may be adjusted if reading in restart files + integer :: timestep=0 + real(kind=dp) :: total_time=0.0_dp + + public contains subroutine box_init !Initialize some box functions box_bd(:) = 0.0_dp + box_bc = 'ppp' end subroutine box_init subroutine alloc_sub_box(n) diff --git a/src/elements.f90 b/src/elements.f90 index 1195e58..22e6eff 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -11,8 +11,8 @@ module elements integer, allocatable :: size_ele(:), lat_ele(:) !Element siz real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array - integer :: ele_num=0 !Number of elements - integer :: node_num=0 !Total number of nodes + integer, save :: ele_num !Number of elements + integer, save :: node_num !Total number of nodes !Data structure used to represent atoms integer, allocatable :: type_atom(:)!atom type @@ -85,6 +85,9 @@ module elements max_basisnum = 0 basisnum(:) = 0 ng_node(:) = 0 + node_num = 0 + ele_num = 0 + atom_num = 0 end subroutine lattice_init subroutine cell_init(lapa,esize,ele_type, orient_mat, cell_mat) diff --git a/src/io.f90 b/src/io.f90 index e6c4645..6e914b2 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -58,7 +58,7 @@ module io cycle end if select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:)) - case('xyz', 'lmp', 'vtk', 'mb') + case('xyz', 'lmp', 'vtk', 'mb', 'restart') outfilenum=outfilenum+1 outfiles(outfilenum) = temp_outfile exit @@ -92,6 +92,8 @@ module io call write_vtk(outfiles(i)) case('mb') call write_mb(outfiles(i)) + case('restart') + call write_pycac(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" @@ -107,16 +109,10 @@ module io !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 + integer :: 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 @@ -273,6 +269,123 @@ module io 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,2) 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 @@ -321,8 +434,10 @@ module io end do end do + close(11) end subroutine write_mb + !!!!!!!!!!!!! Below are subroutines for reading files !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine get_in_file(filename) @@ -481,6 +596,9 @@ module io !Close the file being read close(11) - lattice_types = lattice_types + new_lattice_types + + !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 \ No newline at end of file