diff --git a/src/io.f90 b/src/io.f90 index 853a230..f929e29 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -907,39 +907,46 @@ module io !Read more useless data read(11,*) textholder - read(11,*) textholder !set max values and allocate variables max_basisnum = maxval(basisnum) max_ng_node = maxval(ng_node) call grow_ele_arrays(in_eles, in_atoms) - + !Now start reading the elements - do i = 1, in_eles - read(11,*) j, etype, k, lat_type - 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 + if(in_eles > 0) then + read(11,*) textholder + do i = 1, in_eles + read(11,*) j, etype, k, lat_type + 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 if - !Read useless data - read(11,*) textholder - read(11,*) textholder + if(in_atoms > 0) then - 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) + if (in_eles > 0) then + !Read useless data + read(11,*) textholder + read(11,*) textholder + end if - 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 module io diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 80bdc26..73ea62c 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -8,9 +8,9 @@ module opt_group use box 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 - 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 integer, allocatable :: element_index(:), atom_index(:) @@ -160,14 +160,6 @@ module opt_group call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) stop "Missing remesh element size in group command" 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') max_remesh =.true. 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, & 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), & - 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(:,:,:) character(len=100) :: remesh_ele_type @@ -322,156 +314,218 @@ module opt_group !Right now we just hardcode only remeshing to elements remesh_ele_type = 'fcc' - !Get the orientations, this assumes that the orientation of the subbox for the first atom is the - !same as the rest of the atoms - !If this assumption is false then the code will break and exit - orient = sub_box_ori(:, :, sbox_atom(atom_index(1))) - call matrix_inverse(orient,3,ori_inv) + ! Determine which sub_boxes and lattices types are within in the group + group_sbox_num = 0 + group_lat_num = 0 + do i = 1, group_atom_num + 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 - !First calculate max position in lattice space to be able to allocate lat_points array, also sum the total numbers of - !degrees of freedom which are added - dof = 0 - - select case(trim(adjustl(shape))) - case('block') - group_in_lat = reshape((/ block_bd(1),block_bd(3),block_bd(5), & - block_bd(2),block_bd(3),block_bd(5), & - block_bd(2),block_bd(4),block_bd(5), & - block_bd(1),block_bd(4),block_bd(5), & - block_bd(1),block_bd(3),block_bd(6), & - block_bd(2),block_bd(3),block_bd(6), & - block_bd(2),block_bd(4),block_bd(6), & - block_bd(1),block_bd(4),block_bd(6) /), [3,8]) - - group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/remesh_lat_pam)) - do i = 1, 3 - bd_in_lat(2*i-1) = nint(minval(group_in_lat(i,:))) - bd_in_lat(2*i) = nint(maxval(group_in_lat(i,:))) + 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 - end select + 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 + + do j = 1, group_lat_num + if (lat_list(group_lat_num) == lat_ele(element_index(i))) exit + group_lat_num = group_lat_num + 1 + lat_list(group_lat_num) = lat_ele(element_index(i)) + end do + end do - 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 + do is = 1, group_sbox_num - !Now place all group atoms and group interpolated atoms into lat_points - do i = 1, group_atom_num - r = r_atom(:,atom_index(i))/remesh_lat_pam - r = matmul(fcc_inv,matmul(ori_inv,r)) - do j = 1, 3 - r_lat(j) = nint(r(j)) - 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 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 + orient = sub_box_ori(:, :, sbox_list(is)) + call matrix_inverse(orient,3,ori_inv) + + do ilat = 1, group_lat_num - !Now place interpolated atoms within lat_points array - do i =1, group_ele_num - ie = element_index(i) - 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)) - r = r_interp(:,j)/remesh_lat_pam - r = matmul(fcc_inv,matmul(ori_inv,r)) - do k = 1, 3 - r_lat(k) = nint(r(k)) + !First calculate max position in lattice space to be able to allocate lat_points array, also sum the total number + !of degrees of freedom which are added + dof = 0 + do j = 1, 3 + 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 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 + + 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 - print *, "Group has ", dof, " degrees of freedom to remesh" + !If for some reason there are no dof in this loop then cycle out + if(dof == 0) cycle - !Delete all elements and atoms to make space for new elements and atoms - call delete_atoms(group_atom_num, atom_index) - call delete_elements(group_ele_num, element_index) + 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]) - old_atom = atom_num - old_ele = ele_num + group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/lapa(ilat))) + do i = 1, 3 + bd_in_lat(2*i-1) = nint(minval(group_in_lat(i,:))) + bd_in_lat(2*i) = nint(maxval(group_in_lat(i,:))) + end do + 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 + + !Now place all group atoms and group interpolated atoms into lat_points + do i = 1, group_atom_num + r = r_atom(:,atom_index(i))/lapa(ilat) + r = matmul(fcc_inv,matmul(ori_inv,r)) + do j = 1, 3 + r_lat(j) = nint(r(j)) + 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 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 + !Now place interpolated atoms within lat_points array + do i =1, group_ele_num + ie = element_index(i) + 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)) + r = r_interp(:,j)/lapa(ilat) + r = matmul(fcc_inv,matmul(ori_inv,r)) + do k = 1, 3 + 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 - !Now run remeshing algorithm, not the most optimized or efficient but gets the job done - !Figure out new looping boundaries - 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(3) = bd_in_lat(6) - bd_in_lat(5) + 10 + print *, "Group has ", dof, " degrees of freedom to remesh" + !Delete all elements and atoms to make space for new elements and atoms + call delete_atoms(group_atom_num, atom_index) + call delete_elements(group_ele_num, element_index) - if (max_remesh) then - max_loops = (remesh_size-3)/2 - else - max_loops = 1 - end if - do j = 1, max_loops - working_esize = remesh_size - 2*(j-1) - ele = (working_esize-1)*cubic_cell - zloop: do iz = 1, bd_in_array(3) - yloop: do iy = 1, bd_in_array(2) - xloop: do ix = 1, bd_in_array(1) - if (lat_points(ix, iy,iz)) then - r_new_node(:,:,:) = 0.0_dp - - !Check to see if the element overshoots the bound - if (iz+working_esize-1 > bd_in_array(3)) then - exit zloop - else if (iy+working_esize-1 > bd_in_array(2)) then - cycle zloop - else if (ix+working_esize-1 > bd_in_array(1)) then - cycle yloop - 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(remesh_type) - 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))*remesh_lat_pam - end do + old_atom = atom_num + old_ele = ele_num - 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, remesh_type, sbox_atom(atom_index(1)),r_new_node) - 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 - lat_points(ix,iy,iz) = .false. - r = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam - call add_atom(remesh_type, sbox_atom(atom_index(1)), r) - end if + !Now run remeshing algorithm, not the most optimized or efficient but gets the job done + !Figure out new looping boundaries + 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(3) = bd_in_lat(6) - bd_in_lat(5) + 10 + + + if (max_remesh) then + max_loops = (remesh_size-3)/2 + else + max_loops = 1 + end if + do j = 1, max_loops + working_esize = remesh_size - 2*(j-1) + ele = (working_esize-1)*cubic_cell + zloop: do iz = 1, bd_in_array(3) + yloop: do iy = 1, bd_in_array(2) + xloop: do ix = 1, bd_in_array(1) + if (lat_points(ix, iy,iz)) then + r_new_node(:,:,:) = 0.0_dp + + !Check to see if the element overshoots the bound + if (iz+working_esize-1 > bd_in_array(3)) then + exit zloop + else if (iy+working_esize-1 > bd_in_array(2)) then + cycle zloop + else if (ix+working_esize-1 > bd_in_array(1)) then + cycle yloop + 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 + + !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 + lat_points(ix,iy,iz) = .false. + r = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat) + call add_atom(basis_type(1,ilat), 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." end subroutine remesh_group