Merge pull request #29 from aselimov/master

Accidental push to master

README.md

@@ -214,10 +214,10 @@ This command wraps atoms back into the simulation cell as though periodic bounda
 **Remesh**
 
 ```
-remesh esize 
+remesh esize lattice_parameter lattice_type
 ```
 
-This command remeshes the atoms/elements within the group to the new element size `esize`. Currently only accepts an `esize` of 2 which refines it to full atomistics. When remeshing to atomistics the group can contain any orientations of elements but when remeshing to different finite elements, the group must contain all atoms/elements with the same orientation.
+This command remeshes the atoms/elements within the group to the new element size `esize`. Currently only accepts an `esize` of 2 which refines it to full atomistics. When remeshing to atomistics the group can contain any orientations of elements but when remeshing to different finite elements, the group must contain all atoms/elements with the same orientation. `lattice_parameter` is the lattice parameter for the elements within the group and `lattice_type` is the lattice type (integer) that these new elements will be assigned to. 
 
 **Max**
 
@@ -291,3 +291,20 @@ Specifying positions in cacmb can be done through a variety of ways. Examples of
 
 `rand[val1:val2]` - returns a random position that lies within the range
 
+`rande[facenum]` - Returns a random position in an interelement boundary which is offset of the element face `facenum`. Face numbers are based on the which vertices comprise the face. Vertex numbers are shown in the figure below for the primitive fcc unit cell which is what the fcc rhombohedral element is based from. The face numbers are:
+
+Face 1: [1,2,3,4]
+
+Face 2: [1,2,6,5]
+
+Face 3: [2,3,7,6]
+
+Face 4: [3,4,8,7]
+
+Face 5: [1,4,8,5]
+
+Face 6: [5,6,7,8]
+
+Image for vertex numbers is: 
+
+![](/home/alexselimov/Documents/CACmb/Numbered_element.png)

src/elements.f90

@@ -64,11 +64,11 @@ module elements
                              shape(fcc_cell))
 
         !Now we create a list containing the list of vertices needed to describe the 6 cube faces
-        cubic_faces(:,1) = (/ 1, 4, 8, 5 /)
+        cubic_faces(:,1) = (/ 1, 2, 3, 4 /) 
-        cubic_faces(:,2) = (/ 2, 3, 7, 6 /)
+        cubic_faces(:,2) = (/ 1, 2, 6, 5 /)
-        cubic_faces(:,3) = (/ 1, 2, 6, 5 /)
+        cubic_faces(:,3) = (/ 2, 3, 7, 6 /)
         cubic_faces(:,4) = (/ 3, 4, 8, 7 /)
-        cubic_faces(:,5) = (/ 1, 2, 3, 4 /) 
+        cubic_faces(:,5) = (/ 1, 4, 8, 5 /)
         cubic_faces(:,6) = (/ 5, 6, 7, 8 /)
 
         !!Now initialize the fcc primitive cell
@@ -491,4 +491,144 @@ module elements
         end do 
     end subroutine
 
+    recursive subroutine parse_pos(i, pos_string, pos)
+        !This subroutine parses the pos command allowing for command which include inf
+        integer, intent(in) ::  i !The dimension of the position
+        character(len=100), intent(in) :: pos_string !The position string
+        real(kind=dp), intent(out) :: pos !The output parsed position value
+
+        integer :: iospara, face, ele, randsize
+        real(kind=dp) :: rand, rone, rtwo, rand_ele_pos(3)
+        character(len=100) :: cone, ctwo
+
+        iospara = 0
+        
+        if(trim(adjustl(pos_string)) == 'inf') then 
+            pos=box_bd(2*i)
+
+        else if(trim(adjustl(pos_string)) == '-inf') then
+            pos=box_bd(2*i-1)
+
+        else if (trim(adjustl(pos_string)) == 'rand') then 
+            call random_number(rand)
+            pos = (box_bd(2*i)-box_bd(2*i-1))*rand + box_bd(2*i-1)
+
+        else if (index(pos_string,'rande')>0) then 
+            !First select a random element 
+            call random_number(rand)
+            ele= 1 + floor(ele_num*rand)
+            !Now read the rest of the command which specifies the face we need and check 
+            !to make sure it's an accepted face number
+            cone = pos_string(index(pos_string, '[')+1:index(pos_string, '[')+1)
+            
+            !Check to see if we also pass an element size, if it was passed then make sure our 
+            !random element is the right size 
+            if(index(pos_string,':') > 0) then 
+                ctwo = pos_string(index(pos_string, ':')+1:index(pos_string,']')-1) 
+                read(ctwo, *) randsize
+                do while(randsize /= size_ele(ele))
+                    call random_number(rand)
+                    ele= 1 + floor(ele_num*rand)
+                end do 
+            end if
+
+            read(cone, *) face
+            if ((face < 1).or.(face > 6)) stop "Current face number must be 1,2,3,4,5,6. Please check documentation" 
+            !Now get the position
+            call offset_pos(ele, face, rand_ele_pos)
+            pos = rand_ele_pos(i)
+
+        else if (index(pos_string,'rand')>0) then 
+            call random_number(rand)
+            cone = pos_string(index(pos_string, '[')+1:index(pos_string,':')-1)
+            call parse_pos(i, cone, rone)
+            ctwo = pos_string(index(pos_string, ':')+1:index(pos_string,']')-1)
+            call parse_pos(i, ctwo, rtwo)
+            pos = (rtwo - rone)*rand + rone            
+
+        else if ((index(pos_string,'-') > 0).and.(index(pos_string,'inf')>0)) then 
+            !Now extract the number we are reducing from infinity
+            if(index(pos_string,'inf') < index(pos_string,'-')) then 
+                read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos
+            else
+                read(pos_string(1:index(pos_string,'-')-1), *, iostat=iospara) pos
+            end if
+            pos = box_bd(2*i) - pos
+
+        else if ((index(pos_string,'+') > 0).and.(index(pos_string,'inf')>0)) then 
+            !Now extract the number we are reducing from infinity
+           if(index(pos_string,'inf') < index(pos_string,'+')) then 
+                read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos
+            else
+                read(pos_string(1:index(pos_string,'+')-1), *, iostat=iospara) pos
+            end if
+            pos = box_bd(2*i-1) + pos
+
+        else if ((index(pos_string,'*') > 0).and.(index(pos_string,'inf')>0)) then 
+            !Now extract the number we are reducing from infinity
+           if(index(pos_string,'inf') < index(pos_string,'*')) then 
+                read(pos_string(index(pos_string,'*')+1:), *, iostat=iospara) pos
+            else
+                read(pos_string(1:index(pos_string,'*')-1), *, iostat=iospara) pos
+            end if
+            pos = (box_bd(2*i)-box_bd(2*i-1))*pos + box_bd(2*i-1)
+
+        else
+            read(pos_string, *, iostat=iospara) pos
+
+        end if
+
+        if (iospara > 0) then 
+            print *, "Error reading position argument ", trim(adjustl(pos_string)), ". Please reformat and try again."
+        end if
+    end subroutine parse_pos
+
+    subroutine offset_pos(ie, iface, pos)
+    !This returns a position slightly offset from the center of an element face
+    !This is used to return a random position within an element discontinuity
+
+    integer, intent(in) :: ie !Element index
+    integer, intent(in) :: iface !Face number between 1 and 6 
+    real(kind=dp), dimension(3), intent(out) :: pos !Position vector
+
+    !Other variables we need
+    integer :: esize
+    real(kind=dp) :: orient(3,3), ori_inv(3,3), lp, r_cubic_node(3,8)
+
+    !First find the offset vector for the face in the untransformed cubic cell
+    esize = size_ele(ie)
+    select case(iface)
+        case(1)
+            pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp /)
+        case(2)
+            pos = (/ real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /)
+        case(3)
+            pos = (/ (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /)
+        case(4)
+            pos = (/ real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /)
+        case(5)
+            pos = (/ -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /)
+        case(6)
+            pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp /)
+    end select
+
+    !Now transform it to real space and adjust it to the position of the element in the first node.
+    select case (trim(adjustl(type_ele(ie))))
+
+    case('fcc')
+        !First we have to extract the element lattice parameter
+        call matrix_inverse(sub_box_ori(:,:,sbox_ele(ie)),3,ori_inv)
+        r_cubic_node = r_node(:,1,:,ie)
+        r_cubic_node = matmul(fcc_inv,matmul(ori_inv,r_cubic_node))
+        lp = (maxval(r_cubic_node(1,:))-minval(r_cubic_node(1,:)))/(esize-1)
+        pos = matmul(sub_box_ori(:,:,sbox_ele(ie)),matmul(fcc_mat,pos))*lp
+        pos = pos + r_node(:,1, 1, ie)
+
+    case default
+        print *, trim(adjustl(type_ele(ie))), " is not a currently accepted element type for random element positions"
+        stop 3
+
+    end select
+    end subroutine
+
 end module elements

src/io.f90

@@ -373,7 +373,8 @@ module io
         !NOTE: This code doesn't work for arbitrary number of basis atoms per node. It assumes that the 
         !each element has only 1 atom type at the node.
         character(len=100), intent(in) :: file
-        integer :: interp_max, i, j, inod, ibasis, ip, unique_index(10), unique_num
+        integer :: interp_max, i, j, inod, ibasis, ip, unique_index(50), unique_size(50), unique_num, &
+                   etype
         real(kind=dp) :: box_vec(3)
 
 1  format('time' / i16, f23.15)
@@ -405,6 +406,22 @@ module io
 
         !Below writes the header information for the restart file
 
+
+            !First figure out all of the unique element types
+            unique_num = 0
+            unique_index(:) = 0
+            eleloop:do i = 1, ele_num
+                do j =1 , unique_num
+                    if ( ( size_ele(i) == size_ele( unique_index(j) ) ).and. &
+                        ( lat_ele(i) == lat_ele(unique_index(j)) ) ) then
+                        cycle eleloop
+                    end if
+                end do
+                unique_num = unique_num + 1
+                unique_index(unique_num) = i
+                unique_size(unique_num) = size_ele(i)
+            end do eleloop
+
         !Calculate the max number of atoms per element
         select case(max_ng_node)
         case(8)
@@ -415,7 +432,7 @@ module io
         write(11,20) interp_max
         write(11,3) node_num
         write(11,19) max_ng_node
-        write(11,4) lattice_types
+        write(11,4) unique_num
         write(11,5) atom_num
         write(11,6) 1 !Grain_num is ignored 
         write(11,16) lattice_types, atom_types
@@ -454,26 +471,21 @@ module io
         !write the element information
         if(ele_num > 0) then
             write(11,12)
-            !First figure out all of the unique element types
+           
-            unique_num = 0
-            unique_index(:) = 0
-            eleloop:do i = 1, ele_num
-                do j =1 , unique_num
-                    if ( ( size_ele(i) == size_ele( unique_index(j) ) ).and. &
-                        ( lat_ele(i) == lat_ele(unique_index(j)) ) ) then
-                        cycle eleloop
-                    end if
-                end do
-                unique_num = unique_num + 1
-                unique_index(unique_num) = i
-            end do eleloop
             do i = 1, unique_num
-                write(11,'(3i16)') i, size_ele(i)-1, basis_type(1,i)
+                write(11,'(3i16)') i, size_ele(unique_index(i))-1, basis_type(1,lat_ele(unique_index(i)))
             end do
             ip = 0
             write(11,13)
             do i = 1, ele_num
-                write(11, '(4i16)') i, lat_ele(i), 1, basis_type(1,lat_ele(i))
+                !Figure out the ele type
+                do j = 1, unique_num
+                    if ( unique_size(j) == size_ele(i)) then
+                        etype = j
+                        exit
+                    endif
+                end do
+                write(11, '(4i16)') 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
@@ -633,15 +645,16 @@ module io
 
         print *, "displace", displace
         do i = 1, 3
-            if (displace(i) > lim_zero) then 
+            if (abs(displace(i)) > lim_zero) then 
                 newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
             else
-                newdisplace=displace(i)
+                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 
        
+        print *, "newdisplace", newdisplace 
         call grow_box(temp_box_bd)
         !Read in the number of sub_boxes and allocate the variables
         read(11, *) n
@@ -730,7 +743,7 @@ module io
         do i = 1, in_atoms
             read(11,*) j, type, sbox, r(:)
             r = r+newdisplace
-            call add_atom(new_type_to_type(type), sbox, r)
+            call add_atom(new_type_to_type(type), sbox+sub_box_num, r)
         end do
 
         !Read the elements
@@ -742,7 +755,7 @@ module io
                     r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace
                 end do
             end do
-            call add_element(etype, size, new_lattice_map(type), sbox+n, r_innode)
+            call add_element(etype, size, new_lattice_map(type), sbox+sub_box_num, r_innode)
         end do
 
         !Close the file being read

src/opt_orient.f90

@@ -20,6 +20,7 @@ module opt_orient
         integer, intent(inout) :: arg_pos
 
         integer :: i, ibasis, inod
+        logical :: isortho, isrighthanded
         real(kind=dp) :: inv_sub_box_ori(3,3,sub_box_num)
 
         !First parse the orient command
@@ -70,6 +71,7 @@ module opt_orient
         integer, intent(inout) :: arg_pos
 
         integer :: i, arg_len
+        logical :: isortho, isrighthanded
         character(len=8) :: ori_string
 
         !Pull out the new user orientation
@@ -82,6 +84,15 @@ module opt_orient
 
         !Normalize the orientation matrix
         new_orient = matrix_normal(new_orient,3)
+
+        !Check right hand rule and orthogonality
+        call check_right_ortho(new_orient, isortho, isrighthanded)
+        if (.not.isortho) then 
+            stop "Directions in orient are not orthogonal"
+        else if (.not.isrighthanded) then 
+            stop "Directions in orient are not righthanded"
+        end if
+
         arg_pos = arg_pos + 1
 
     end subroutine parse_orient

src/subroutines.f90

@@ -172,66 +172,6 @@ module subroutines
         return
     end subroutine parse_ori_vec
     
-    recursive subroutine parse_pos(i, pos_string, pos)
-        !This subroutine parses the pos command allowing for command which include inf
-        integer, intent(in) ::  i !The dimension of the position
-        character(len=100), intent(in) :: pos_string !The position string
-        real(kind=dp), intent(out) :: pos !The output parsed position value
-
-        integer :: iospara
-        real(kind=dp) :: rand, rone, rtwo
-        character(len=100) :: cone, ctwo
-
-        iospara = 0
-        if(trim(adjustl(pos_string)) == 'inf') then 
-            pos=box_bd(2*i)
-        else if(trim(adjustl(pos_string)) == '-inf') then
-            pos=box_bd(2*i-1)
-
-        else if (trim(adjustl(pos_string)) == 'rand') then 
-            call random_number(rand)
-            pos = (box_bd(2*i)-box_bd(2*i-1))*rand + box_bd(2*i-1)
-
-        else if (index(pos_string,'rand')>0) then 
-            call random_number(rand)
-            cone = pos_string(index(pos_string, '[')+1:index(pos_string,':')-1)
-            call parse_pos(i, cone, rone)
-            ctwo = pos_string(index(pos_string, ':')+1:index(pos_string,']')-1)
-            call parse_pos(i, ctwo, rtwo)
-            pos = (rtwo - rone)*rand + rone            
-        else if ((index(pos_string,'-') > 0).and.(index(pos_string,'inf')>0)) then 
-            !Now extract the number we are reducing from infinity
-            if(index(pos_string,'inf') < index(pos_string,'-')) then 
-                read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos
-            else
-                read(pos_string(1:index(pos_string,'-')-1), *, iostat=iospara) pos
-            end if
-            pos = box_bd(2*i) - pos
-        else if ((index(pos_string,'+') > 0).and.(index(pos_string,'inf')>0)) then 
-            !Now extract the number we are reducing from infinity
-           if(index(pos_string,'inf') < index(pos_string,'+')) then 
-                read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos
-            else
-                read(pos_string(1:index(pos_string,'+')-1), *, iostat=iospara) pos
-            end if
-            pos = box_bd(2*i-1) + pos
-        else if ((index(pos_string,'*') > 0).and.(index(pos_string,'inf')>0)) then 
-            !Now extract the number we are reducing from infinity
-           if(index(pos_string,'inf') < index(pos_string,'*')) then 
-                read(pos_string(index(pos_string,'*')+1:), *, iostat=iospara) pos
-            else
-                read(pos_string(1:index(pos_string,'*')-1), *, iostat=iospara) pos
-            end if
-            pos = (box_bd(2*i)-box_bd(2*i-1))*pos + box_bd(2*i-1)
-
-        else
-            read(pos_string, *, iostat=iospara) pos
-        end if
-
-        if (iospara > 0) then 
-            print *, "Error reading position argument ", trim(adjustl(pos_string)), ". Please reformat and try again."
-        end if
-    end subroutine parse_pos
 
     subroutine heapsort(a)