Added writing to pycac restart file format
This commit is contained in:
Alex Selimov
parent
c291ec65b4
commit
b3a9577cc2
@ -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)
|
||||
|
||||
@ -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)
|
||||
|
||||
136
src/io.f90
136
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
|
||||
Loading…
x
Reference in New Issue
Block a user