aselimov
/
CACmb
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Current working changes

Browse Source
master
Alex Selimov 5 years ago
parent 9a484c86f6
commit babf4e176d
2 changed files with 226 additions and 165 deletions
Show Stats Download Patch File Download Diff File

49
src/io.f90
Unescape View File

@ -907,7 +907,6 @@ module io
!Read more useless data !Read more useless data
read(11,*) textholder read(11,*) textholder
read(11,*) textholder
!set max values and allocate variables !set max values and allocate variables
max_basisnum = maxval(basisnum) max_basisnum = maxval(basisnum)
@ -915,31 +914,39 @@ module io
call grow_ele_arrays(in_eles, in_atoms) call grow_ele_arrays(in_eles, in_atoms)
!Now start reading the elements !Now start reading the elements
do i = 1, in_eles if(in_eles > 0) then
read(11,*) j, etype, k, lat_type read(11,*) textholder
do inod = 1, 8 do i = 1, in_eles
read(11, *) j, k, r_in(:,1,inod) read(11,*) j, etype, k, lat_type
r_in(:,1,inod) = r_in(:,1,inod) + newdisplace do inod = 1, 8
read(11, *) j, k, r_in(:,1,inod)
r_in(:,1,inod) = r_in(:,1,inod) + newdisplace
end do
call add_element(in_lattype_map(lat_type), etype_map(etype), new_lattice_map(lat_type), sub_box_num + 1, r_in)
end do end do
call add_element(in_lattype_map(lat_type), etype_map(etype), new_lattice_map(lat_type), sub_box_num + 1, r_in) end if
end do
!Read useless data if(in_atoms > 0) then
read(11,*) textholder
read(11,*) textholder
do i = 1, in_atoms if (in_eles > 0) then
read(11,*) j, k, atom_type, r_in_atom(:) !Read useless data
r_in_atom = r_in_atom + newdisplace read(11,*) textholder
call add_atom(atom_type_map(atom_type), sub_box_num + 1, r_in_atom) read(11,*) textholder
end do end if
!Close file
close(11)
lattice_types = maxval(new_lattice_map) do i = 1, in_atoms
read(11,*) j, k, atom_type, r_in_atom(:)
r_in_atom = r_in_atom + newdisplace
call add_atom(atom_type_map(atom_type), sub_box_num + 1, r_in_atom)
end do
!Close file
close(11)
sub_box_num = sub_box_num + 1 lattice_types = maxval(new_lattice_map)
call set_max_esize sub_box_num = sub_box_num + 1
call set_max_esize
end if
end subroutine read_pycac end subroutine read_pycac
end module io end module io

328
src/opt_group.f90
Unescape View File

@ -8,9 +8,9 @@ module opt_group
use box use box
implicit none implicit none
integer :: group_ele_num, group_atom_num, remesh_size,normal, remesh_type, dim1, dim2 integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2
character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape
real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), remesh_lat_pam real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3)
logical :: displace, delete, max_remesh, refine logical :: displace, delete, max_remesh, refine
integer, allocatable :: element_index(:), atom_index(:) integer, allocatable :: element_index(:), atom_index(:)
@ -160,14 +160,6 @@ module opt_group
call get_command_argument(arg_pos, textholder, arglen) call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) stop "Missing remesh element size in group command" if (arglen==0) stop "Missing remesh element size in group command"
read(textholder, *) remesh_size read(textholder, *) remesh_size
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) stop "Missing remesh lattice parameter in group command"
read(textholder, *) remesh_lat_pam
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) stop "Missing remesh type in group command"
read(textholder, *) remesh_type
case('max') case('max')
max_remesh =.true. max_remesh =.true.
case('delete') case('delete')
@ -311,10 +303,10 @@ module opt_group
integer :: i, j, k, ix, iy, iz, inod, ibasis, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num, & integer :: i, j, k, ix, iy, iz, inod, ibasis, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num, &
current_esize, dof, max_lat(3), r_lat(3), ele(3,8), vlat(3), bd_in_lat(6), bd_in_array(3), old_ele, old_atom, & current_esize, dof, max_lat(3), r_lat(3), ele(3,8), vlat(3), bd_in_lat(6), bd_in_array(3), old_ele, old_atom, &
max_loops, working_esize max_loops, working_esize, group_lat_num, lat_list(10), group_sbox_num, sbox_list(100), is, ilat
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), ori_inv(3,3), r(3), & real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), ori_inv(3,3), r(3), &
r_new_node(3,max_basisnum, max_ng_node), orient(3,3), group_in_lat(3,8) r_new_node(3,max_basisnum, max_ng_node), orient(3,3), group_in_lat(3,8), group_bd(6)
logical, allocatable :: lat_points(:,:,:) logical, allocatable :: lat_points(:,:,:)
character(len=100) :: remesh_ele_type character(len=100) :: remesh_ele_type
@ -322,157 +314,219 @@ module opt_group
!Right now we just hardcode only remeshing to elements !Right now we just hardcode only remeshing to elements
remesh_ele_type = 'fcc' remesh_ele_type = 'fcc'
!Get the orientations, this assumes that the orientation of the subbox for the first atom is the ! Determine which sub_boxes and lattices types are within in the group
!same as the rest of the atoms group_sbox_num = 0
!If this assumption is false then the code will break and exit group_lat_num = 0
orient = sub_box_ori(:, :, sbox_atom(atom_index(1))) do i = 1, group_atom_num
call matrix_inverse(orient,3,ori_inv) do j = 1, group_sbox_num
if (sbox_list(j) == sbox_atom(atom_index(i))) exit
group_sbox_num = group_sbox_num + 1
sbox_list(group_sbox_num) = sbox_atom(atom_index(i))
end do
do j = 1, group_lat_num
if (basis_type(1,lat_list(j)) == type_atom(atom_index(i))) exit
group_lat_num = group_lat_num + 1
do k = 1, lattice_types
if (basis_type(1,k) == type_atom(atom_index(i))) lat_list(group_lat_num) = k
end do
end do
end do
do i = 1, group_ele_num
do j = 1, group_sbox_num
if (sbox_list(j) == sbox_ele(element_index(i))) exit
group_sbox_num = group_sbox_num + 1
sbox_list(group_sbox_num) = sbox_ele(element_index(i))
end do
!First calculate max position in lattice space to be able to allocate lat_points array, also sum the total numbers of do j = 1, group_lat_num
!degrees of freedom which are added if (lat_list(group_lat_num) == lat_ele(element_index(i))) exit
dof = 0 group_lat_num = group_lat_num + 1
lat_list(group_lat_num) = lat_ele(element_index(i))
end do
end do
select case(trim(adjustl(shape))) do is = 1, group_sbox_num
case('block')
group_in_lat = reshape((/ block_bd(1),block_bd(3),block_bd(5), & orient = sub_box_ori(:, :, sbox_list(is))
block_bd(2),block_bd(3),block_bd(5), & call matrix_inverse(orient,3,ori_inv)
block_bd(2),block_bd(4),block_bd(5), &
block_bd(1),block_bd(4),block_bd(5), & do ilat = 1, group_lat_num
block_bd(1),block_bd(3),block_bd(6), &
block_bd(2),block_bd(3),block_bd(6), & !First calculate max position in lattice space to be able to allocate lat_points array, also sum the total number
block_bd(2),block_bd(4),block_bd(6), & !of degrees of freedom which are added
block_bd(1),block_bd(4),block_bd(6) /), [3,8]) dof = 0
do j = 1, 3
group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/remesh_lat_pam)) group_bd(2*j) = -huge(1.0_dp)
group_bd(2*j-1) = huge(1.0_dp)
end do
do i = 1, group_atom_num
if ((type_atom(atom_index(i)) == basis_type(1,ilat)).and.(sbox_atom(atom_index(i)) == is)) then
do j =1 ,3
if (r_atom(j,atom_index(i)) > group_bd(2*j)) group_bd(2*j) = r_atom(j,atom_index(i))
if (r_atom(j,atom_index(i)) < group_bd(2*j-1)) group_bd(2*j-1) = r_atom(j,atom_index(i))
end do
dof = dof + 1
end if
end do
do i = 1, group_ele_num
if ((lat_ele(element_index(i)) == ilat).and.(sbox_ele(element_index(i)) == is)) then
do inod =1, ng_node(ilat)
do ibasis = 1, basisnum(ilat)
do j = 1, 3
r =r_node(j,ibasis,inod,element_index(i))
if (r(j) > group_bd(2*j)) group_bd(2*j) = r(j)
if (r(j) < group_bd(2*j-1)) group_bd(2*j-1) = r(j)
end do
end do
end do
dof = dof + size_ele(element_index(i))**3
end if
end do
!If for some reason there are no dof in this loop then cycle out
if(dof == 0) cycle
group_in_lat = reshape((/ group_bd(1),group_bd(3),group_bd(5), &
group_bd(2),group_bd(3),group_bd(5), &
group_bd(2),group_bd(4),group_bd(5), &
group_bd(1),group_bd(4),group_bd(5), &
group_bd(1),group_bd(3),group_bd(6), &
group_bd(2),group_bd(3),group_bd(6), &
group_bd(2),group_bd(4),group_bd(6), &
group_bd(1),group_bd(4),group_bd(6) /), [3,8])
group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/lapa(ilat)))
do i = 1, 3 do i = 1, 3
bd_in_lat(2*i-1) = nint(minval(group_in_lat(i,:))) bd_in_lat(2*i-1) = nint(minval(group_in_lat(i,:)))
bd_in_lat(2*i) = nint(maxval(group_in_lat(i,:))) bd_in_lat(2*i) = nint(maxval(group_in_lat(i,:)))
end do end do
end select 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.
dof = 0
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)) !Now place all group atoms and group interpolated atoms into lat_points
lat_points(:,:,:) = .false. do i = 1, group_atom_num
dof = 0 r = r_atom(:,atom_index(i))/lapa(ilat)
r = matmul(fcc_inv,matmul(ori_inv,r))
!Now place all group atoms and group interpolated atoms into lat_points do j = 1, 3
do i = 1, group_atom_num r_lat(j) = nint(r(j))
r = r_atom(:,atom_index(i))/remesh_lat_pam end do
r = matmul(fcc_inv,matmul(ori_inv,r)) !Do a check to make sure the code is working and that lattice points aren't being written on top of each other.
do j = 1, 3 !This is primarily a debugging statement
r_lat(j) = nint(r(j)) if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then
end do stop "Multiple atoms share same position in lat point array, this shouldn't happen"
!Do a check to make sure the code is working and that lattice points aren't being written on top of each other. else
!This is primarily a debugging statement lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true.
if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then dof = dof + 1
stop "Multiple atoms share same position in lat point array, this shouldn't happen" end if
else end do
lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true.
dof = dof + 1
end if
end do
!Now place interpolated atoms within lat_points array !Now place interpolated atoms within lat_points array
do i =1, group_ele_num do i =1, group_ele_num
ie = element_index(i) ie = element_index(i)
call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
do j = 1, size_ele(ie)**3 * basisnum(lat_ele(ie)) do j = 1, size_ele(ie)**3 * basisnum(lat_ele(ie))
r = r_interp(:,j)/remesh_lat_pam r = r_interp(:,j)/lapa(ilat)
r = matmul(fcc_inv,matmul(ori_inv,r)) r = matmul(fcc_inv,matmul(ori_inv,r))
do k = 1, 3 do k = 1, 3
r_lat(k) = nint(r(k)) r_lat(k) = nint(r(k))
end do
!Do a check to make sure the code is working and that lattice points aren't being written on top of each
!other. This is primarily a debugging statement
if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then
stop "Multiple atoms/interpolated atoms share same position in lat point array, this shouldn't happen"
else
lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true.
dof = dof + 1
end if
end do
end do end do
!Do a check to make sure the code is working and that lattice points aren't being written on top of each other.
!This is primarily a debugging statement
if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then
stop "Multiple atoms/interpolated atoms share same position in lat point array, this shouldn't happen"
else
lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true.
dof = dof + 1
end if
end do
end do
print *, "Group has ", dof, " degrees of freedom to remesh" print *, "Group has ", dof, " degrees of freedom to remesh"
!Delete all elements and atoms to make space for new elements and atoms !Delete all elements and atoms to make space for new elements and atoms
call delete_atoms(group_atom_num, atom_index) call delete_atoms(group_atom_num, atom_index)
call delete_elements(group_ele_num, element_index) call delete_elements(group_ele_num, element_index)
old_atom = atom_num old_atom = atom_num
old_ele = ele_num old_ele = ele_num
!Now run remeshing algorithm, not the most optimized or efficient but gets the job done !Now run remeshing algorithm, not the most optimized or efficient but gets the job done
!Figure out new looping boundaries !Figure out new looping boundaries
bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10 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(2) = bd_in_lat(4) - bd_in_lat(3) + 10
bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10 bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10
if (max_remesh) then if (max_remesh) then
max_loops = (remesh_size-3)/2 max_loops = (remesh_size-3)/2
else else
max_loops = 1 max_loops = 1
end if end if
do j = 1, max_loops do j = 1, max_loops
working_esize = remesh_size - 2*(j-1) working_esize = remesh_size - 2*(j-1)
ele = (working_esize-1)*cubic_cell ele = (working_esize-1)*cubic_cell
zloop: do iz = 1, bd_in_array(3) zloop: do iz = 1, bd_in_array(3)
yloop: do iy = 1, bd_in_array(2) yloop: do iy = 1, bd_in_array(2)
xloop: do ix = 1, bd_in_array(1) xloop: do ix = 1, bd_in_array(1)
if (lat_points(ix, iy,iz)) then if (lat_points(ix, iy,iz)) then
r_new_node(:,:,:) = 0.0_dp r_new_node(:,:,:) = 0.0_dp
!Check to see if the element overshoots the bound !Check to see if the element overshoots the bound
if (iz+working_esize-1 > bd_in_array(3)) then if (iz+working_esize-1 > bd_in_array(3)) then
exit zloop exit zloop
else if (iy+working_esize-1 > bd_in_array(2)) then else if (iy+working_esize-1 > bd_in_array(2)) then
cycle zloop cycle zloop
else if (ix+working_esize-1 > bd_in_array(1)) then else if (ix+working_esize-1 > bd_in_array(1)) then
cycle yloop cycle yloop
end if end if
if (all(lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1))) then
do inod = 1, ng_node(ilat)
vlat = ele(:,inod) + (/ix, iy, iz /)
do i = 1, 3
vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5
end do
r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat)
end do
lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1) = .false.
!Add the element, for the sbox we just set it to the same sbox that we get the orientation
!from. In this case it is from the sbox of the first atom in the group.
call add_element(remesh_ele_type, working_esize, ilat, sbox_atom(atom_index(1)),r_new_node)
end if
end if
end do xloop
end do yloop
end do zloop
end do
if (all(lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1))) then !Now we have to add any leftover lattice points as atoms
do inod = 1, ng_node(remesh_type) do iz = 1, bd_in_array(3)
vlat = ele(:,inod) + (/ix, iy, iz /) do iy=1, bd_in_array(2)
do i = 1, 3 do ix = 1, bd_in_array(1)
vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5 if(lat_points(ix,iy,iz)) then
end do vlat = (/ ix, iy, iz /)
r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam do i = 1, 3
vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5
end do end do
lat_points(ix,iy,iz) = .false.
lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1) = .false. r = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat)
!Add the element, for the sbox we just set it to the same sbox that we get the orientation from. call add_atom(basis_type(1,ilat), sbox_atom(atom_index(1)), r)
!In this case it is from the sbox of the first atom in the group.
call add_element(remesh_ele_type, working_esize, remesh_type, sbox_atom(atom_index(1)),r_new_node)
end if end if
end if
end do xloop
end do yloop
end do zloop
end do
!Now we have to add any leftover lattice points as atoms
do iz = 1, bd_in_array(3)
do iy=1, bd_in_array(2)
do ix = 1, bd_in_array(1)
if(lat_points(ix,iy,iz)) then
vlat = (/ ix, iy, iz /)
do i = 1, 3
vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5
end do end do
lat_points(ix,iy,iz) = .false. end do
r = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam
call add_atom(remesh_type, sbox_atom(atom_index(1)), r)
end if
end do end do
end do end do
end do end do
print *, "Remeshing has created ", ele_num-old_ele, " elements and ", atom_num-old_atom, " atoms." print *, "Remeshing has created ", ele_num-old_ele, " elements and ", atom_num-old_atom, " atoms."
end subroutine remesh_group end subroutine remesh_group

Write Preview
Loading…
Cancel
Save