Merge bcc to development branch

src/elements.f90

@@ -28,7 +28,7 @@ module elements
     integer :: atom_types = 0
 
     !Variables for creating elements based on primitive cells
-    real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3)
+    real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3), bcc_cell(3,8), bcc_mat(3,3), bcc_inv(3,3)
     integer :: cubic_faces(4,6)
 
     !Below are lattice type arrays which provide information on the general form of the elements.
@@ -87,7 +87,27 @@ module elements
                              0.0_dp, 0.5_dp, 0.5_dp, &
                              0.5_dp, 0.0_dp, 0.5_dp /), &
                              shape(fcc_mat))
+
+        !Initialize the bcc primitive cell
+        bcc_cell = reshape((/ 0.0_dp,  0.0_dp,  0.0_dp, &
+                              0.5_dp, -0.5_dp,  0.5_dp, &
+                              1.0_dp,  0.0_dp,  1.0_dp, &
+                              0.5_dp,  0.5_dp,  0.5_dp, &
+                             -0.5_dp,  0.5_dp,  0.5_dp, &
+                              0.0_dp,  0.0_dp,  1.0_dp, &
+                              0.5_dp,  0.5_dp,  1.5_dp, &
+                              0.0_dp,  1.0_dp,  1.0_dp /),  &
+                              shape(bcc_cell))
+
+        bcc_mat = reshape((/ 0.5_dp, -0.5_dp,  0.5_dp, &
+                             0.5_dp,  0.5_dp,  0.5_dp, &
+                            -0.5_dp,  0.5_dp,  0.5_dp /), &
+                             shape(bcc_mat))
+                             
+                              
+        
         call matrix_inverse(fcc_mat,3,fcc_inv)
+        call matrix_inverse(bcc_mat,3,bcc_inv)
 
         max_basisnum = 0
         basisnum(:) = 0
@@ -130,6 +150,9 @@ module elements
             end if
             cell_mat(:, 1:8) = fcc_cell + adjustVar(:,1:8)
             cell_mat(:,1:8) = lapa * ((esize-1)*matmul(orient_mat, cell_mat(:,1:8)))
+        case('bcc')
+            cell_mat(:,1:8) = bcc_cell
+            cell_mat(:,1:8) = lapa* ((esize-1)*matmul(orient_mat, cell_mat(:,1:8)))
         case default
             print *, "Element type ", trim(ele_type), " currently not accepted"
             stop
@@ -327,6 +350,8 @@ module elements
             select case(trim(adjustl(element_types(i))))
             case('fcc')
                 ng_node(i) = 8
+            case('bcc')
+                ng_node(i) = 8
             end select
 
             if(ng_node(i) > max_ng_node) max_ng_node = ng_node(i)
@@ -372,7 +397,7 @@ module elements
         end select
 
         select case(trim(adjustl(type)))
-        case('fcc')
+        case('fcc','bcc')
             allocate(a_shape(8))
             !Now loop over all the possible sites
             do it = 1, esize

src/io.f90

@@ -161,8 +161,10 @@ module io
         !Calculate total atom number
         write_num = atom_num
         do i = 1,ele_num
-            if(type_ele(i) == 'fcc') write_num = write_num + basisnum(lat_ele(i))*size_ele(i)**3
+            if(type_ele(i) == 'fcc') write_num = write_num + size_ele(i)**3
+            if(type_ele(i) == 'bcc') write_num = write_num + size_ele(i)**3
         end do
+
         !Write total number of atoms + elements
         write(11, '(i16, a)') write_num, ' atoms'
         !Write number of atom types 
@@ -339,6 +341,7 @@ module io
         end do
         close(11)
 
+        !Now we write the vtk file for the elements
         open(unit=11, file='cg_'//trim(adjustl(file)), action='write', status='replace',position='rewind')
         write(11,1)
         write(11,2)
@@ -357,6 +360,7 @@ module io
         write(11,5) ele_num
         do i = 1, ele_num
             if(trim(adjustl(type_ele(i))) == 'fcc') write(11, '(i16)') 12
+            if(trim(adjustl(type_ele(i))) == 'bcc') write(11, '(i16)') 12
         end do
         write(11,12) ele_num
         write(11,20)

src/mode_create.f90

@@ -64,10 +64,17 @@ module mode_create
             do i = 1, 8
                 box_vert(:,i) = duplicate(:)*esize*lattice_space(:)*cubic_cell(:,i) + (origin(:)/lattice_parameter) 
             end do
-            !Now get the rotated box vertex positions in lattice space. Should be integer units
+            !Now get the rotated box vertex positions in lattice space. Should be integer units and get the new maxlen
-            box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
+            select case(trim(adjustl(element_type)))
+            case('fcc')
+                box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
+                maxbd = maxval(matmul(orient,matmul(fcc_mat,box_lat_vert)),2)
+            case('bcc')
+                box_lat_vert = int(matmul(bcc_inv, matmul(orient_inv, box_vert)))+1
+                maxbd = maxval(matmul(orient,matmul(bcc_mat,box_lat_vert)),2)
+            end select
+                
             !Find the new maxlen
-            maxbd = maxval(matmul(orient,matmul(fcc_mat,box_lat_vert)),2)
             do i = 1, 3
                 box_bd(2*i) = maxval(box_vert(i,:)) - 0.25_dp*lattice_space(i)
                 box_bd(2*i-1) = origin(i)-0.25_dp*lattice_space(i)
@@ -84,7 +91,12 @@ module mode_create
                 box_vert(:,i) = (cubic_cell(:,i)*box_len(:) + origin(:))/lattice_parameter
             end do
             !Now get the rotated box vertex positions in lattice space. Should be integer units
-            box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
+            select case(trim(adjustl(element_type)))
+            case('fcc')
+                box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
+            case('bcc')
+                box_lat_vert = int(matmul(bcc_inv, matmul(orient_inv, box_vert)))+1
+            end select
 
             !Now get the box_bd in lattice units
             do i = 1, 3
@@ -118,8 +130,9 @@ module mode_create
             !Call the build function with the correct transformation matrix
             select case(trim(adjustl(element_type)))
             case('fcc')
-
                 call build_with_rhomb(box_lat_vert, fcc_mat)
+            case('bcc')
+                call build_with_rhomb(box_lat_vert, bcc_mat)
             case default
                 print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", &
                          "element type"
@@ -276,13 +289,20 @@ module mode_create
                     lattice_space(i) = 0.5_dp * lattice_space(i)
 
                 !Check if one direction is 112
-                else if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(abs(orient(i,3)),2.0_dp*abs(orient(i,1))))).or.&
+                else if((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(abs(orient(i,3)),2.0_dp*abs(orient(i,1))))).or.&
-                         (is_equal(abs(orient(i,2)), abs(orient(i,3))).and.(is_equal(abs(orient(i,1)),2.0_dp*abs(orient(i,2))))).or.&
+                        (is_equal(abs(orient(i,2)), abs(orient(i,3))).and.(is_equal(abs(orient(i,1)),2.0_dp*abs(orient(i,2))))).or.&
-                         (is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(abs(orient(i,2)),2.0_dp*abs(orient(i,3))))))&
+                        (is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(abs(orient(i,2)),2.0_dp*abs(orient(i,3))))))&
-                          then 
+                        then 
 
                     lattice_space(i) = 0.5_dp * lattice_space(i)
+                end if
+            end do
 
+        case('bcc')
+            do i = 1, 3
+                !Check if the direction is 111
+                if((is_equal(abs(orient(i,1)),abs(orient(i,2)))).and.(is_equal(abs(orient(i,2)),abs(orient(i,3))))) then 
+                    lattice_space(i) = 0.5_dp * lattice_space(i)
                 end if
             end do
         end select