From acd902db4ba0c72fe3c5ab99132ba95bcbd381fa Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Tue, 21 Apr 2020 13:25:33 -0400 Subject: [PATCH 1/2] working changes to bcc files --- src/elements.f90 | 27 ++++++++++++++++++++++++--- src/mode_create.f90 | 40 ++++++++++++++++++++++++++++++---------- 2 files changed, 54 insertions(+), 13 deletions(-) diff --git a/src/elements.f90 b/src/elements.f90 index e9fc928..5240640 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -27,7 +27,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. @@ -86,7 +86,26 @@ 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 /), & + 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 @@ -300,6 +319,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) @@ -345,7 +366,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 @@ -632,4 +653,4 @@ module elements end select end subroutine -end module elements \ No newline at end of file +end module elements diff --git a/src/mode_create.f90 b/src/mode_create.f90 index ffdfc05..545cea4 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -63,10 +63,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 - box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1 + !Now get the rotated box vertex positions in lattice space. Should be integer units and get the new maxlen + 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) @@ -83,7 +90,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 @@ -117,8 +129,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" @@ -270,13 +283,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.& - (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))))))& - then + 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,3)), abs(orient(i,1))).and.(is_equal(abs(orient(i,2)),2.0_dp*abs(orient(i,3))))))& + 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 @@ -518,4 +538,4 @@ module mode_create end subroutine error_message -end module mode_create \ No newline at end of file +end module mode_create From 3a59b23be786278c6bca02ec0de2f27a240849d0 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sat, 25 Apr 2020 14:42:57 -0400 Subject: [PATCH 2/2] Working bcc crystal structure when viewing in .lmp or .vtk format --- src/elements.f90 | 14 +++++++++----- src/io.f90 | 6 +++++- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/src/elements.f90 b/src/elements.f90 index 5240640..7e7f1c5 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -92,14 +92,15 @@ module elements 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.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.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 /), & + 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)) @@ -148,6 +149,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 diff --git a/src/io.f90 b/src/io.f90 index f929e29..e46f941 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -162,7 +162,9 @@ module io write_num = atom_num do i = 1,ele_num 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 @@ -196,7 +198,7 @@ module io do i = 1, ele_num call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp) select case(trim(adjustl(type_ele(i)))) - case('fcc') + case('fcc','bcc') do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 interp_num = interp_num+1 call apply_periodic(r_interp(:,iatom)) @@ -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)