diff --git a/src/elements.f90 b/src/elements.f90 index 6b41e57..b4e59be 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -9,7 +9,7 @@ module elements !Data structures used to represent the CAC elements. Each index represents an element character(len=100), allocatable :: type_ele(:) !Element type - integer, allocatable :: size_ele(:), lat_ele(:) !Element siz + integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:) !Element size real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array integer, save :: ele_num !Number of elements @@ -39,6 +39,7 @@ module elements !This can be easily increased with no change to efficiency integer :: max_basisnum, basisnum(10) !Max basis atom number, number of basis atoms in each lattice type integer :: basis_type(10,10) + real(kind=dp) :: lapa(10) public contains @@ -140,7 +141,7 @@ module elements !Allocate element arrays if(n > 0) then - allocate(type_ele(n), size_ele(n), lat_ele(n), r_node(3,max_basisnum, max_ng_node,n), & + allocate(type_ele(n), size_ele(n), lat_ele(n), sbox_ele(n), r_node(3,max_basisnum, max_ng_node,n), & stat=allostat) if(allostat > 0) then print *, "Error allocating element arrays in elements.f90 because of: ", allostat @@ -179,12 +180,17 @@ module elements allocate(temp_int(n+ele_num+buffer_size)) temp_int(1:ele_size) = lat_ele temp_int(ele_size+1:) = 0 - call move_alloc(temp_int(1:ele_size), lat_ele) + call move_alloc(temp_int, lat_ele) allocate(temp_int(n+ele_num+buffer_size)) temp_int(1:ele_size) = size_ele temp_int(ele_size+1:) = 0 - call move_alloc(temp_int(1:ele_size), size_ele) + call move_alloc(temp_int, size_ele) + + allocate(temp_int(n+ele_num+buffer_size)) + temp_int(1:ele_size) = lat_ele + temp_int(ele_size+1:) = 0 + call move_alloc(temp_int, sbox_ele) allocate(char_temp(n+ele_num+buffer_size)) char_temp(1:ele_size) = type_ele @@ -210,9 +216,9 @@ module elements end if end subroutine - subroutine add_element(type, size, lat, r) + subroutine add_element(type, size, lat, sbox, r) !Subroutine which adds an element to the element arrays - integer, intent(in) :: size, lat + integer, intent(in) :: size, lat, sbox character(len=100), intent(in) :: type real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node) @@ -220,6 +226,7 @@ module elements type_ele(ele_num) = type size_ele(ele_num) = size lat_ele(ele_num) = lat + sbox_ele(ele_num) = sbox r_node(:,:,:,ele_num) = r(:,:,:) node_num = node_num + ng_node(lat) diff --git a/src/io.f90 b/src/io.f90 index 4f9b74c..e5b098d 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -204,7 +204,7 @@ module io !This subroutine writes out a .lmp style dump file character(len=100), intent(in) :: file integer :: write_num, i, inod, ibasis - real(kind=dp) :: mass + real(kind=dp) :: mass, fcc_adjust(3,8), local_adjust(3,8), rout(3) 1 format(i16, ' Eight_Node', 4i16) 2 format(i16, ' Atom', 4i16) @@ -242,16 +242,32 @@ module io write(11, '(a)') 'CAC Elements' write(11, '(a)') ' ' + !Set up the nodal adjustment variables for all the different element types. This adjusts the node centers + !from the center of the unit cell (as formulated in this code) to the corners of the unit cells + do inod = 1, 8 + do i = 1,3 + if(is_equal(cubic_cell(i, inod),0.0_dp)) then + fcc_adjust(i,inod) = -0.5_dp + else + fcc_adjust(i, inod) = 0.5_dp + end if + end do + end do + fcc_adjust = matmul(fcc_mat, fcc_adjust) + !Write element nodal positions do i = 1, ele_num select case(trim(adjustl(type_ele(i)))) case('fcc') + !Now orient the current adjustment vector to the correct orientation + local_adjust = matmul(sub_box_ori(:,:,sbox_ele(i)), fcc_adjust) * lapa(lat_ele(i)) !The first entry is the element specifier write(11,1) i, basisnum(lat_ele(i)), size_ele(i), size_ele(i), size_ele(i) do ibasis = 1, basisnum(lat_ele(i)) do inod = 1, 8 !Nodal information for every node - write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i) + rout = r_node(:,ibasis,inod,i) + local_adjust(:,inod) + write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), rout end do end do end select @@ -486,6 +502,8 @@ module io write(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types) !Now for every lattice type write the basis atom types write(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = 1, lattice_types) + !Now for every lattice type write the lattice parameters + write(11,*) (lapa(i), i = 1, lattice_types) !Now write the numbers of elements and atoms write(11,*) atom_num, ele_num @@ -498,7 +516,7 @@ module io !Write out the elements, this is written in two stages, one line for the element and then 1 line for !every basis at every node do i = 1, ele_num - write(11, *) i, lat_ele(i), size_ele(i), type_ele(i) + write(11, *) i, lat_ele(i), size_ele(i), sbox_ele(i), type_ele(i) do inod = 1, ng_node(lat_ele(i)) do ibasis =1, basisnum(lat_ele(i)) write(11,*) inod, ibasis, r_node(:, ibasis, inod, i) @@ -582,7 +600,7 @@ module io real(kind = dp), dimension(6), intent(out) :: temp_box_bd integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles, new_atom_types, & - new_type_to_type(10), new_lattice_types + new_type_to_type(10), new_lattice_types, sbox character(len=100) :: etype real(kind=dp) :: r(3), newdisplace(3) real(kind=dp), allocatable :: r_innode(:,:,:) @@ -623,8 +641,6 @@ module io sub_box_array_bd(:,1,sub_box_num+1:) = sub_box_array_bd(:,1,sub_box_num+1:) + atom_num sub_box_array_bd(:,2,sub_box_num+1:) = sub_box_array_bd(:,2,sub_box_num+1:) + ele_num - sub_box_num = sub_box_num + n - !Read in the number of atom types and all their names read(11, *) new_atom_types, (new_type_to_name(i), i = 1, new_atom_types) !Now fit these into the global list of atom types, after this new_type_to_type is the actual global @@ -646,6 +662,8 @@ module io basis_type(i,j) = new_type_to_type(basis_type(i,j)) end do end do + !Read the lattice parameters for every lattice type + read(11,*) (lapa(i), i = lattice_types+1, lattice_types+new_lattice_types) !Read number of elements and atoms and allocate arrays read(11, *) in_atoms, in_eles call grow_ele_arrays(in_eles, in_atoms) @@ -659,7 +677,7 @@ module io !Read the elements do i = 1, in_eles - read(11, *) n, type, size, etype + read(11, *) n, type, size, sbox, etype do inod = 1, ng_node(type) do ibasis =1, basisnum(type) read(11,*) j, k, r_innode(:, ibasis, inod) @@ -667,7 +685,7 @@ module io end do end do type = type + lattice_types - call add_element(etype, size, type, r_innode) + call add_element(etype, size, type, sbox+n, r_innode) end do !Close the file being read @@ -676,5 +694,9 @@ module io !Only increment the lattice types if there are elements, if there are no elements then we !just overwrite the arrays if(in_eles > 0) lattice_types = lattice_types + new_lattice_types + + + sub_box_num = sub_box_num + n + end subroutine read_mb end module io diff --git a/src/main.f90 b/src/main.f90 index ccb4d31..ed34c53 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -28,17 +28,6 @@ program main ! Command line parsing arg_num = command_argument_count() - !First check if we are writing out to lammpscac format by looping over all arguments - do i = 1, arg_num - call get_command_argument(i, argument) - select case(argument(scan(argument,'.',.true.)+1:)) - case('cac') - lmpcac = .true. - case default - continue - end select - end do - !Determine if a mode is being used and what it is. The first argument has to be the mode !if a mode is being used call get_command_argument(1, argument) diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 9b2cd73..9463e92 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -12,7 +12,7 @@ module mode_create 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), duplicate(3), adjustVar(3,8) + orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3) integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), & basis_pos(3,10) logical :: dup_flag, dim_flag @@ -105,7 +105,7 @@ module mode_create 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) + call add_element(element_type, esize, 1, 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 @@ -141,7 +141,7 @@ module mode_create r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis) end do end do - call add_element(element_type, esize, 1, r_node_temp) + call add_element(element_type, esize, 1, 1, r_node_temp) end do end if end if @@ -199,7 +199,20 @@ module mode_create 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 - call parse_ori_vec(orient_string, orient(i,:)) + 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 @@ -234,6 +247,7 @@ module mode_create exit end select end do + !Calculate the lattice periodicity length in lattice units do i = 1, 3 lattice_space(i) = norm2(orient(i,:)) @@ -264,8 +278,9 @@ module mode_create !Now normalize the orientation matrix orient = matrix_normal(orient,3) - !Set lattice_num to 1 + !Set lattice_num to 1 and add the lattice_parameter to the elements module lattice paramter variable lattice_types = 1 + lapa(1) = lattice_parameter !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 @@ -286,7 +301,7 @@ module mode_create !Internal variables integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,8), rzero(3), & vlat(3), temp_lat(3,8), m, n, o - real(kind=dp) :: v(3), temp_nodes(3,1,8), adjustVar(3,8) + real(kind=dp) :: v(3), temp_nodes(3,1,8) real(kind=dp), allocatable :: resize_lat_array(:,:) logical, allocatable :: lat_points(:,:,:) logical :: node_in_bd(8) @@ -294,23 +309,6 @@ module mode_create !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.5_dp - else - adjustVar(i, inod) = 0.5_dp - end if - end do - end do - - adjustVar(:,1:8) = matmul(orient,matmul(fcc_mat,adjustVar(:,1:8))) - 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 @@ -328,6 +326,7 @@ module mode_create 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 @@ -440,11 +439,9 @@ module mode_create 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 - + lat_ele_num = lat_ele_num+1 + r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) + !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,:)) @@ -510,4 +507,4 @@ module mode_create end subroutine error_message -end module mode_create +end module mode_create \ No newline at end of file