Working version of slip_plane ooption with effill

development
Alex Selimov 4 years ago
parent b9ce916e42
commit 84a84e578a

@ -20,7 +20,8 @@ module opt_slip_plane
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, 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)
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node), ratom(3,max_basisnum), &
maxp, minp
integer :: type_interp(max_basisnum*max_esize**3)
logical :: lat_points(max_esize,max_esize, max_esize)
@ -77,23 +78,29 @@ module opt_slip_plane
ele = cubic_cell*(esize -1)
do o = 1, size_ele(ie) - esize
do n = 1, size_ele(ie) - esize
do m = 1, size_ele(ie) - esize
latloop:do m = 1, size_ele(ie) - esize
do inod = 1, ng_node(lat_ele(ie))
vlat = ele(:,inod) + (/ m, n, o /)
if (.not.lat_points(vlat(1), vlat(2),vlat(3))) cycle latloop
call get_interp_pos(vlat(1), vlat(2), vlat(3), ie, rfill(:,:,inod))
end do
!Check to make sure all lattice points exist for the current element
if(any(lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1) == 0)) cycle latloop
!Check to see if the plane intersects this element if not then add it
if((spos < maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))).and. &
(spos > minval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)))))) then
maxp = maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))
minp = minval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))
if(.not.(spos < maxp).and.(spos > minp))then
nump_ele = nump_ele - esize**3
lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1) = .false.
call add_element(0, type_ele(ie), esize, lat_ele(ie), sbox_ele(ie), rfill)
new_ele_num = new_ele_num + 1
end if
end do
end do latloop
end do
end do
end if
esize= esize-2
end do
! Now add the leftover lattice points as atoms
do o = 1, size_ele(ie)
@ -102,6 +109,7 @@ module opt_slip_plane
if(lat_points(m,n,o)) then
call get_interp_pos(m,n,o, ie, ratom(:,:))
do ibasis = 1, basisnum(lat_ele(ie))
call apply_periodic(r_atom(:,ibasis))
call add_atom(0, basis_type(ibasis,lat_ele(ie)), sbox_ele(ie), ratom(:,ibasis))
end do
end if

Loading…
Cancel
Save