Merge pull request #22 from aselimov/ft--option-orient

Ft  option orient

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
 #FFLAGS=-mcmodel=large -Ofast -no-wrap-margin
 MODES=mode_create.o mode_merge.o mode_convert.o
-OPTIONS=opt_disl.o opt_group.o
+OPTIONS=opt_disl.o opt_group.o opt_orient.o
 OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o  $(MODES) $(OPTIONS) call_option.o
 
 .SUFFIXES:
@@ -40,4 +40,4 @@ call_mode.o : $(MODES)
 call_option.o : $(OPTIONS)
 $(MODES) $(OPTIONS) subroutines.o io.o : atoms.o box.o
 $(MODES) main.o : io.o
-testfuncs.o elements.o mode_create.o $(OPTIONS): subroutines.o
+testfuncs.o elements.o mode_create.o $(OPTIONS) $(MODES): subroutines.o

src/call_option.f90

@@ -2,6 +2,7 @@ subroutine call_option(option, arg_pos)
     use parameters
     use opt_disl
     use opt_group
+    use opt_orient
     use box
     implicit none
 
@@ -15,14 +16,16 @@ subroutine call_option(option, arg_pos)
         call group(arg_pos)
     case('-ow')
         arg_pos = arg_pos + 1
-        continue
     case('-wrap')
         arg_pos = arg_pos + 1
-        continue
+    case('-orient')
+        call orient(arg_pos)
+    case('-unorient')
+        call unorient
+        arg_pos = arg_pos + 1
     case('-boundary')
         arg_pos=arg_pos+1
         call get_command_argument(arg_pos, box_bc)
-        print *, box_bc
         arg_pos=arg_pos+1
     case default
         print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.'

src/elements.f90

@@ -5,6 +5,7 @@ module elements
     use parameters
     use functions
     use subroutines
+    use box
     implicit none
 
     !Data structures used to represent the CAC elements. Each index represents an element
@@ -17,6 +18,7 @@ module elements
 
     !Data structure used to represent atoms 
     integer, allocatable ::  type_atom(:)!atom type 
+    integer, allocatable :: sbox_atom(:)
     real(kind =dp),allocatable :: r_atom(:,:) !atom position
     integer :: atom_num=0 !Number of atoms
 
@@ -153,7 +155,7 @@ module elements
 
         if(m > 0) then 
             !Allocate atom arrays
-            allocate(type_atom(m), r_atom(3,m), stat=allostat)
+            allocate(type_atom(m), sbox_atom(m), r_atom(3,m), stat=allostat)
             if(allostat > 0) then 
                 print *, "Error allocating atom arrays in elements.f90 because of: ", allostat
                 stop
@@ -211,6 +213,11 @@ module elements
             temp_int(atom_size+1:) = 0
             call move_alloc(temp_int, type_atom)
 
+            allocate(temp_int(m+atom_num+buffer_size))
+            temp_int(1:atom_size) = sbox_atom
+            temp_int(atom_size+1:) = 0
+            call move_alloc(temp_int, sbox_atom)
+
             allocate(temp_real(3,m+atom_num+buffer_size))
             temp_real(:,1:atom_size) = r_atom
             temp_real(:, atom_size+1:) = 0.0_dp
@@ -237,9 +244,9 @@ module elements
 
     end subroutine add_element
 
-    subroutine add_atom(type, r)
+    subroutine add_atom(type, sbox, r)
         !Subroutine which adds an atom to the atom arrays
-        integer, intent(in) :: type
+        integer, intent(in) :: type, sbox
         real(kind=dp), intent(in), dimension(3) :: r
 
         atom_num = atom_num+1
@@ -247,6 +254,7 @@ module elements
         call grow_ele_arrays(0,1)
         type_atom(atom_num) = type
         r_atom(:,atom_num) = r(:)
+        sbox_atom(atom_num) = sbox
 
     end subroutine add_atom
 
@@ -449,4 +457,38 @@ module elements
         end do 
     end subroutine wrap_atoms
 
+    subroutine def_new_box
+        !This subroutine calculates new box boundaries based on minimum and maximum nodal/atomic positions
+        integer :: i, j, inod, ibasis
+
+        real(kind=dp) :: max_bd(3), min_bd(3)
+
+        max_bd(:) = -huge(1.0_dp)
+        min_bd(:) = huge(1.0_dp)
+
+        do i = 1, atom_num
+            do j = 1, 3
+                if (r_atom(j,i) > max_bd(j)) max_bd(j) = r_atom(j,i) + lim_zero
+                if (r_atom(j,i) < min_bd(j)) min_bd(j) = r_atom(j,i) - lim_zero
+            end do 
+        end do 
+
+        do i = 1, ele_num
+            do inod = 1, ng_node(lat_ele(i))
+                do ibasis = 1, basisnum(lat_ele(i))
+                    do j = 1, 3
+                        if (r_node(j,ibasis,inod,i) > max_bd(j)) max_bd(j) = r_node(j,ibasis,inod,i) + lim_zero
+                        if (r_node(j,ibasis,inod,i) < min_bd(j)) min_bd(j) = r_node(j,ibasis,inod,i) -lim_zero
+
+                    end do 
+                end do 
+            end do 
+        end do 
+
+        do j = 1, 3
+            box_bd(2*j) = max_bd(j)
+            box_bd(2*j-1) = min_bd(j)
+        end do 
+    end subroutine
+
 end module elements

src/io.f90

@@ -8,7 +8,7 @@ module io
     implicit none
 
     integer :: outfilenum = 0, infilenum = 0
-    character(len=100) :: outfiles(10), infiles(10)
+    character(len=100) :: outfiles(100), infiles(100)
     logical :: force_overwrite
 
    public 
@@ -130,14 +130,14 @@ module io
         do i = 1, ele_num
             do inod = 1, ng_node(lat_ele(i))
                 do ibasis = 1, basisnum(lat_ele(i))
-                        write(11, '(a, 3f23.15)') basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i)
+                        write(11, '(i16, 3f23.15)') basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i)
                 end do 
             end do 
         end do 
 
         !Write atom positions
         do i = 1, atom_num
-            write(11, '(a, 3f23.15)') type_atom(i), r_atom(:,i) 
+            write(11, '(i16, 3f23.15)') type_atom(i), r_atom(:,i) 
         end do
 
         !Finish writing 
@@ -520,7 +520,7 @@ module io
 
         !Write out atoms first
         do i = 1, atom_num
-            write(11,*) i, type_atom(i), r_atom(:,i)
+            write(11,*) i, type_atom(i), sbox_atom(i), r_atom(:,i)
         end do
 
         !Write out the elements, this is written in two stages, one line for the element and then 1 line for 
@@ -713,8 +713,8 @@ module io
 
         !Read the atoms 
         do i = 1, in_atoms
-            read(11,*) j, type, r(:)
-            call add_atom(new_type_to_type(type), r+newdisplace)
+            read(11,*) j, type, sbox, r(:)
+            call add_atom(new_type_to_type(type), sbox, r+newdisplace )
         end do
 
         !Read the elements

src/mode_create.f90

@@ -135,7 +135,7 @@ module mode_create
             if(lat_atom_num > 0) then 
                 do i = 1, lat_atom_num
                     do ibasis = 1, basisnum(1)
-                        call add_atom(basis_type(ibasis, 1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis))
+                        call add_atom(basis_type(ibasis, 1), 1, (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis))
                     end do
                 end do
                 deallocate(r_atom_lat)
@@ -205,21 +205,23 @@ module mode_create
                 do i = 1, 3
                     call get_command_argument(arg_pos, orient_string, arglen)
                     if (arglen==0) STOP "Missing orientation in orient command of mode create"
+                    call parse_ori_vec(orient_string, orient(i,:))
                     arg_pos = arg_pos+1
-                    ori_pos=2
-                    do j = 1,3
-                        if (orient_string(ori_pos:ori_pos) == '-') then 
-                            ori_pos = ori_pos + 1
-                            read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
-                            if (stat>0) STOP "Error reading orient value"
-                            orient(i,j) = -orient(i,j)
-                            ori_pos = ori_pos + 1
-                        else
-                            read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
-                            if(stat>0) STOP "Error reading orient value"
-                            ori_pos=ori_pos + 1
-                        end if
-                    end do
+
+                    ! ori_pos=2
+                    ! do j = 1,3
+                    !     if (orient_string(ori_pos:ori_pos) == '-') then 
+                    !         ori_pos = ori_pos + 1
+                    !         read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
+                    !         if (stat>0) STOP "Error reading orient value"
+                    !         orient(i,j) = -orient(i,j)
+                    !         ori_pos = ori_pos + 1
+                    !     else
+                    !         read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
+                    !         if(stat>0) STOP "Error reading orient value"
+                    !         ori_pos=ori_pos + 1
+                    !     end if
+                    ! end do
                 end do
 
 

src/opt_group.f90

@@ -237,7 +237,7 @@ module opt_group
                     !here as well to make sure they are in the box
                     do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3
                         call apply_periodic(r_interp(:,j))
-                        call add_atom(type_interp(j), r_interp(:,j))
+                        call add_atom(type_interp(j), sbox_ele(ie), r_interp(:,j))
                     end do 
                 end do
 

src/opt_orient.f90

@@ -0,0 +1,121 @@
+module opt_orient
+    !This module contains the orient option which allows for the reorientation 
+    !of simulation cells. This can be used to create arbitrarily oriented dislocation or loops.
+    use parameters
+    use subroutines
+    use elements
+    use box
+    implicit none
+
+
+    real(kind=dp), save :: new_orient(3,3)
+    real(kind=dp), dimension(6) :: orig_box_bd
+    real(kind=dp), allocatable :: orig_sub_box_ori(:,:,:)
+
+    public 
+    contains
+
+    subroutine orient(arg_pos)
+
+        integer, intent(inout) :: arg_pos
+
+        integer :: i, ibasis, inod
+        real(kind=dp) :: inv_sub_box_ori(3,3,sub_box_num)
+
+        !First parse the orient command
+        call parse_orient(arg_pos)
+
+        !Now rotate the basis. To do this we transform the basis to [100] [010] [001] and then 
+        !transform to user specified basis.
+
+        !Find all inverse orientation matrices for all sub_boxes
+        do i = 1, sub_box_num
+            call matrix_inverse(sub_box_ori, 3, inv_sub_box_ori)
+        end do
+
+        !Now transform all atoms
+        do i = 1, atom_num
+            r_atom(:,i) = matmul(new_orient,matmul(inv_sub_box_ori(:,:,sbox_atom(i)),r_atom(:,i)))
+        end do 
+
+        !Now transform all elements
+        do i = 1, ele_num
+            do inod =1, ng_node(lat_ele(i))
+                do ibasis = 1, basisnum(lat_ele(i))
+                    r_node(:,ibasis,inod,i) = matmul(new_orient,matmul(inv_sub_box_ori(:,:,sbox_ele(i)),r_node(:,ibasis,inod,i)))
+                end do
+            end do  
+        end do
+
+        !Now save the original sub_box_ori and overwrite them
+        if(allocated(orig_sub_box_ori)) deallocate(orig_sub_box_ori)
+
+        allocate(orig_sub_box_ori(3,3,sub_box_num))
+        orig_sub_box_ori = sub_box_ori
+
+        !Now overwrite the orientations
+        do i = 1, sub_box_num
+            sub_box_ori(:,:,i) = new_orient
+        end do 
+
+        !Save original box boundaries
+        orig_box_bd = box_bd
+
+        !Now find new box boundaries
+        call def_new_box
+    end subroutine orient
+
+    subroutine parse_orient(arg_pos)
+        !This parses the orient option 
+        integer, intent(inout) :: arg_pos
+
+        integer :: i, arg_len
+        character(len=8) :: ori_string
+
+        !Pull out the new user orientation
+        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,:))
+        end do
+
+        !Normalize the orientation matrix
+        new_orient = matrix_normal(new_orient,3)
+        arg_pos = arg_pos + 1
+
+    end subroutine parse_orient
+
+    subroutine unorient
+
+        integer :: i, ibasis, inod
+        real(kind=dp) :: inv_ori(3,3)
+
+        !Now rotate the basis. To do this we transform the basis to [100] [010] [001] and then 
+        !transform to the original sbox_ele
+
+        !Find the inverse for the new orientation matrix
+        call matrix_inverse(new_orient, 3, inv_ori)
+
+        !Recover original sub_box_ori
+        sub_box_ori = orig_sub_box_ori 
+
+        !Now transform all atoms
+        do i = 1, atom_num
+            r_atom(:,i) = matmul(sub_box_ori(:,:,sbox_atom(i)),matmul(inv_ori(:,:),r_atom(:,i)))
+        end do 
+
+        !Now transform all elements
+        do i = 1, ele_num
+            do inod =1, ng_node(lat_ele(i))
+                do ibasis = 1, basisnum(lat_ele(i))
+                    r_node(:,ibasis,inod,i) = matmul(sub_box_ori(:,:,sbox_ele(i)),matmul(inv_ori,r_node(:,ibasis,inod,i)))
+                end do
+            end do  
+        end do
+
+        !Restore original box boundaries
+        box_bd = orig_box_bd
+    end subroutine unorient
+
+end module opt_orient