Added deform box option and allowed for reading time information from restart files

src/Makefile

@@ -2,7 +2,7 @@ FC=ifort
 FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays
 #FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays
 MODES=mode_create.o mode_merge.o mode_convert.o
-OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o 
+OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o opt_deform.o
 OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o  $(MODES) $(OPTIONS) call_option.o sorts.o
 
 .SUFFIXES:

src/box.f90

@@ -14,7 +14,6 @@ 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

src/call_option.f90

@@ -3,6 +3,7 @@ subroutine call_option(option, arg_pos)
     use opt_disl
     use opt_group
     use opt_orient
+    use opt_deform
     use opt_delete
     use box
     implicit none
@@ -31,6 +32,8 @@ subroutine call_option(option, arg_pos)
         bound_called = .true.
     case('-sbox_ori')
         call sbox_ori(arg_pos)
+    case('-deform')
+        call deform(arg_pos)
     case('-delete')
         call run_delete(arg_pos)
     case('-set_cac')

src/elements.f90

@@ -503,6 +503,7 @@ module elements
                 size_ele(sorted_index(i)) = 0
                 lat_ele(sorted_index(i)) = 0
                 sbox_ele(sorted_index(i)) = 0
+                tag_ele(sorted_index(i)) = 0
             else
                 node_num = node_num - ng_node(lat_ele(sorted_index(i)))
                 r_node(:,:,:,sorted_index(i)) = r_node(:,:,:,ele_num)
@@ -510,6 +511,7 @@ module elements
                 size_ele(sorted_index(i)) = size_ele(ele_num)
                 lat_ele(sorted_index(i)) = lat_ele(ele_num)
                 sbox_ele(sorted_index(i)) = sbox_ele(ele_num)
+                tag_ele(sorted_index(i)) = tag_ele(ele_num)
             end if            
             ele_num = ele_num - 1
         end do 

src/io.f90

@@ -491,7 +491,7 @@ module io
                         exit
                     endif
                 end do
-                write(11, '(4i16)') i, etype, 1, basis_type(1,lat_ele(i))
+                write(11, '(4i16)') tag_ele(i), etype, 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
@@ -505,7 +505,7 @@ module io
         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)
+                write(11, '(3i16, 3f23.15)') tag_atom(i), 1, type_atom(i), r_atom(:,i)
             end do
         end if
 
@@ -781,19 +781,19 @@ module io
 
         integer :: i, inod, ibasis, j, k, l, in_eles, in_atoms, ele_types, in_lat_num, in_atom_types, &
                    atom_type_map(100), etype_map(100), etype, lat_type, new_lattice_map(100), &
-                   atom_type
+                   atom_type, stat
         real(kind=dp) :: newdisplace(3), r_in(3,1,8), r_in_atom(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 the timestep information
+        read(11,*) textholder
+        read(11,*) timestep, total_time
 
         !Read element number
+        read(11,*) textholder
         read(11,*) in_eles 
 
         !Discard info and read ng_max_node
@@ -941,7 +941,11 @@ module io
             end if
 
             do i = 1, in_atoms
-                read(11,*) j, k, atom_type, r_in_atom(:)
+                read(11,*, iostat=stat) j, k, atom_type, r_in_atom(:)
+                if(stat > 0) then 
+                    print *, j
+                    stop
+                end if
                 r_in_atom = r_in_atom + newdisplace
                 call add_atom(j,atom_type_map(atom_type), sub_box_num + 1, r_in_atom)
             end do 

src/opt_deform.f90

@@ -0,0 +1,98 @@
+module opt_deform
+    !This module constains the deform option which applies a uniaxial strain to the system 
+    use parameters
+    use subroutines
+    use elements
+    use box
+
+    implicit none
+
+
+    real(kind=dp), save :: applied_strain
+    integer, save :: sdim
+
+    public 
+    contains
+
+    subroutine deform(arg_pos)
+        !This subroutine applies the simulation box deformation to the system
+        
+        integer, intent(inout) :: arg_pos
+
+        character(len=1) :: dims(3)
+        integer :: i, j, k
+        real(kind=dp) :: frac_atom(atom_num), frac_node(max_basisnum, max_ng_node, ele_num)
+
+        !initialize some variables
+        dims(1) = 'x'
+        dims(2) = 'y'
+        dims(3) = 'z'
+
+        !First parse the input command
+        call parse_deform(arg_pos)
+        print *, '-----------------------Option Deform------------------------'
+        !Now convert all positions in the specified dimension to fractional coordinates
+        do i = 1, atom_num
+            frac_atom(i) = (r_atom(sdim, i) - box_bd(2*sdim-1))/(box_bd(2*sdim)-box_bd(2*sdim-1))
+        end do
+        do i = 1, ele_num
+            do j = 1, ng_node(lat_ele(i))
+                do k = 1, basisnum(lat_ele(i))
+                    frac_node(k,j,i) = (r_node(sdim,k,j,i) - box_bd(2*sdim-1))/(box_bd(2*sdim)-box_bd(2*sdim-1))
+                end do
+            end do
+        end do
+
+        print *, "Original box bounds in ", dims(sdim), " are ", box_bd(2*sdim-1:2*sdim)
+        box_bd(2*sdim) = box_bd(2*sdim) + applied_strain
+        print *, "New box bounds are ", box_bd(2*sdim-1:2*sdim)
+
+        !Now reassign the positions
+        do i = 1, atom_num
+            r_atom(sdim,i) = frac_atom(i)*(box_bd(2*sdim)-box_bd(2*sdim-1)) + box_bd(2*sdim-1)
+        end do
+        do i = 1, ele_num
+            do j = 1, ng_node(lat_ele(i))
+                do k = 1, basisnum(lat_ele(i))
+                    r_node(sdim,k,j,i) = frac_node(k,j,i)*(box_bd(2*sdim)-box_bd(2*sdim-1)) + box_bd(2*sdim-1)
+                end do
+            end do
+        end do
+
+    end subroutine deform
+
+    subroutine parse_deform(arg_pos)
+
+        integer, intent(inout) :: arg_pos
+        integer :: arg_len
+        character(len=100) :: textholder
+
+        !Pull out the dimension to be strained 
+        arg_pos=arg_pos+1
+        call get_command_argument(arg_pos, textholder, arg_len)
+        if (arg_len == 0) stop "Missing dim in deform command"
+
+        select case(trim(adjustl(textholder)))
+        case('x','X') 
+            sdim = 1
+        case('y','Y')
+            sdim = 2
+        case('z','Z')
+            sdim = 3
+        case default
+            print *, "Dimension ", trim(adjustl(textholder)), " is not accepted. Please select either x, y, or z"
+        end select 
+
+        !Now pull out the strain vector, currently only accepts a real number by which to reduce the simulation cell size by in
+        !that dimension
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, textholder, arg_len)
+        if (arg_len == 0) stop "Missing strain in deform command"
+        read(textholder, *) applied_strain
+
+        arg_pos = arg_pos + 1
+
+    end subroutine parse_deform
+
+
+end module opt_deform