Added read_pycac subroutine

src/io.f90

@@ -389,7 +389,7 @@ module io
 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')
+13 format('ie ele_type grain_ele lat_type_ele'/ 'ip ibasis 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)
@@ -521,11 +521,10 @@ module io
         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
+        !For every subbox write the orientation, sub box boundar
         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
@@ -590,13 +589,13 @@ module io
             end if
 
             select case(temp_infile(scan(temp_infile,'.',.true.)+1:))
-            case('xyz', 'lmp', 'vtk', 'mb')
+            case('restart', '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."
+                         "please input a filename with extension from following list: mb, restart."
                 read(*,*) temp_infile
 
             end select
@@ -615,6 +614,8 @@ module io
         select case(trim(adjustl(infiles(i)(scan(infiles(i),'.',.true.)+1:))))
         case('mb')
             call read_mb(infiles(i), displace, temp_box_bd)
+        case('restart')
+            call read_pycac(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"
@@ -643,7 +644,6 @@ module io
         !Read in the box boundary and grow the current active box bd 
         read(11, *) temp_box_bd(:)
 
-        print *, "displace", displace
         do i = 1, 3
             if (abs(displace(i)) > lim_zero) then 
                 newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
@@ -654,7 +654,6 @@ module io
             temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i)
         end do 
        
-        print *, "newdisplace", newdisplace 
         call grow_box(temp_box_bd)
         !Read in the number of sub_boxes and allocate the variables
         read(11, *) n
@@ -676,14 +675,8 @@ module io
                 sub_box_bd(2*j-1,sub_box_num+i) = sub_box_bd(2*j-1, sub_box_num+i) + displace(j)
                 sub_box_bd(2*j,sub_box_num+i) = sub_box_bd(2*j, sub_box_num+i) + displace(j)
             end do
-            !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
-
         !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 
@@ -715,8 +708,8 @@ module io
         new_loop:do i = 1, new_lattice_types
             old_loop:do j = 1, lattice_types
                 !First check all the lattice level variables
-                if ((basisnum(i) == basisnum(j)).and. &
+                if ((basisnum(lattice_types + i) == basisnum(j)).and. &
-                    (ng_node(i) == ng_node(j))) then 
+                    (ng_node(lattice_types + i) == ng_node(j))) then 
                     !Now check the basis  level variables
                     do ibasis =1, basisnum(i)
                         if(basis_type(ibasis,lattice_types+i) /= basis_type(ibasis,j)) then 
@@ -770,4 +763,181 @@ module io
     
         call set_max_esize
     end subroutine read_mb
+
+    subroutine read_pycac(file, displace, temp_box_bd)
+        !Add subroutine for reading in restart files from PyCAC. This code currently only 
+        !works for 8 node elements.
+        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, inod, ibasis, j, k, in_eles, in_atoms, ele_types, in_lat_num, in_atom_types, &
+                   atom_type_map(10), etype_map(10), etype, lat_type, new_lattice_map(10), &
+                   atom_type
+        real(kind=dp) :: newdisplace(3), r_in(3,1,8), r_in_atom(3), new_displace(3)
+        character(len=100) :: textholder, in_lattype_map(10) 
+        character(len=2) :: atomic_element
+        !First open the file
+        open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') 
+
+        !Disregard unneeded information
+        do i = 1, 3
+            read(11,*) textholder
+        end do 
+
+        !Read element number
+        read(11,*) in_eles 
+
+        !Discard info and read ng_max_node
+        do i = 1, 5 
+            read(11,*) textholder
+        end do 
+        read(11,*) max_ng_node
+
+        !Read element types (only needed inside this subroutine)
+        read(11,*) textholder
+        read(11,*) ele_types
+
+        !Read in atom num
+        read(11,*) textholder
+        read(11,*) in_atoms
+
+        !read in lattice_types and atom types
+        do i = 1,3 
+            read(11,*) textholder
+        end do
+        read(11,*) in_lat_num, in_atom_types
+
+
+        !Read define atom_types by name
+        read(11,*) textholder
+        do i = 1, in_atom_types
+            read(11,*) j, atomic_element
+            call add_atom_type(atomic_element, atom_type_map(i))
+        end do 
+
+        !Read in the boundary
+        read(11,*) textholder
+        read(11,*) box_bc
+
+        !Disregard useless info
+        do i = 1, 3
+            read(11,*) textholder
+        end do 
+
+        !Read in lat_type map
+        do i = 1, in_lat_num
+            read(11,*) j, in_lattype_map(i)
+            ng_node(lattice_types+i) = 8 !Only cubic elements are currently supported in pyCAC
+        end do 
+
+        !Disregard useless info
+        do i =1 , 3 
+            read(11,*) textholder
+        end do 
+
+        !Read box boundaries and displace them if necessary 
+        read(11,*) temp_box_bd(:)
+        do i = 1, 3
+            if (abs(displace(i)) > lim_zero) then 
+                newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
+            else
+                newdisplace(i)=displace(i)
+            end if 
+            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 
+
+        call grow_box(temp_box_bd)
+
+        !Allocate sub_box 
+        if (sub_box_num == 0) then 
+            call alloc_sub_box(1)
+        else
+            call grow_sub_box(1)
+        end if
+
+        !Because orientations and other needed sub_box information isn't really 
+        !present within the restart file we just default a lot of it.
+        sub_box_ori(:,:,sub_box_num+1) = identity_mat(3)
+        sub_box_bd(:, sub_box_num+1) = temp_box_bd
+
+        !Read in more useless info
+        do i = 1, 10
+            read(11,*) textholder
+        end do  
+
+        !Now read the ele_type to size and lat map
+        do i = 1, ele_types
+            read(11,*) j, etype_map(i)
+            etype_map(i) = etype_map(i) + 1
+        end do 
+
+
+        !Now set up the lattice types. In this code it assumes that lattice_type 1 maps to 
+        !atom type 1 because it only allows 1 atom per basis in pyCAC at the moment.
+        do i = 1, in_lat_num
+            basisnum(lattice_types+i) = 1
+            basis_type(1,lattice_types+i) = atom_type_map(i)
+        end do 
+
+        !Figure out the lattice type maps in case we have repeated lattice_types
+        k = lattice_types + 1
+        new_lattice_map(:) = 0
+        new_loop:do i = 1, in_lat_num
+            old_loop:do j = 1, lattice_types
+                !First check all the lattice level variables
+                if ((basisnum(lattice_types+i) == basisnum(j)).and. &
+                    (ng_node(lattice_types+i) == ng_node(j))) then 
+                    !Now check the basis  level variables
+                    do ibasis =1, basisnum(i)
+                        if(basis_type(ibasis,lattice_types+i) /= basis_type(ibasis,j)) then 
+                            cycle old_loop
+                        end if
+                    end do
+                    new_lattice_map(i) = j
+                    cycle new_loop
+                end if
+            end do old_loop
+            new_lattice_map(i) = k 
+            k = k+1
+        end do new_loop
+
+        !Read more useless data
+        read(11,*) textholder
+        read(11,*) textholder
+
+        !set max values and allocate variables 
+        max_basisnum = maxval(basisnum)
+        max_ng_node = maxval(ng_node)
+        call grow_ele_arrays(in_eles, in_atoms)
+        
+        !Now start reading the elements 
+        do i = 1, in_eles 
+            read(11,*) j, etype, k, lat_type 
+            do inod = 1, 8
+                read(11, *) j, k, r_in(:,1,inod)
+                r_in(:,1,inod) = r_in(:,1,inod) + newdisplace
+            end do 
+            call add_element(in_lattype_map(lat_type), etype_map(etype), new_lattice_map(lat_type), sub_box_num + 1, r_in)
+        end do
+
+        !Read useless data 
+        read(11,*)  textholder
+        read(11,*) textholder
+
+        do i = 1, in_atoms
+            read(11,*) j, k, atom_type, r_in_atom(:)
+            r_in_atom = r_in_atom + newdisplace
+            call add_atom(atom_type_map(atom_type), sub_box_num + 1, r_in_atom)
+        end do 
+        !Close file
+        close(11)
+
+        lattice_types = maxval(new_lattice_map)
+
+        sub_box_num = sub_box_num + 1
+
+        call set_max_esize
+    end subroutine read_pycac
 end module io