diff --git a/src/io.f90 b/src/io.f90 index e5b098d..99f9e96 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -194,6 +194,7 @@ module io select case(trim(adjustl(type_ele(i)))) case('fcc') do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 + call apply_periodic(r_interp(:,iatom)) write(11, '(2i16, 3f23.15)') atom_num+iatom, type_interp(iatom), r_interp(:,iatom) end do end select diff --git a/src/mode_merge.f90 b/src/mode_merge.f90 index cfb4273..567b19b 100644 --- a/src/mode_merge.f90 +++ b/src/mode_merge.f90 @@ -8,7 +8,9 @@ module mode_merge use elements character(len=4) :: dim - integer :: in_num + integer :: in_num, new_starts(2) + real(kind=dp) :: shift_vec(3) + logical :: wrap, shift_flag public contains @@ -18,12 +20,23 @@ module mode_merge integer :: i real(kind=dp) :: displace(3), temp_box_bd(6) + wrap = .false. + shift_flag = .false. + + shift_vec(:) = 0.0_dp !First we parse the merge command call parse_command(arg_pos) !Now loop over all files and stack them do i = 1, in_num displace(:) = 0.0_dp + + !The new starts variable dictate where in the atom and element array each new + !file starts. This is used for additional options that can be applied to solely + !these new atoms/elements that are read in. + new_starts(1) = atom_num + 1 + new_starts(2) = ele_num + 1 + if ((i==1).or.(trim(adjustl(dim)) == 'none')) then call read_in(i, displace, temp_box_bd) call grow_box(temp_box_bd) @@ -40,6 +53,8 @@ module mode_merge call read_in(i, displace, temp_box_bd) call grow_box(temp_box_bd) end if + + if(shift_flag) call shift(new_starts, i) end do return @@ -89,7 +104,16 @@ module mode_merge !Choose what to based on what the option string is select case(trim(textholder)) - + case('shift') + shift_flag = .true. + do i = 1,3 + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) stop "Missing vector component for shift command" + read(textholder, *) shift_vec(i) + end do + case('wrap') + wrap = .true. case default !If it isn't an available option to mode merge then we just exit exit @@ -97,4 +121,52 @@ module mode_merge end do end subroutine parse_command + + subroutine shift(array_start, filenum) + !This subroutine applies a shift to newly added atoms and elements. It also wraps the atoms + !if the user provides the wrap flag + + integer, dimension(2), intent(in) :: array_start + integer, intent(in) :: filenum + + integer :: i, j, ibasis, inod + real(kind=dp), dimension(3) :: current_shift + + + !Calculate the current shift which is the filenum-1 multiplied by the user specified shift + current_shift = (filenum-1)*shift_vec + + !First shift all the atoms + do i = array_start(1), atom_num + r_atom(:,i) = r_atom(:,i) + current_shift + end do + + !Now shift all the elements + do i = array_start(2), ele_num + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + r_node(:,ibasis, inod, i) = r_node(:,ibasis, inod, i) + current_shift + end do + end do + end do + + !Now we check if we have to wrap the atoms, nodes are not wrapped. For elements the periodic + !boundary conditions are applied in the actual CAC codes + if(wrap) then + do i = array_start(1), atom_num + call apply_periodic(r_atom(:,i)) + end do + + !If we don't include the wrap command then we have to increase the size of the box + else + do i = 1,3 + if (current_shift(i) < -lim_zero) then + box_bd(2*i-1) = box_bd(2*i-1) - current_shift(i) + else if (current_shift(i) > lim_zero) then + box_bd(2*i) = box_bd(2*i) + current_shift(i) + end if + end do + end if + + end subroutine shift end module mode_merge \ No newline at end of file diff --git a/src/subroutines.f90 b/src/subroutines.f90 index 330fbca..7085df1 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -249,4 +249,23 @@ module subroutines end do end subroutine siftdown + + + subroutine apply_periodic(r) + !This function checks if an atom is outside the box and wraps it back in. This is generally used when some + !kind of displacement is applied but the simulation cell is desired to be maintained as the same size. + + real(kind=dp), dimension(3), intent(inout) :: r + + integer :: j + + do j = 1, 3 + if (r(j) > box_bd(2*j)) then + r(j) = r(j) - box_bd(2*j) + else if (r(j) < box_bd(2*j-1)) then + r(j) = r(j) + box_bd(2*j-1) + end if + end do + end subroutine + end module subroutines \ No newline at end of file