module mode_create !This mode is intended for creating element/atom regions and writing them to specific files. use parameters use atoms use io use subroutines use elements implicit none character(len=100) :: name, element_type real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), & orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), adjustVar(3,8) integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6) logical :: dup_flag, dim_flag real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:) public contains subroutine create() ! Main subroutine which controls execution character(len=100) :: textholder integer :: i, ibasis, inod real(kind=dp), allocatable :: r_node_temp(:,:,:) !Initialize default parameters orient = reshape((/ 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp /), shape(orient)) cell_mat(:,:)=0.0_dp name ='' element_type = '' esize=0 lattice_parameter=0.0_dp duplicate(:) = 0 box_len(:) = 0.0_dp dup_flag = .false. dim_flag = .false. basisnum = 0 lat_ele_num = 0 lat_atom_num = 0 !First we parse the command call parse_command() ! Before building do a check on the file if (outfilenum == 0) then textholder = 'none' call get_out_file(textholder) end if !Now we setup the unit element and call other init subroutines call def_ng_node(1, element_type) allocate(r_node_temp(3,max_basisnum,max_ng_node)) if(dup_flag) then !We initialize the cell with a lattice_parameter of 1 because we will add the lattice parameter later call cell_init(1.0_dp, esize, element_type, orient, cell_mat) do i = 1, 8 box_vert(:,i) = duplicate(:)*esize*lattice_space(:)*cubic_cell(:,i) + origin(:) end do call matrix_inverse(orient,3,orient_inv) !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 !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) end do !and then 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 default print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", & "element type" stop 3 end select !Now that it is multiply by the lattice parameter box_bd = box_bd*lattice_parameter else if(dim_flag) then continue else if(lmpcac) then do inod = 1, 8 do i = 1,3 if(is_equal(cubic_cell(i, inod),0.0_dp)) then adjustVar(i,inod) = -0.5_dp else adjustVar(i, inod) = 0.5_dp end if end do end do else adjustVar(:,:)=0.0_dp end if ! call cell_init(lattice_parameter, esize, element_type, orient, cell_mat) !If the user doesn't pass any build instructions than we just put the cell mat into the element_array call alloc_ele_arrays(1,0) !Add the basis atoms to the unit cell do inod = 1, max_ng_node do ibasis = 1, basisnum(1) r_node_temp(:,ibasis,inod) = lattice_parameter*matmul(orient, & matmul(fcc_mat, (esize+1)*cubic_cell(:,inod)+adjustVar(:,inod))) & + basis_pos(:,ibasis,1) end do end do do i = 1,3 box_bd(2*i) = maxval(r_node_temp(i,:,:)) box_bd(2*i-1) = origin(i) end do call add_element(element_type, esize, 1, r_node_temp) end if !If we passed the dup_flag or dim_flag then we have to convert the lattice points and add them to the atom/element arrays if(dup_flag.or.dim_flag) then !Allocate variables call alloc_ele_arrays(lat_ele_num, lat_atom_num*basisnum(1)) 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,1)) end do end do deallocate(r_atom_lat) end if if(lat_ele_num > 0) then do i = 1, lat_ele_num do inod= 1, ng_node(1) do ibasis = 1, basisnum(1) r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis,1) end do end do call add_element(element_type, esize, 1, r_node_temp) end do end if end if end subroutine create !This subroutine parses the command and pulls out information needed for mode_create subroutine parse_command() integer :: arg_pos, ori_pos, i, j, arglen, stat character(len=100) :: textholder character(len=8) :: orient_string !Pull out all required positional arguments call get_command_argument(2, name, arglen) if(arglen==0) STOP "Name is missing in mode create" call get_command_argument(3,element_type, arglen) if(arglen==0) STOP "Element_type is missing in mode create" call get_command_argument(4,textholder, arglen) if(arglen==0) STOP "Lattice Parameter is missing in mode create" read(textholder, *, iostat=stat) lattice_parameter if(stat > 0) STOP "Error reading lattice parameter" call get_command_argument(5, textholder, arglen) if(arglen==0) STOP "Esize missing in mode create" read(textholder, *, iostat=stat) esize if(stat > 0) STOP "Error reading esize" arg_pos = 6 !Check for optional keywords do while(.true.) if(arg_pos > command_argument_count()) exit !Pull out the next argument which should either be a keyword or an option call get_command_argument(arg_pos, textholder) textholder=adjustl(textholder) arg_pos=arg_pos+1 !Choose what to based on what the option string is select case(trim(textholder)) !If orient command is passed extract the orientation to numeric array format case('orient') 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" 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 end do !If the duplicate command is passed then we extract the information on the new bounds. case('duplicate') dup_flag = .true. do i = 1, 3 call get_command_argument(arg_pos, textholder) read(textholder, *) duplicate(i) arg_pos = arg_pos + 1 end do case('origin') do i = 1, 3 call get_command_argument(arg_pos, textholder) read(textholder, *) origin(i) arg_pos = arg_pos + 1 end do !If a filetype is passed then we add name.ext to the outfiles list case('xyz') textholder = trim(adjustl(name)) //'.xyz' call get_out_file(textholder) case default !Check to see if it is an option command, if so then mode_create must be finished if(textholder(1:1) == '-') then exit !Check to see if a filename was passed elseif(scan(textholder,'.',.true.) > 0) then call get_out_file(textholder) end if end select end do !Calculate the lattice periodicity length in lattice units do i = 1, 3 lattice_space(i) = norm2(orient(i,:)) end do !Check special periodicity relations select case(trim(adjustl(element_type))) case('fcc') do i = 1,3 !Check if one of the directions is 110 if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(orient(i,3),0.0_dp))).or.& (is_equal(abs(orient(i,2)), abs(orient(i,3))).and.(is_equal(orient(i,1),0.0_dp))).or.& (is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(orient(i,2),0.0_dp)))) then 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 lattice_space(i) = 0.5_dp * lattice_space(i) end if end do end select !Now normalize the orientation matrix orient = matrix_normal(orient,3) !If we haven't defined a basis then define the basis and add the default basis atom type and position if (basisnum(1) == 0) then max_basisnum = 1 basisnum(1) = 1 call add_atom_type(name, basis_type(1,1)) !If basis command not defined then we use name as the atom_name basis_pos(:,1,1) = 0.0_dp end if end subroutine subroutine build_with_rhomb(box_in_lat, transform_matrix) !This subroutine returns all the lattice points in the box in r_lat !Inputs integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space !Internal variables integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, templatpoints, ele(3,8), rzero(3), ilat, & type_interp(basisnum(1)*esize**3), vlat(3), temp_lat(3,8), m, n, o real(kind=dp) :: v(3), temp_nodes(3,1,8), ele_atoms(3,esize**3), r_interp(3,basisnum(1)*esize**3), adjustVar(3,8) real(kind=dp), allocatable :: resize_lat_array(:,:) logical, allocatable :: lat_points(:,:,:) logical :: node_in_bd(8) !Do some value initialization max_esize = esize !If we are writing out to lammpscac format then we have to adjust the nodal positions if(lmpcac) then do inod = 1, 8 do i = 1,3 if(is_equal(cubic_cell(i, inod),0.0_dp)) then adjustVar(i,inod) = -0.25_dp else adjustVar(i, inod) = 0.25_dp end if end do end do else adjustVar(:,:)=0.0_dp end if !First find the bounding lattice points (min and max points for the box in each dimension) numlatpoints = 1 do i = 1, 3 bd_in_lat(2*i-1) = minval(box_in_lat(i,:)) bd_in_lat(2*i) = maxval(box_in_lat(i,:)) numlatpoints = numlatpoints*(bd_in_lat(2*i)-bd_in_lat(2*i-1)) end do !Allocate the correct lat variable select case(esize) !Atomistics case(2) allocate(r_atom_lat(3,numlatpoints)) case default continue end select !Loop over all of lattice points within the boundary, we choose between two loops. One for the atomistic case !and one for the regular case if (esize==2) then !atomistics do iz = bd_in_lat(5)-5, bd_in_lat(6)+5 do iy = bd_in_lat(3)-5, bd_in_lat(4)+5 do ix = bd_in_lat(1)-5, bd_in_lat(2)+5 v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /) !Transform point back to real space for easier checking v = matmul(orient, matmul(transform_matrix,v)) !If within the boundaries if(in_block_bd(v, box_bd)) then lat_atom_num = lat_atom_num + 1 r_atom_lat(:,lat_atom_num) = v end if end do end do end do else !If we are working with elements we have to use more complex code !Initialize finite element ele(:,:) = (esize-1) * cubic_cell(:,:) !Make a 3 dimensional array of lattice points. This array is indexed by the integer lattice position. !A value of true means that the coordinate is a lattice point which is within the box_bd allocate(lat_points(bd_in_lat(2)-bd_in_lat(1)+10,bd_in_lat(4)-bd_in_lat(3)+10,bd_in_lat(6)-bd_in_lat(5)+10)) lat_points(:,:,:) = .false. do iz = bd_in_lat(5)-5, bd_in_lat(6)+5 do iy = bd_in_lat(3)-5, bd_in_lat(4)+5 do ix = bd_in_lat(1)-5, bd_in_lat(2)+5 v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /) !Transform point back to real space for easier checking v = matmul(orient, matmul(transform_matrix,v)) !If within the boundaries if(in_block_bd(v, box_bd)) then lat_points(ix-bd_in_lat(1)+5,iy-bd_in_lat(3)+5,iz-bd_in_lat(5) + 5) = .true. end if end do end do end do !Now we redefine bd_in_lat The first 3 indices contains limits for the lat_points array bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10 bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10 bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10 !Figure out where the starting point is. This is the first piont which fully contains the finite element outerloop: do iz = 1, bd_in_array(3) do iy = 1, bd_in_array(2) do ix = 1, bd_in_array(1) node_in_bd(8) = .false. do inod = 1, 8 vlat = ele(:,inod) + (/ ix, iy, iz /) !Check to see if the node resides at a position containing a lattice point within the box if(any(vlat > shape(lat_points))) then continue else if(lat_points(vlat(1),vlat(2),vlat(3))) then node_in_bd(inod) = .true. end if end do if(all(node_in_bd)) then rzero = (/ ix, iy, iz /) exit outerloop end if end do end do end do outerloop !Now build the finite element region lat_ele_num = 0 lat_atom_num = 0 allocate(r_lat(3,8,numlatpoints/esize)) !Redefined the second 3 indices as the number of elements that fit in the bounds do i = 1, 3 bd_in_array(3+i) = bd_in_array(i)/esize end do !Now start the element at rzero do inod=1, 8 ele(:,inod) = ele(:,inod) + rzero end do do iz = -bd_in_array(6), bd_in_array(6) do iy = -bd_in_array(5), bd_in_array(5) do ix = -bd_in_array(4), bd_in_array(4) node_in_bd(:) = .false. temp_nodes(:,:,:) = 0.0_dp temp_lat(:,:) = 0 do inod = 1, 8 vlat= ele(:,inod) + (/ ix*(esize), iy*(esize), iz*(esize) /) !Transform point back to real space for easier checking ! v = matmul(orient, matmul(transform_matrix,v)) do i = 1,3 v(i) = real(vlat(i) + bd_in_lat(2*i-1) - 5) end do temp_nodes(:,1, inod) = matmul(orient, matmul(transform_matrix, v)) temp_lat(:,inod) = vlat !Check to see if the lattice point values are greater then the array limits if(any(vlat > shape(lat_points)).or.any(vlat < 1)) then continue !If within array boundaries check to see if it is a lattice point else if(lat_points(vlat(1),vlat(2),vlat(3))) then node_in_bd(inod) = .true. end if end do if(all(node_in_bd)) then lat_ele_num = lat_ele_num+1 do inod = 1, 8 r_lat(:,inod,lat_ele_num) = temp_nodes(:,1,inod) + adjustVar(:,inod) end do !Now set all the lattice points contained within an element to false do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:)) do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:)) do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:)) lat_points(m,n,o) = .false. end do end do end do end if end do end do end do !Now figure out how many lattice points could not be contained in elements allocate(r_atom_lat(3,count(lat_points))) lat_atom_num = 0 do ix = 1, bd_in_array(3) do iy = 1, bd_in_array(2) do iz = 1, bd_in_array(1) !If this point is a lattice point then save the lattice point as an atom if (lat_points(ix,iy,iz)) then v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /) do i = 1,3 v(i) = v(i) + real(bd_in_lat(2*i-1) - 5) end do !Transform point back to real space v = matmul(orient, matmul(transform_matrix,v)) lat_atom_num = lat_atom_num + 1 r_atom_lat(:,lat_atom_num) = v end if end do end do end do end if end subroutine build_with_rhomb subroutine error_message(errorid) integer, intent(in) :: errorid select case(errorid) case(1) STOP "Name is missing." case(2) print *, "Element_type is missing" case(3) print *, "Lattice_parameter is missing" case(4) print *, "Lattice parameter is not in float form" case(5) print *, "Esize is missing" case(6) print *, "Esize is not in integer form" case(7) print *, "Cx(1) should equal Cx(2) in plane centroid finding algorithm. Please double check implementation" end select STOP 3 end subroutine error_message end module mode_create