Merge branch 'ft--read_pycac' into 'development'

Ft  read pycac

See merge request aselimov/cacmb!1

README.md

@@ -277,6 +277,14 @@ Delete requires the usage of additional keywords to specify which delete action
 
 This command will delete all overlapping atoms within a specific cutoff radius `rc_off`. This currently does not affect elements. 
 
+### Option sbox_ori
+
+```
+-sbox_ori sbox [hkl] [hkl] [hkl]
+```
+
+This option is primarily used when reading data from non .mb formats. This code simply sets the orientation variable for the specified sub box `sbox`.
+
 ****
 
 ## Position Specification

src/box.f90

@@ -11,7 +11,6 @@ module box
     !command. Currently only mode_merge will require sub_boxes, for mode_create it will always 
     !allocate to only 1 sub_box
     integer :: sub_box_num = 0
-    integer, allocatable :: sub_box_array_bd(:,:,:)!Boundaries in the atom and element arrays for each sub_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
 
@@ -38,11 +37,10 @@ module box
 
         integer :: i
 
-        allocate(sub_box_ori(3,3,n), sub_box_bd(6,n), sub_box_array_bd(2,2,n))
+        allocate(sub_box_ori(3,3,n), sub_box_bd(6,n))
         do i = 1, n
             sub_box_ori(:,:,i) = identity_mat(3)
             sub_box_bd(:,i) = 0.0_dp
-            sub_box_array_bd(:,:,i) = 1
         end do
     end subroutine alloc_sub_box
 
@@ -65,10 +63,6 @@ module box
         temp_bd(:, sub_box_num+1:) = 0.0_dp
         call move_alloc(temp_bd, sub_box_bd)
 
-        temp_array_bd(:,:,1:sub_box_num) = sub_box_array_bd
-        temp_array_bd(:,:,sub_box_num+1:) = 1
-        call move_alloc(temp_array_bd, sub_box_array_bd)
-
         return
     end subroutine grow_sub_box
 

src/call_option.f90

@@ -29,6 +29,8 @@ subroutine call_option(option, arg_pos)
         call get_command_argument(arg_pos, box_bc)
         arg_pos=arg_pos+1
         bound_called = .true.
+    case('-sbox_ori')
+        call sbox_ori(arg_pos)
     case('-delete')
         call run_delete(arg_pos)
     case default

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

src/mode_create.f90

@@ -158,9 +158,6 @@ module mode_create
         sub_box_num = 1
         sub_box_ori(:,:,1) = orient
         sub_box_bd(:,1) = box_bd
-        sub_box_array_bd(1,:,1) = 1
-        sub_box_array_bd(2,1,1) = atom_num
-        sub_box_array_bd(2,2,1) = ele_num
     end subroutine create
     !This subroutine parses the command and pulls out information needed for mode_create
     subroutine parse_command(arg_pos)

src/opt_group.f90

@@ -364,7 +364,7 @@ module opt_group
 
 
         if (max_remesh) then 
-            max_loops = (remesh_size-2)/2
+            max_loops = (remesh_size-3)/2
         else 
             max_loops = 1
         end if

src/opt_orient.f90

@@ -133,4 +133,29 @@ module opt_orient
         box_bd = orig_box_bd
     end subroutine unorient
 
+    subroutine sbox_ori(arg_pos)
+        integer, intent(inout) :: arg_pos
+
+        integer :: i, sbox_in, arg_len
+        real(kind = dp) :: new_orient(3,3)
+        character(len=100) :: textholder 
+        character(len=8) :: ori_string
+
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, textholder,arg_len)
+        if (arg_len== 0) stop 'Missing sbox in sbox_ori command'
+        read(textholder,*) sbox_in
+
+        do i = 1, 3
+            arg_pos = arg_pos + 1
+            call get_command_argument(arg_pos, ori_string, arg_len)
+            if (arg_len == 0) print *, "Missing orientation vector in -orient option"
+            call parse_ori_vec(ori_string, new_orient(i,:))
+            new_orient(i,:) = new_orient(i,:) / norm2(new_orient(i,:))
+        end do
+
+        sub_box_ori(:,:,sbox_in) = new_orient
+        arg_pos = arg_pos + 1
+
+    end subroutine sbox_ori
 end module opt_orient