Current working version of insert group

src/main.f90

@@ -37,7 +37,7 @@ program main
     !Call initialization functions
     call lattice_init
     call box_init
-    call random_seed
+    call init_random_seed
     force_overwrite=.false.
     wrap_flag = .false.
 

src/opt_group.f90

@@ -8,9 +8,12 @@ module opt_group
     use box 
     implicit none
 
-    integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize
+    integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize, insert_type, &
-    character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape
+               insert_site
-    real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth, shell_thickness
+    character(len=15) :: type, gshape!Type indicates what element type is selected and shape is the group shape
+    real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth, shell_thickness, insert_conc, &
+    insert_lattice
+
     logical :: displace, delete, max_remesh, refine, group_nodes, flip, efill, refinefill
 
     integer, allocatable :: element_index(:), atom_index(:) 
@@ -29,6 +32,7 @@ module opt_group
         group_ele_num = 0
         group_atom_num = 0
         remesh_size=0
+        insert_type = 0
         random_num=0
         group_type=0
         notsize=0
@@ -75,6 +79,11 @@ module opt_group
             call change_group_type
         end if
 
+        if(insert_type > 0) then 
+            call get_group
+            call insert_group
+        end if
+
     end subroutine group
 
     subroutine parse_group(arg_pos)
@@ -98,11 +107,11 @@ module opt_group
         end select
 
         arg_pos = arg_pos + 1
-        call get_command_argument(arg_pos, shape, arglen)
+        call get_command_argument(arg_pos, gshape, arglen)
         if (arglen==0) STOP "Missing group_shape in group command"
 
         !Now parse the arguments required by the user selected shape
-        select case(trim(adjustl(shape)))
+        select case(trim(adjustl(gshape)))
         case('block')
             do i= 1, 6
                 arg_pos  = arg_pos + 1
@@ -146,7 +155,7 @@ module opt_group
                         vertices(i,2) = box_bd(2*i-1)
                         vertices(i,3) = box_bd(2*i-1)
                     else
-                        print *, "bwidth cannot be 0 in wedge shaped group"
+                        print *, "bwidth cannot be 0 in wedge gshaped group"
                         stop 3
                     end if
                 else if (i == dim2) then 
@@ -369,7 +378,7 @@ module opt_group
             continue
 
         case default
-            print *, "Group shape ", trim(adjustl(shape)), " is not currently accepted. Please check documentation ", &
+            print *, "Group shape ", trim(adjustl(gshape)), " is not currently accepted. Please check documentation ", &
                      "for accepted group shapes."
         end select
 
@@ -426,8 +435,34 @@ module opt_group
                 if(arglen ==0) stop "Missing notsize size"
                 read(textholder, *) notsize
                 print *, "Ignoring elements with size ", notsize
+            case('insert')
+                arg_pos=arg_pos+1
+                call get_command_argument(arg_pos, textholder, arglen)
+                if (arglen==0) stop "Missing element type for insert command"
+                call add_atom_type(textholder, insert_type)
+                arg_pos=arg_pos + 1 
+                call get_command_argument(arg_pos, textholder, arglen)
+                select case(trim(adjustl(textholder)))
+                case('tetra')
+                    insert_site=1
+                case('octa')
+                    insert_site=2
+                case('mixed')
+                    insert_site=3
+                case default
+                    print *, "site value must be tetra, octa, or mixed and not ", trim(adjustl(textholder))
+                    stop 3
+                end select
+                arg_pos=arg_pos+1
+                call get_command_argument(arg_pos, textholder, arglen)
+                if(arglen ==0) stop "Missing lattice_type in insert command"
+                read(textholder, *) insert_lattice
+                arg_pos=arg_pos+1
+                call get_command_argument(arg_pos, textholder, arglen)
+                if(arglen ==0) stop "Missing concentration in insert command"
+                read(textholder, *) insert_conc
             case default
-                !If it isn't an available option to opt_disl then we just exit
+                !If it isn't an available option to opt_group then we just exit
                 exit 
             end select
         end do
@@ -440,7 +475,7 @@ module opt_group
         integer, allocatable :: resize_array(:) 
         real(kind=dp) :: r_center(3), rand
 
-        select case(trim(adjustl(shape)))
+        select case(trim(adjustl(gshape)))
         case('block')
             print *, "Group has block shape with boundaries: ", block_bd
         case ('wedge')
@@ -1007,10 +1042,71 @@ module opt_group
 
     end subroutine change_group_type
     
+    subroutine insert_group
+        !This code inserts atoms into interstitial sites. This only works on atoms within the group due to the limitations with the
+        !Coarse-graining methodology which doesn't allow for concentration fluctuations.
+        real(kind=dp)  interstitial_sites(3,14), rand, rinsert(3)
+        integer :: add_num, i, j, rindex, sindex, ia, rlo, rhi, sbox
+        integer, allocatable :: used_sites(:,:)
+
+        !First save all of the displacement vectors from a lattice site to interstitial site
+        !The first 6 are the octohedral sites and the next 8 are the tetrahedral sites
+        interstitial_sites= reshape( (/ -0.5_dp, 0.0_dp, 0.0_dp, &
+                                         0.5_dp, 0.0_dp, 0.0_dp, &
+                                         0.0_dp,-0.5_dp, 0.0_dp, &
+                                         0.0_dp, 0.5_dp, 0.0_dp, &
+                                         0.0_dp, 0.0_dp,-0.5_dp, &
+                                         0.0_dp, 0.0_dp, 0.5_dp, &
+                                        -0.25_dp,-0.25_dp,-0.25_dp, &
+                                        -0.25_dp,-0.25_dp, 0.25_dp, &
+                                        -0.25_dp, 0.25_dp,-0.25_dp, &
+                                        -0.25_dp, 0.25_dp, 0.25_dp, &
+                                         0.25_dp,-0.25_dp,-0.25_dp, &
+                                         0.25_dp,-0.25_dp, 0.25_dp, &
+                                         0.25_dp, 0.25_dp,-0.25_dp, &
+                                         0.25_dp, 0.25_dp, 0.25_dp /), &
+                                         shape(interstitial_sites))
+
+        !First we calculate the number of atoms needed based on the number of atoms in the group and the concentration
+        interstitial_sites=interstitial_sites*insert_lattice
+
+        add_num = (insert_conc*group_atom_num)/(1-insert_conc)
+        allocate(used_sites(2,add_num))
+
+        print *, "Inserting ", add_num, " atoms as atom type ", insert_type
+
+        !Now set up the random number generator for the desired interstitial type
+        select case(insert_site)
+        case(1)
+            rlo=1
+            rhi=6
+        case(2)
+            rlo=7 
+            rhi = 14
+        case(3)
+            rlo=1
+            rhi=14
+        end select
 
-    subroutine split_group_elements
+        !Now add the atoms
-        !
+        i = 1
-    end subroutine split_group_elements
+        addloop:do while ( i < add_num)
+            call random_number(rand)
+            rindex = int(1+rand*(group_atom_num-1))
+            ia=atom_index(rindex)
+            call random_number(rand)
+            sindex = int(rlo+rand*(rhi-rlo))
+            do j = 1, i
+                if((ia == used_sites(1,i)).and.(sindex == used_sites(2,i))) cycle addloop
+            end do
+            rinsert = r_atom(:,ia) + matmul(sub_box_ori(:,:,sbox_atom(ia)),interstitial_sites(:,sindex))
+            sbox = sbox_atom(ia)
+            call add_atom(0, insert_type, sbox, rinsert)
+            used_sites(1,i) = ia
+            used_sites(2,i) = sindex
+            i = i + 1
+        end do addloop
+    end subroutine insert_group
 
     function in_group(r)
         !This subroutine determines if a point is within the group boundaries
@@ -1019,7 +1115,7 @@ module opt_group
 
         integer :: dim3, i
         logical :: in_group
-        select case(trim(adjustl(shape)))
+        select case(trim(adjustl(gshape)))
         case('block')
             in_group=in_block_bd(r,block_bd)
         case('wedge')
@@ -1068,7 +1164,7 @@ module opt_group
 
         in_group_ele=.false.
 
-        if(trim(adjustl(shape)) == 'shell') then 
+        if(trim(adjustl(gshape)) == 'shell') then 
             node_in_out(:) = -1
             !First calculate whether each element node is within the shell region, inside the shell sphere, or outside the
             !shell region
@@ -1083,13 +1179,13 @@ module opt_group
                     exit nodeloop 
                 end if
                 
-                shape ='sphere'
+                gshape ='sphere'
                 if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(esize/=notsize)) then 
                     node_in_out(inod) = 1
                 else
                     node_in_out(inod) = 0
                 end if
-                shape='shell'
+                gshape='shell'
             end do nodeloop
 
             !If any nodes are within the shell region, or if the shell region interescts an element then add it to the group

src/subroutines.f90

@@ -241,5 +241,52 @@ module subroutines
         return
     end subroutine check_right_ortho
 
+    subroutine init_random_seed()
+        implicit none
+        integer, allocatable :: seed(:)
+        integer :: i, n, un, istat, dt(8), pid, t(2), s
+        integer(8) :: count, tms
+      
+        call random_seed(size = n)
+        allocate(seed(n))
+        ! First try if the OS provides a random number generator
+        open(newunit=un, file="/dev/urandom", access="stream", &
+             form="unformatted", action="read", status="old", iostat=istat)
+        if (istat == 0) then
+           read(un) seed
+           close(un)
+        else
+           ! Fallback to XOR:ing the current time and pid. The PID is
+           ! useful in case one launches multiple instances of the same
+           ! program in parallel.
+           call system_clock(count)
+           if (count /= 0) then
+              t = transfer(count, t)
+           else
+              call date_and_time(values=dt)
+              tms = (dt(1) - 1970) * 365_8 * 24 * 60 * 60 * 1000 &
+                   + dt(2) * 31_8 * 24 * 60 * 60 * 1000 &
+                   + dt(3) * 24 * 60 * 60 * 60 * 1000 &
+                   + dt(5) * 60 * 60 * 1000 &
+                   + dt(6) * 60 * 1000 + dt(7) * 1000 &
+                   + dt(8)
+              t = transfer(tms, t)
+           end if
+           s = ieor(t(1), t(2))
+           pid = getpid() + 1099279 ! Add a prime
+           s = ieor(s, pid)
+           if (n >= 3) then
+              seed(1) = t(1) + 36269
+              seed(2) = t(2) + 72551
+              seed(3) = pid
+              if (n > 3) then
+                 seed(4:) = s + 37 * (/ (i, i = 0, n - 4) /)
+              end if
+           else
+              seed = s + 37 * (/ (i, i = 0, n - 1 ) /)
+           end if
+        end if
+        call random_seed(put=seed)
+      end subroutine init_random_seed
 
 end module subroutines