Merge pull request #5 from aselimov/ft--pycac-restart

Added writing to pycac restart file format

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)

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, save :: ele_num !Number of elements
-    integer :: node_num=0 !Total number of nodes
+    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)

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