diff --git a/src/Makefile b/src/Makefile index d26a8e0..b31870b 100644 --- a/src/Makefile +++ b/src/Makefile @@ -33,7 +33,7 @@ $(OBJDIR)/%.o: %.f90 deps: - @fortdepend -o Makefile.dep -i mpi -b obj -w + @makedepf90 -b obj *.f90 > Makefile.dep #----------------- CLEAN UP -----------------# diff --git a/src/caller.f90 b/src/caller.f90 index a95fbb3..63221a2 100644 --- a/src/caller.f90 +++ b/src/caller.f90 @@ -60,6 +60,8 @@ module caller call group(arg_pos) case('-ow') arg_pos = arg_pos + 1 + case('-novel') + arg_pos = arg_pos+1 case('-wrap') arg_pos = arg_pos + 1 case('-norefine') diff --git a/src/elements.f90 b/src/elements.f90 index 8a14b86..eedcd95 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -437,7 +437,7 @@ module elements real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: r_interp !Interpolated atomic positions real(kind = dp), optional, intent(in) :: eng(:,:), f(:,:,:), & v(:,:,:) - real(kind=dp), dimension(10, max_basisnum*max_esize**3), optional,intent(out) :: data_interp !Interpolated atomic positions + real(kind=dp), dimension(:,:), optional,intent(out) :: data_interp !Interpolated atomic positions !Internal variables integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp, adjust @@ -528,6 +528,68 @@ module elements return end subroutine interpolate_atoms + subroutine interpolate_vel(type, esize, lat_type, vel_in, vel_interp) + !This subroutine returns the interpolated atoms from the elements. + + !Arguments + character(len=*), intent(in) :: type !The type of element that it is + integer, intent(in) :: esize !The number of atoms per side + integer, intent(in) :: lat_type !The integer lattice type of the element + real(kind=dp), dimension(:,:,:), intent(in) :: vel_in !Nodal positions + real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: vel_interp !Interpolated atomic positions + + !Internal variables + integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp, adjust + real(kind=dp) :: a_shape(max_ng_node) + real(kind=dp) :: t, s, r + + !Initialize some variables + vel_interp(:,:) = 0.0_dp + ia = 0 + + !Define bnum based on the input lattice type. If lat_type=0 then we are interpolating lattice points which means + !the basis is 0,0,0, and the type doesn't matter + + select case(lat_type) + case(0) + bnum=1 + lat_type_temp = 1 + case default + bnum = basisnum(lat_type) + lat_type_temp = lat_type + end select + + select case(type) + case('fcc','bcc') + !Now loop over all the possible sites + do it = 1, esize + t=-1.0_dp +(it-1)*(2.0_dp/(esize-1)) + do is =1, esize + s=-1.0_dp +(is-1)*(2.0_dp/(esize-1)) + do ir = 1, esize + r=-1.0_dp +(ir-1)*(2.0_dp/(esize-1)) + call rhombshape(r,s,t,a_shape) + + do ibasis = 1, bnum + ia = ia + 1 + do inod = 1, 8 + vel_interp(:,ia) = vel_interp(:,ia) + a_shape(inod) * vel_in(:,ibasis,inod) + + end do + end do + + end do + end do + end do + end select + if (ia /= bnum*esize**3) then + print *, "Incorrect interpolation" + stop 3 + end if + return + + end subroutine interpolate_vel + subroutine rhombshape(r,s,t, shape_fun) !Shape function for rhombohedral elements real(kind=8), intent(in) :: r, s, t @@ -661,6 +723,7 @@ module elements end do end subroutine wrap_atoms + subroutine wrap_elements integer :: i, inod, ibasis, j, node_in_bd(3,max_ng_node) real(kind =dp) :: box_len(3) @@ -1018,11 +1081,36 @@ do i = 1, atom_num subroutine add_atom_data(ia, eng, force, virial) !Function which sets the atom data arrays integer, intent(in) :: ia - real(kind=dp), intent(in) :: eng, force(3), virial(6) + real(kind=dp), intent(in) :: eng, force(:), virial(:) + + integer :: size_atom_array + real(kind=dp), allocatable :: tmp_eng(:), tmp_force(:,:), tmp_virial(:,:) + + size_atom_array=size(tag_atom) + if(ia > atom_num) then + stop "ia in add_atom_dat shouldn't be larger than atom_num" + else if(ia > size(energy_atom)) then + allocate(tmp_eng(size(energy_atom)+1024)) + allocate(tmp_force(3, size(energy_atom)+1024)) + allocate(tmp_virial(6, size(energy_atom)+1024)) + + tmp_eng=0.0_dp + tmp_eng(1:size(energy_atom)) = energy_atom + call move_alloc(tmp_eng, energy_atom) + + tmp_force=0.0_dp + tmp_force(:, 1:size(force_atom, 2)) = force_atom + call move_alloc(tmp_force, force_atom) + + tmp_virial=0.0_dp + tmp_virial(:, 1:size(virial_atom,2)) = virial_atom + call move_alloc(tmp_virial, virial_atom) + end if energy_atom(ia) = eng force_atom(:,ia) = force(:) virial_atom(:,ia) = virial(:) + vflag = .true. return end subroutine add_atom_data diff --git a/src/functions.f90 b/src/functions.f90 index e0bdd16..8d120ad 100644 --- a/src/functions.f90 +++ b/src/functions.f90 @@ -168,11 +168,11 @@ END FUNCTION StrDnCase do i =1 ,3 !Check upper bound - if(v(i) > (box_bd(2*i)+10.0_dp**(-10)) ) then + if(v(i) > (box_bd(2*i))) then in_block_bd =.false. exit !Check lower bound - else if (v(i) < (box_bd(2*i-1)-10.0_dp**(-10))) then + else if (v(i) < (box_bd(2*i-1))+1d-6) then in_block_bd = .false. exit end if diff --git a/src/io.f90 b/src/io.f90 index 098cd42..f8640fb 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -139,7 +139,7 @@ module io do i = 1, ele_num do inod = 1, ng_node(lat_ele(i)) do ibasis = 1, basisnum(lat_ele(i)) - write(11, '(2i16, 3f23.15)') basis_type(ibasis,lat_ele(i)), 1, r_node(:,ibasis,inod,i) + write(11, *) basis_type(ibasis,lat_ele(i)), 1, r_node(:,ibasis,inod,i) outn = outn + 1 end do end do @@ -151,7 +151,7 @@ module io !Write atom positions do i = 1, atom_num - write(11, '(2i16, 3f23.15)') type_atom(i), 0, r_atom(:,i) + write(11, *) type_atom(i), 0, r_atom(:,i) outn = outn + 1 end do @@ -207,7 +207,7 @@ module io write(11, '(a)') 'Atoms' write(11, '(a)') ' ' do i = 1, atom_num - write(11, '(2i16, 3f23.15)') tag_atom(i), type_atom(i), r_atom(:,i) + write(11, *) tag_atom(i), type_atom(i), r_atom(:,i) end do !Write refined element atomic positions @@ -219,7 +219,7 @@ module io do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 interp_num = interp_num+1 call apply_periodic(r_interp(:,iatom)) - write(11, '(2i16, 3f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom) + write(11, *) atom_num+interp_num, type_interp(iatom), r_interp(:,iatom) end do end select end do @@ -277,12 +277,12 @@ module io if (write_dat) then do i = 1, atom_num - write(11, '(2i16, 13f23.15)') tag_atom(i), type_atom(i), r_atom(:,i), energy_atom(i), & + write(11, *) tag_atom(i), type_atom(i), r_atom(:,i), energy_atom(i), & force_atom(:,i), virial_atom(:,i) end do else do i = 1, atom_num - write(11, '(2i16, 3f23.15)') tag_atom(i), type_atom(i), r_atom(:,i) + write(11, *) tag_atom(i), type_atom(i), r_atom(:,i) end do end if @@ -294,12 +294,12 @@ module io do inod =1, ng_node(lat_ele(i)) do ibasis=1, basisnum(lat_ele(i)) if(write_dat) then - write(11, '(2i16, 13f23.15)') atom_num+interp_num, basis_type(ibasis,lat_ele(i)), & + write(11, *) atom_num+interp_num, basis_type(ibasis,lat_ele(i)), & r_node(:, ibasis, inod, i), energy_node(ibasis,inod,i), force_node(:, ibasis, inod, i), & virial_node(:, ibasis, inod, i) else - write(11, '(2i16, 3f23.15)') atom_num+interp_num, basis_type(ibasis,lat_ele(i)), & + write(11, *) atom_num+interp_num, basis_type(ibasis,lat_ele(i)), & r_node(:, ibasis, inod, i) end if interp_num = interp_num +1 @@ -321,14 +321,14 @@ module io do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 interp_num = interp_num+1 call apply_periodic(r_interp(:,iatom)) - write(11, '(2i16, 13f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom), & + write(11, *) atom_num+interp_num, type_interp(iatom), r_interp(:,iatom), & data_interp(:,iatom) end do else do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 interp_num = interp_num+1 call apply_periodic(r_interp(:,iatom)) - write(11, '(2i16, 3f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom) + write(11, *) atom_num+interp_num, type_interp(iatom), r_interp(:,iatom) end do end if end select @@ -998,7 +998,7 @@ module io !Internal Variables integer :: i, j, in_eles, in_atoms, inbtypes(10), lat_type, ia, ie, inod, & - id, type_node, ilat, esize, tag, type, bnum, n, ibasis, ip, atom_type_map(100) + id, type_node, ilat, esize, tag, atype, bnum, n, ibasis, ip, atom_type_map(100) real(kind=dp) :: newdisplace(3), ra(3), in_lapa, ea, fa(3), va(6), vela(3),& ee(10,20), fe(3,10,20), ve(6,10,20), re(3,10,20), vele(3,10,20), atomic_masses(10) character(len=100) :: textholder, fcc @@ -1064,11 +1064,11 @@ module io do ia = 1, in_atoms read(11,'(a)') line(:) if(read_vel) then - read(line,*) tag, type, ra(:), ea, fa(:), va(:), vela(:) + read(line,*) tag, atype, ra(:), ea, fa(:), va(:), vela(:) else - read(line,*) tag, type, ra(:), ea, fa(:), va(:) + read(line,*) tag, atype, ra(:), ea, fa(:), va(:) end if - call add_atom(tag, atom_type_map(type), ra) + call add_atom(tag, atom_type_map(atype), ra) call add_atom_data(atom_num, ea, fa, va) if(read_vel) vel_atom(:,atom_num) = vela end do diff --git a/src/mode_calc.f90 b/src/mode_calc.f90 index 25086a3..007fd17 100644 --- a/src/mode_calc.f90 +++ b/src/mode_calc.f90 @@ -25,6 +25,9 @@ module mode_calc case('tot_virial') allocate(calculated(6)) call calc_tot_virial + case('tot_energy') + allocate(calculated(1)) + call calc_tot_energy case default print *, trim(adjustl(calc_opt)), " is not accepted as a calc option in mode_calc" stop 3 @@ -87,5 +90,26 @@ module mode_calc print *, "Total virial is calculated as : (v11, v22, v33, v32, v31, v21)" print *, calculated - end subroutine + end subroutine calc_tot_virial + + subroutine calc_tot_energy + integer :: i, inod, ibasis, j + calculated=0 + !Sum atom energies + do i = 1, atom_num + calculated(1) = calculated(1) + energy_atom(i) + end do + + !Sum element energies + do i =1, ele_num + do inod=1, ng_node(lat_ele(i)) + do ibasis=1, basisnum(lat_ele(i)) + calculated(1) = calculated(1) + energy_node(ibasis, inod, i)*(size_ele(i)**3)/ng_node(lat_ele(i)) + end do + end do + end do + + print *, 'Total energy in eV is:' + print *, calculated(1) + end subroutine calc_tot_energy end module mode_calc diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 1f21393..23c5d02 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -12,9 +12,9 @@ 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) + orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3), basis_pos(3,10) integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), & - basis_pos(3,10), esize_nums, esize_index(10) + esize_nums, esize_index(10) logical :: dup_flag, dim_flag, efill, crossb(3) real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:) diff --git a/src/neighbors.f90 b/src/neighbors.f90 index c415eae..ca59185 100644 --- a/src/neighbors.f90 +++ b/src/neighbors.f90 @@ -82,7 +82,7 @@ module neighbors real(kind=dp), intent(in) :: rc_off integer, intent(in) :: n - real(kind=dp), dimension(3,n) :: r_list + real(kind=dp), dimension(:, :), intent(in) :: r_list integer :: i, j, c(3),cn(3), ci, cj, ck, num_nei, nei, v(3), period_dir(3) !Variables for cell list code @@ -92,7 +92,7 @@ module neighbors real(kind=dp) :: r(3), box_len(3) logical :: period_bd(3) - !First reallocate the neighbor list codes + !First reallocate the neighbor list variables if (allocated(nei_list)) then deallocate(nei_list,nei_num, periodvec) end if diff --git a/src/opt_bubble.f90 b/src/opt_bubble.f90 index 57d3fa9..c507979 100644 --- a/src/opt_bubble.f90 +++ b/src/opt_bubble.f90 @@ -32,7 +32,7 @@ module opt_bubble type='all' gshape='sphere' group_nodes = .true. - + group_atom_types=0 call get_group call refine_group call get_group @@ -80,6 +80,7 @@ module opt_bubble integer, intent(inout) :: arg_pos integer :: i, arglen + real(kind=dp) :: mass character(len=100) :: tmptxt @@ -111,7 +112,34 @@ module opt_bubble print *, species if(arglen == 0) stop "Missing bubble species" - arg_pos = arg_pos +1 + + !OPtional arguments + do while(.true.) + if(arg_pos > command_argument_count()) exit + arg_pos=arg_pos+1 + + call get_command_argument(arg_pos, tmptxt) + tmptxt=adjustl(tmptxt) + select case(trim(tmptxt)) + + case('excludetypes') + arg_pos=arg_pos + 1 + call get_command_argument(arg_pos, tmptxt, arglen) + if(arglen==0) stop "Missing number of atoms to exclude" + read(tmptxt, *) exclude_num + do i=1,exclude_num + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, tmptxt, arglen) + if(arglen==0) stop "Missing exclude atom types" + call atommass(tmptxt, mass) + call add_atom_type(mass, exclude_types(i)) + end do + + case default + exit + end select + end do + return end subroutine parse_bubble end module opt_bubble diff --git a/src/opt_delete.f90 b/src/opt_delete.f90 index e693510..35463cb 100644 --- a/src/opt_delete.f90 +++ b/src/opt_delete.f90 @@ -87,7 +87,7 @@ module opt_delete allocate(which_cell(3,atom_num)) !First pass the atom list and atom num to the algorithm which builds the cell list - call build_cell_list(atom_num, r_atom, rc_off, cell_num, num_in_cell, cell_list, which_cell) + call build_cell_list(atom_num, r_atom, 5*rc_off, cell_num, num_in_cell, cell_list, which_cell) !Now loop over every atom and figure out if it has neighbors within the rc_off delete_num = 0 diff --git a/src/opt_disl.f90 b/src/opt_disl.f90 index 7e5d13e..c88073d 100644 --- a/src/opt_disl.f90 +++ b/src/opt_disl.f90 @@ -520,7 +520,7 @@ module opt_disl end do return - end subroutine + end subroutine disloop !******************************************************** ! SOLIDANGLE diff --git a/src/opt_group.f90 b/src/opt_group.f90 index b8d6fe9..f3babae 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -12,7 +12,7 @@ module opt_group implicit none integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize, insert_type, & - insert_site, num_species + insert_site, num_species, group_atom_types, exclude_num, exclude_types(10) character(len=15) :: type, gshape!Type indicates what element type is selected and shape is the group shape character(len=2) :: species_type(10) real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth, shell_thickness, insert_conc, & @@ -39,6 +39,7 @@ module opt_group insert_type = 0 random_num=0 group_type=0 + group_atom_types=0 notsize=0 displace=.false. delete=.false. @@ -47,6 +48,7 @@ module opt_group flip = .false. group_nodes=.false. num_species = 0 + exclude_num = 0 if(allocated(element_index)) deallocate(element_index) if(allocated(atom_index)) deallocate(atom_index) @@ -87,10 +89,10 @@ module opt_group call displace_group end if if(insert_type > 0) then - print *, "Insert command has been dropped" - stop 1 + !print *, "Insert command has been dropped" + !stop 1 call get_group -! call insert_group + call insert_group end if if(num_species > 0) then @@ -470,6 +472,24 @@ module opt_group refine=.true. case('refinefill') refinefill=.true. + case('selecttype') + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) stop "Missing select type" + call atommass(textholder, mass) + call add_atom_type(mass, group_atom_types) + case('excludetypes') + arg_pos=arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if(arglen==0) stop "Missing number of atoms to exclude" + read(textholder, *) exclude_num + do i=1,exclude_num + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if(arglen==0) stop "Missing exclude atom types" + call atommass(textholder, mass) + call add_atom_type(mass, exclude_types(i)) + end do case('remesh') arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) @@ -619,7 +639,7 @@ module opt_group select case(trim(adjustl(type))) case('atoms','all') do i = 1, atom_num - if(in_group(r_atom(:,i)).neqv.flip) then + if(in_group(type_atom(i), r_atom(:,i)).neqv.flip) then group_atom_num = group_atom_num + 1 if (group_atom_num > size(atom_index)) then allocate(resize_array(size(atom_index) + 1024)) @@ -687,8 +707,19 @@ module opt_group !This command is used to refine the group to full atomistics integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num - real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3) + real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), data_interp(10, max_basisnum*max_esize**3), & + vel_interp(3, max_basisnum*max_esize**3) + real(kind=dp), allocatable :: vel_tmp(:,:) + logical :: refine_dat + + + if (allocated(force_node)) then + refine_dat=.true. + else + refine_dat=.false. + end if + print *, size(data_interp) !Refining to atoms if(group_ele_num > 0) then orig_atom_num = atom_num @@ -698,17 +729,38 @@ module opt_group do i = 1, group_ele_num ie = element_index(i) !Get the interpolated atom positions - call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) + if (refine_dat) then + call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp, & + energy_node(:,:,i), force_node(:,:,:,i), virial_node(:,:,:,i), data_interp) + else + call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) + end if + if(allocated(vel_atom)) then + call interpolate_vel(type_ele(i), size_ele(i), lat_ele(i), vel_node(:,:,:,i), vel_interp) + end if !Loop over all interpolated atoms and add them to the system, we apply periodic boundaries !here as well to make sure they are in the box do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3 call apply_periodic(r_interp(:,j)) call add_atom(0,type_interp(j), r_interp(:,j)) + if(refine_dat) then + call add_atom_data(atom_num, data_interp(1, j), data_interp(2:4, j), data_interp(5:10, j)) + end if + if(allocated(vel_atom)) then + if(size(vel_atom, 2) < size(type_atom)) then + allocate(vel_tmp(3, size(type_atom))) + vel_tmp=0.0_dp + vel_tmp(:, 1:size(vel_atom,2)) = vel_atom + call move_alloc(vel_tmp, vel_atom) + end if + vel_atom(:, atom_num) = vel_interp(:,j) + end if end do end do !Once all atoms are added we delete all of the elements call delete_elements(group_ele_num, element_index) print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms." + end if end subroutine refine_group @@ -745,7 +797,7 @@ module opt_group ravg = ravg + rfill(:,ibasis, 1)/basisnum(lat_ele(ie)) end do - if( in_group(ravg)) then + if( in_group(0,ravg)) then nump_ele = nump_ele - 1 end if end do @@ -1118,77 +1170,79 @@ module opt_group end subroutine change_group_type -! subroutine insert_group -! !This code inserts atoms into interstitial sites. This only works on atoms within the group due to the limitations with the -! !Coarse-graining methodology which doesn't allow for concentration fluctuations. -! real(kind=dp) interstitial_sites(3,14), rand, rinsert(3) -! integer :: add_num, i, j, rindex, sindex, ia, rlo, rhi -! integer, allocatable :: used_sites(:,:) -! -! !First save all of the displacement vectors from a lattice site to interstitial site -! !The first 6 are the octohedral sites and the next 8 are the tetrahedral sites -! interstitial_sites= reshape( (/ -0.5_dp, 0.0_dp, 0.0_dp, & -! 0.5_dp, 0.0_dp, 0.0_dp, & -! 0.0_dp,-0.5_dp, 0.0_dp, & -! 0.0_dp, 0.5_dp, 0.0_dp, & -! 0.0_dp, 0.0_dp,-0.5_dp, & -! 0.0_dp, 0.0_dp, 0.5_dp, & -! -0.25_dp,-0.25_dp,-0.25_dp, & -! -0.25_dp,-0.25_dp, 0.25_dp, & -! -0.25_dp, 0.25_dp,-0.25_dp, & -! -0.25_dp, 0.25_dp, 0.25_dp, & -! 0.25_dp,-0.25_dp,-0.25_dp, & -! 0.25_dp,-0.25_dp, 0.25_dp, & -! 0.25_dp, 0.25_dp,-0.25_dp, & -! 0.25_dp, 0.25_dp, 0.25_dp /), & -! shape(interstitial_sites)) -! -! !First we calculate the number of atoms needed based on the number of atoms in the group and the concentration -! interstitial_sites=interstitial_sites*insert_lattice -! -! add_num = (insert_conc*group_atom_num)/(1-insert_conc) -! allocate(used_sites(2,add_num)) -! -! print *, "Inserting ", add_num, " atoms as atom type ", insert_type -! -! !Now set up the random number generator for the desired interstitial type -! select case(insert_site) -! case(1) -! rlo=1 -! rhi=6 -! case(2) -! rlo=7 -! rhi = 14 -! case(3) -! rlo=1 -! rhi=14 -! end select -! -! !Now add the atoms -! i = 1 -! addloop:do while ( i < add_num) -! call random_number(rand) -! rindex = int(1+rand*(group_atom_num-1)) -! ia=atom_index(rindex) -! call random_number(rand) -! sindex = int(rlo+rand*(rhi-rlo)) -! do j = 1, i -! if((ia == used_sites(1,i)).and.(sindex == used_sites(2,i))) cycle addloop -! end do -! rinsert = r_atom(:,ia) + matmul(sub_box_ori(:,:,sbox_atom(ia)),interstitial_sites(:,sindex)) -! sbox = sbox_atom(ia) -! call add_atom(0, insert_type, sbox, rinsert) -! used_sites(1,i) = ia -! used_sites(2,i) = sindex -! i = i + 1 -! end do addloop -! end subroutine insert_group + subroutine insert_group + !This code inserts atoms into interstitial sites. This only works on atoms within the group due to the limitations with the + !Coarse-graining methodology which doesn't allow for concentration fluctuations. + real(kind=dp) interstitial_sites(3,14), rand, rinsert(3) + integer :: add_num, i, j, rindex, sindex, ia, rlo, rhi + integer, allocatable :: used_sites(:,:) + + !First save all of the displacement vectors from a lattice site to interstitial site + !The first 6 are the octohedral sites and the next 8 are the tetrahedral sites + interstitial_sites= reshape( (/ -0.5_dp, 0.0_dp, 0.0_dp, & + 0.5_dp, 0.0_dp, 0.0_dp, & + 0.0_dp,-0.5_dp, 0.0_dp, & + 0.0_dp, 0.5_dp, 0.0_dp, & + 0.0_dp, 0.0_dp,-0.5_dp, & + 0.0_dp, 0.0_dp, 0.5_dp, & + -0.25_dp,-0.25_dp,-0.25_dp, & + -0.25_dp,-0.25_dp, 0.25_dp, & + -0.25_dp, 0.25_dp,-0.25_dp, & + -0.25_dp, 0.25_dp, 0.25_dp, & + 0.25_dp,-0.25_dp,-0.25_dp, & + 0.25_dp,-0.25_dp, 0.25_dp, & + 0.25_dp, 0.25_dp,-0.25_dp, & + 0.25_dp, 0.25_dp, 0.25_dp /), & + shape(interstitial_sites)) + + !First we calculate the number of atoms needed based on the number of atoms in the group and the concentration + interstitial_sites=interstitial_sites*insert_lattice + + + !Now set up the random number generator for the desired interstitial type + select case(insert_site) + case(1) + rlo=1 + rhi=6 + case(2) + rlo=7 + rhi = 14 + case(3) + rlo=1 + rhi=14 + end select + + if ((insert_conc > 1).or.(insert_conc < 0)) stop "Insert concentration should be between 0 and 1" + add_num = (insert_conc*group_atom_num) + allocate(used_sites(2,add_num)) + + print *, "Inserting ", add_num, " atoms as atom type ", insert_type + + !Now add the atoms + i = 1 + addloop:do while ( i < add_num) + call random_number(rand) + rindex = int(1+rand*(group_atom_num-1)) + ia=atom_index(rindex) + call random_number(rand) + sindex = int(rlo+rand*(rhi-rlo)) + do j = 1, i + if((ia == used_sites(1,i)).and.(sindex == used_sites(2,i))) cycle addloop + end do + rinsert = r_atom(:,ia) + matmul(box_ori,interstitial_sites(:,sindex)) + call add_atom(0, insert_type, rinsert) + used_sites(1,i) = ia + used_sites(2,i) = sindex + i = i + 1 + end do addloop + end subroutine insert_group subroutine alloy_group !This subroutine randomizes the atom types to reach desired concentrations, this only operates on atoms integer :: i, j, ia, type_map(10), added_types(num_species) real(kind=dp) :: rand, mass, old_fractions(num_species) character(len=2) :: sorted_species_type(10) + logical :: species_mapped(10) print *, "Alloying group with desired fractions", s_fractions(1:num_species) @@ -1196,10 +1250,13 @@ module opt_group !First we sort the fractions call quicksort(s_fractions(1:num_species)) + species_mapped = .false. do i = 1,num_species do j=1, num_species - if (is_equal(s_fractions(i), old_fractions(j))) then + if (is_equal(s_fractions(i), old_fractions(j)) .and. .not. species_mapped(j)) then sorted_species_type(i)=species_type(j) + species_mapped(j) = .true. + exit end if end do end do @@ -1230,8 +1287,9 @@ module opt_group return end subroutine alloy_group - function in_group(r) + function in_group(atype, r) !This subroutine determines if a point is within the group boundaries + integer, intent(in) :: atype real(kind=dp), intent(in) :: r(3) real(kind=dp) :: rnorm @@ -1276,6 +1334,19 @@ module opt_group case('all') in_group = .true. end select + + if(group_atom_types> 0) then + in_group = in_group.and.(atype== group_atom_types) + elseif(exclude_num > 0) then + if(in_group) then + do i=1,exclude_num + if (atype == exclude_types(i)) then + in_group=.false. + exit + end if + end do + end if + end if end function in_group function in_group_ele(esize, elat, rn) @@ -1299,13 +1370,13 @@ module opt_group r_center(:)= r_center(:) + rn(:,ibasis,inod)/basisnum(elat) end do - if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(size_ele(i)/=notsize)) then + if((in_group(0,rn(:, ibasis, inod)).neqv.flip).and.(size_ele(i)/=notsize)) then node_in_out(inod) = 2 exit nodeloop end if gshape ='sphere' - if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(esize/=notsize)) then + if((in_group(0,rn(:, ibasis, inod)).neqv.flip).and.(esize/=notsize)) then node_in_out(inod) = 1 else node_in_out(inod) = 0 @@ -1327,7 +1398,7 @@ module opt_group end do end do - if ((in_group(r_center).neqv.flip).and.(esize/= notsize)) then + if ((in_group(0,r_center).neqv.flip).and.(esize/= notsize)) then in_group_ele=.true. return end if @@ -1336,7 +1407,7 @@ module opt_group r_center(:) = 0.0_dp do inod = 1, ng_node(elat) do ibasis = 1, basisnum(elat) - if ((in_group(rn(:,ibasis,inod)).neqv.flip).and.(esize/=notsize)) then + if ((in_group(0,rn(:,ibasis,inod)).neqv.flip).and.(esize/=notsize)) then in_group_ele=.true. return end if diff --git a/src/opt_slip_plane.f90 b/src/opt_slip_plane.f90 index e000460..88b9409 100644 --- a/src/opt_slip_plane.f90 +++ b/src/opt_slip_plane.f90 @@ -17,11 +17,13 @@ module opt_slip_plane !Main calling function for the slip_plane option integer, intent(inout) :: arg_pos - integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3), ibasis + integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3), ibasis,& + starting_anum integer, allocatable :: slip_eles(:), temp_int(:) real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node), ratom(3,max_basisnum), & - maxp, minp + maxp, minp, vel_interp(3, max_basisnum*max_esize**3) + real(kind=dp), allocatable :: vel_tmp(:,:) integer :: type_interp(max_basisnum*max_esize**3) logical :: lat_points(max_esize,max_esize, max_esize) @@ -58,11 +60,25 @@ module opt_slip_plane !If we aren't efilling then just refine the element if(.not.efill) then + starting_anum=atom_num call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3 call apply_periodic(r_interp(:,ia)) call add_atom(0, type_interp(ia), r_interp(:,ia)) end do + + if(allocated(vel_atom)) then + call interpolate_vel(type_ele(ie), size_ele(ie), lat_ele(ie), vel_node(:,:,:,ie), vel_interp) + if(size(vel_atom,2) < size(type_atom)) then + allocate(vel_tmp(3, size(type_atom))) + vel_tmp=0.0_dp + vel_tmp(:, 1:size(vel_atom,2)) = vel_atom + call move_alloc(vel_tmp, vel_atom) + end if + do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3 + vel_atom(:, starting_anum+ia) = vel_interp(:, ia) + end do + end if !If we are efilling then the code is slightly more complex else !First populate the lat points array