diff --git a/src/Makefile b/src/Makefile index 4705b1e..d26a8e0 100644 --- a/src/Makefile +++ b/src/Makefile @@ -2,8 +2,8 @@ .DEFAULT_GOAL := all FC=mpif90 -FFLAGS=-Wall -mcmodel=large -O0 -g -fbacktrace -fcheck=all -ffpe-trap=invalid,zero,overflow,underflow,denormal -FFLAGS=-O2 -mcmodel=large +#FFLAGS=-Wall -mcmodel=large -O0 -g -fbacktrace -fcheck=all -ffpe-trap=invalid,zero,overflow,underflow,denormal +FFLAGS=-O3 -g OBJDIR=obj SRCS := $(wildcard *.f90) diff --git a/src/elements.f90 b/src/elements.f90 index 3b3fc28..cb64943 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -33,7 +33,8 @@ module elements integer :: atom_types = 0 !Variables for creating elements based on primitive cells - real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3), bcc_cell(3,8), bcc_mat(3,3), bcc_inv(3,3) + real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3), bcc_cell(3,8), bcc_mat(3,3), bcc_inv(3,3), & + cube20(3,20) integer :: cubic_faces(4,6) !Below are lattice type arrays which provide information on the general form of the elements. @@ -72,6 +73,29 @@ module elements 0.0_dp, 1.0_dp, 1.0_dp /), & shape(fcc_cell)) + !Now initialize the 20 node cube element + cube20 = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, & + 1.0_dp, 0.0_dp, 0.0_dp, & + 1.0_dp, 1.0_dp, 0.0_dp, & + 0.0_dp, 1.0_dp, 0.0_dp, & + 0.0_dp, 0.0_dp, 1.0_dp, & + 1.0_dp, 0.0_dp, 1.0_dp, & + 1.0_dp, 1.0_dp, 1.0_dp, & + 0.0_dp, 1.0_dp, 1.0_dp, & + 0.5_dp, 0.0_dp, 0.0_dp, & + 1.0_dp, 0.5_dp, 0.0_dp, & + 0.5_dp, 1.0_dp, 0.0_dp, & + 0.0_dp, 0.5_dp, 0.0_dp, & + 0.0_dp, 0.0_dp, 0.5_dp, & + 1.0_dp, 0.0_dp, 0.5_dp, & + 1.0_dp, 1.0_dp, 0.5_dp, & + 0.0_dp, 1.0_dp, 0.5_dp, & + 0.5_dp, 0.0_dp, 1.0_dp, & + 1.0_dp, 0.5_dp, 1.0_dp, & + 0.5_dp, 1.0_dp, 1.0_dp, & + 0.0_dp, 0.5_dp, 1.0_dp /), & + shape(cube20)) + !Now we create a list containing the list of vertices needed to describe the 6 cube faces cubic_faces(:,1) = (/ 1, 2, 3, 4 /) cubic_faces(:,2) = (/ 1, 2, 6, 5 /) @@ -384,6 +408,8 @@ module elements ng_node(i) = 8 case('bcc') ng_node(i) = 8 + case('20fcc') + ng_node(i) = 20 end select if(ng_node(i) > max_ng_node) max_ng_node = ng_node(i) @@ -462,6 +488,34 @@ module elements end do end do + end do + end do + end do + case('20fcc') + !Now loop over all the possible sites + do it = 1, esize + t = (1.0_dp*(it-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2) + do is =1, esize + s = (1.0_dp*(is-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2) + do ir = 1, esize + r = (1.0_dp*(ir-1) - (esize-1)/2)/(1.0_dp*(esize-1)/2) + call quad_rhombshape(r,s,t,a_shape) + + do ibasis = 1, bnum + ia = ia + 1 + do inod = 1, 20 + type_interp(ia) = basis_type(ibasis,lat_type_temp) + r_interp(:,ia) = r_interp(:,ia) + a_shape(inod) * r_in(:,ibasis,inod) + + if(present(data_interp)) then + !If data is present then interpolate data arrays as well + data_interp(1,ia) = data_interp(1,ia) + eng(ibasis, inod)*a_shape(inod) + data_interp(2:4,ia) = data_interp(2:4,ia) + f(:, ibasis, inod)*a_shape(inod) + data_interp(5:10,ia) = data_interp(5:10,ia) + v(:, ibasis, inod)*a_shape(inod) + end if + end do + end do + end do end do end do @@ -492,6 +546,42 @@ module elements return end subroutine rhombshape + subroutine quad_rhombshape(r,s,t, shape_fun) + !Shape function for 20 node quadratic rhombohedral elements + real(kind=8), intent(in) :: r, s, t + real(kind=8), intent(out) :: shape_fun(:) + + !Corner nodes + shape_fun(1) = (1.0_dp-r)*(1.0_dp-s)*(1.0_dp-t)*(-r-s-t-2)/8.0_dp + shape_fun(2) = (1.0_dp+r)*(1.0_dp-s)*(1.0_dp-t)*(r-s-t-2)/8.0_dp + shape_fun(3) = (1.0_dp+r)*(1.0_dp+s)*(1.0_dp-t)*(r+s-t-2)/8.0_dp + shape_fun(4) = (1.0_dp-r)*(1.0_dp+s)*(1.0_dp-t)*(-r+s-t-2)/8.0_dp + shape_fun(5) = (1.0_dp-r)*(1.0_dp-s)*(1.0_dp+t)*(-r-s+t-2)/8.0_dp + shape_fun(6) = (1.0_dp+r)*(1.0_dp-s)*(1.0_dp+t)*(r-s+t-2)/8.0_dp + shape_fun(7) = (1.0_dp+r)*(1.0_dp+s)*(1.0_dp+t)*(r+s+t-2)/8.0_dp + shape_fun(8) = (1.0_dp-r)*(1.0_dp+s)*(1.0_dp+t)*(-r+s+t-2)/8.0_dp + + !Side nodes, first node r is zero + shape_fun(9) = (1-r*r)*(1-s)*(1-t)/4.0_dp + shape_fun(11) = (1-r*r)*(1+s)*(1-t)/4.0_dp + shape_fun(17) = (1-r*r)*(1-s)*(1+t)/4.0_dp + shape_fun(19) = (1-r*r)*(1+s)*(1+t)/4.0_dp + + !node s is zero + shape_fun(10) = (1+r)*(1-s*s)*(1-t)/4.0_dp + shape_fun(12) = (1-r)*(1-s*s)*(1-t)/4.0_dp + shape_fun(18) = (1+r)*(1-s*s)*(1+t)/4.0_dp + shape_fun(20) = (1-r)*(1-s*s)*(1+t)/4.0_dp + + !node t is zero + shape_fun(13) = (1-r)*(1-s)*(1-t*t)/4.0_dp + shape_fun(14) = (1+r)*(1-s)*(1-t*t)/4.0_dp + shape_fun(15) = (1+r)*(1+s)*(1-t*t)/4.0_dp + shape_fun(16) = (1-r)*(1+s)*(1-t*t)/4.0_dp + + + end subroutine quad_rhombshape + subroutine delete_atoms(num, in_index) !This subroutine deletes atoms from the atom arrays integer, intent(in) :: num @@ -624,6 +714,9 @@ do i = 1, atom_num else if(trim(adjustl(pos_string)) == '-inf') then pos=box_bd(2*i-1) + else if(trim(adjustl(pos_string)) == 'mid') then + pos=(box_bd(2*i)+box_bd(2*i-1))/2 + else if (trim(adjustl(pos_string)) == 'rand') then call random_number(rand) pos = (box_bd(2*i)-box_bd(2*i-1))*rand + box_bd(2*i-1) @@ -647,7 +740,6 @@ do i = 1, atom_num end do end if - read(cone, *) face if ((face < 1).or.(face > 6)) stop "Current face number must be 1,2,3,4,5,6. Please check documentation" !Now get the position call offset_pos(ele, face, rand_ele_pos) @@ -670,6 +762,15 @@ do i = 1, atom_num end if pos = box_bd(2*i) - pos + else if ((index(pos_string,'-') > 0).and.(index(pos_string,'mid')>0)) then + !Now extract the number we are reducing from midinity + if(index(pos_string,'mid') < index(pos_string,'-')) then + read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos + else + read(pos_string(1:index(pos_string,'-')-1), *, iostat=iospara) pos + end if + pos = (box_bd(2*i)+box_bd(2*i-1))/2.0_dp - pos + else if ((index(pos_string,'+') > 0).and.(index(pos_string,'inf')>0)) then !Now extract the number we are reducing from infinity if(index(pos_string,'inf') < index(pos_string,'+')) then @@ -679,6 +780,15 @@ do i = 1, atom_num end if pos = box_bd(2*i-1) + pos + else if ((index(pos_string,'+') > 0).and.(index(pos_string,'mid')>0)) then + !Now extract the number we are reducing from midinity + if(index(pos_string,'mid') < index(pos_string,'+')) then + read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos + else + read(pos_string(1:index(pos_string,'+')-1), *, iostat=iospara) pos + end if + pos = (box_bd(2*i)+box_bd(2*i-1))/2.0_dp + pos + else if ((index(pos_string,'*') > 0).and.(index(pos_string,'inf')>0)) then !Now extract the number we are reducing from infinity if(index(pos_string,'inf') < index(pos_string,'*')) then @@ -809,6 +919,7 @@ do i = 1, atom_num integer, intent(in) :: n,m !n-size of element arrays, m-size of atom arrays integer :: allostat + print *, max_ng_node !Allocate element arrays if (n > 0) then if (allocated(force_node)) then diff --git a/src/io.f90 b/src/io.f90 index e78fef7..62a2687 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -177,6 +177,7 @@ module io do i = 1,ele_num if(type_ele(i) == 'fcc') write_num = write_num + size_ele(i)**3 if(type_ele(i) == 'bcc') write_num = write_num + size_ele(i)**3 + if(type_ele(i) == '20fcc') write_num = write_num + size_ele(i)**3 end do !Write total number of atoms + elements @@ -212,7 +213,7 @@ module io do i = 1, ele_num call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp) select case(trim(adjustl(type_ele(i)))) - case('fcc','bcc') + case('fcc','bcc', '20fcc') do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 interp_num = interp_num+1 call apply_periodic(r_interp(:,iatom)) @@ -255,6 +256,7 @@ module io do i = 1,ele_num if(type_ele(i) == 'fcc') write_num = write_num + size_ele(i)**3 if(type_ele(i) == 'bcc') write_num = write_num + size_ele(i)**3 + if(type_ele(i) == '20fcc') write_num = write_num + size_ele(i)**3 end do end if !Write total number of atoms @@ -311,7 +313,7 @@ module io call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp) end if select case(trim(adjustl(type_ele(i)))) - case('fcc','bcc') + case('fcc','bcc', '20fcc') if(write_dat) then do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 interp_num = interp_num+1 @@ -562,7 +564,12 @@ module io ip = 0 write(11,13) do i = 1, ele_num - write(11, '(4i16)') i, basisnum(lat_ele(i)), 2, (size_ele(i)-1) + select case(type_ele(i)) + case('fcc','bcc') + write(11, '(4i16)') i, basisnum(lat_ele(i)), 2, (size_ele(i)-1) + case('20fcc') + write(11, '(4i16)') i, basisnum(lat_ele(i)), 3, (size_ele(i)-1) + end select do inod = 1, ng_node(lat_ele(i)) do ibasis = 1, basisnum(lat_ele(i)) ip = ip + 1 @@ -624,7 +631,7 @@ module io !Write out the elements, this is written in two stages, one line for the element and then 1 line for !every basis at every node do i = 1, ele_num - write(11, *) tag_ele(i), lat_ele(i), size_ele(i), sbox_ele(i), type_ele(i) + write(11, *) tag_ele(i), lat_ele(i), size_ele(i), sbox_ele(i), trim(adjustl(type_ele(i))) do inod = 1, ng_node(lat_ele(i)) do ibasis =1, basisnum(lat_ele(i)) write(11,*) inod, ibasis, r_node(:, ibasis, inod, i) @@ -1003,7 +1010,7 @@ module io 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) real(kind=dp) :: newdisplace(3), ra(3), in_lapa, ea, fa(3), va(6), & - ee(10,8), fe(3,10,8), ve(6,10,8), re(3,10,8), atomic_masses(10) + ee(10,20), fe(3,10,20), ve(6,10,20), re(3,10,20), atomic_masses(10) character(len=100) :: textholder, fcc character(len=1000) :: line character(len=2) :: atomic_element @@ -1013,7 +1020,7 @@ module io !Now initialize some important variables if they aren't defined if (max_basisnum==0) max_basisnum = 10 - if (max_ng_node==0) max_ng_node=8 + if (max_ng_node==0) max_ng_node=20 fcc="fcc" !Skip header comment lines @@ -1096,6 +1103,11 @@ module io inbtypes(ibasis) = atom_type_map(inbtypes(ibasis)) end do call lattice_map(bnum, inbtypes, n, 1.0_dp, lat_type) + if (n == 8)then + fcc="fcc" + else if (n == 20) then + fcc = "20fcc" + end if call add_element(tag, fcc, esize+1, lat_type, sub_box_num, re) call add_element_data(ele_num, ee, fe, ve) end do @@ -1135,6 +1147,8 @@ module io read(11,*) temp_box_bd(3:4), textholder, textholder read(11,*) temp_box_bd(5:6), textholder, textholder + call grow_ele_arrays(0,atom_in) + print *, "Read in ", atom_in, " atoms from ", trim(adjustl(file)) !Shift the box boundaries if needed do i = 1, 3 diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 29e4d05..12bb2b8 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -67,7 +67,7 @@ module mode_create end do !Now get the rotated box vertex positions in lattice space. Should be integer units and get the new maxlen select case(trim(adjustl(element_type))) - case('fcc') + case('fcc', '20fcc') box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1 maxbd = maxval(matmul(orient,matmul(fcc_mat,box_lat_vert)),2) case('bcc') @@ -93,7 +93,7 @@ module mode_create end do !Now get the rotated box vertex positions in lattice space. Should be integer units select case(trim(adjustl(element_type))) - case('fcc') + case('fcc', '20fcc') box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1 case('bcc') box_lat_vert = int(matmul(bcc_inv, matmul(orient_inv, box_vert)))+1 @@ -131,9 +131,11 @@ module mode_create !Call the build function with the correct transformation matrix select case(trim(adjustl(element_type))) case('fcc') - call build_with_rhomb(box_lat_vert, fcc_mat) + call build_with_rhomb(box_lat_vert, fcc_mat, 8) case('bcc') - call build_with_rhomb(box_lat_vert, bcc_mat) + call build_with_rhomb(box_lat_vert, bcc_mat, 8) + case('20fcc') + call build_with_rhomb(box_lat_vert, fcc_mat, 20) case default print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ",& "element type" @@ -280,7 +282,7 @@ module mode_create !Check special periodicity relations select case(trim(adjustl(element_type))) - case('fcc') + case('fcc', '20fcc') do i = 1,3 !Check if one of the directions is 110 if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(orient(i,3),0.0_dp))).or.& @@ -331,19 +333,20 @@ module mode_create ! end subroutine - subroutine build_with_rhomb(box_in_lat, transform_matrix) + subroutine build_with_rhomb(box_in_lat, transform_matrix, nn) !This subroutine returns all the lattice points in the box in r_lat !Inputs integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space + integer, intent(in) :: nn !Internal variables - integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,8), rzero(3), efill_size, & - vlat(3), temp_lat(3,8), m, n, o, j, k, nump_ele, efill_temp_lat(3,8), filzero(3), bd_ele_lat(6),& - efill_ele(3,8), ebd(6) - real(kind=dp) :: v(3), temp_nodes(3,1,8), r(3), centroid_bd(6) + integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,nn), rzero(3), efill_size, & + vlat(3), temp_lat(3,nn), m, n, o, j, k, nump_ele, efill_temp_lat(3,nn), filzero(3), & + bd_ele_lat(6), efill_ele(3,nn), ebd(6) + real(kind=dp) :: v(3), temp_nodes(3,1,nn), r(3), centroid_bd(6) logical, allocatable :: lat_points(:,:,:) - logical :: node_in_bd(8), add, lat_points_ele(esize,esize,esize), intersect_bd(3) + logical :: node_in_bd(nn), add, lat_points_ele(esize,esize,esize), intersect_bd(3) !Do some value initialization max_esize = esize @@ -393,7 +396,12 @@ module mode_create !If we are working with elements we have to use more complex code !Initialize finite element - ele(:,:) = (esize-1) * cubic_cell(:,:) + select case(element_type) + case('fcc','bcc') + ele(:,:) = (esize-1) * cubic_cell(:,:) + case('20fcc') + ele(:,:) = (esize-1)*cube20 + end select !Make a 3 dimensional array of lattice points. This array is indexed by the integer lattice position. !A value of true means that the coordinate is a lattice point which is within the box_bd @@ -423,7 +431,7 @@ module mode_create do iy = 1, bd_in_array(2) do ix = 1, bd_in_array(1) node_in_bd(:) = .false. - do inod = 1, 8 + do inod = 1, nn vlat = ele(:,inod) + (/ ix, iy, iz /) !Check to see if the node resides at a position containing a lattice point within the box @@ -445,16 +453,22 @@ module mode_create !Now build the finite element region lat_ele_num = 0 lat_atom_num = 0 - allocate(r_lat(3,8,numlatpoints/esize), elat(numlatpoints/esize)) + allocate(r_lat(3,nn,numlatpoints/esize), elat(numlatpoints/esize)) - ele(:,:) = (esize-1) * cubic_cell(:,:) + select case(element_type) + case('fcc','bcc') + ele(:,:) = (esize-1) * cubic_cell(:,:) + case('20fcc') + ele(:,:) = (esize-1)*cube20 + end select + !Redefined the second 3 indices as the number of elements that fit in the bounds do i = 1, 3 bd_in_array(3+i) = bd_in_array(i)/esize end do !Now start the element at rzero - do inod=1, 8 + do inod=1, nn ele(:,inod) = ele(:,inod) + rzero end do do iz = -bd_in_array(6), bd_in_array(6) @@ -463,7 +477,7 @@ module mode_create node_in_bd(:) = .false. temp_nodes(:,:,:) = 0.0_dp temp_lat(:,:) = 0 - do inod = 1, 8 + do inod = 1, nn vlat= ele(:,inod) + (/ ix*(esize), iy*(esize), iz*(esize) /) !Transform point back to real space for easier checking ! v = matmul(orient, matmul(transform_matrix,v)) @@ -537,7 +551,7 @@ module mode_create ze: do k = 1, (esize-efill_size) ye: do j = 1, (esize-efill_size) xe: do i = 1, (esize-efill_size) - do inod = 1,8 + do inod = 1,nn vlat = efill_ele(:,inod) + (/ i, j, k /) if (.not.lat_points_ele(vlat(1),vlat(2),vlat(3))) cycle xe do o = 1,3 diff --git a/src/mode_metric.f90 b/src/mode_metric.f90 index 79f1ee9..8d0b6a5 100644 --- a/src/mode_metric.f90 +++ b/src/mode_metric.f90 @@ -193,7 +193,7 @@ module mode_metric integer :: i, inod, ibasis npoints=0 - + rout = -huge(1.0_dp) if(atom_num > 0) then do i = 1, atom_num rout(:,tag_atom(i)) = r_atom(:,i) diff --git a/src/neighbors.f90 b/src/neighbors.f90 index d3683e6..27cd608 100644 --- a/src/neighbors.f90 +++ b/src/neighbors.f90 @@ -41,7 +41,6 @@ module neighbors box_len(i) = box_bd(2*i) - box_bd(2*i-1) cell_num(i) = int(box_len(i)/(rc_off/2))+1 end do - print *, box_len, cell_num !Initialize/allocate variables cell_lim = 10 @@ -49,13 +48,15 @@ module neighbors !Now place points within cell num_in_cell = 0 + do i = 1, numinlist !Check to see if the current point is a filler point and if so just skip it - if(r_list(1,i) < -huge(1.0_dp)+1) cycle + if(r_list(1,i) < box_bd(1)) cycle !c is the position of the cell that the point belongs to + if(i == 11651203) print *, r_list(:,i), (r_list(1,i) < box_bd(1)) do j = 1, 3 - c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off/2)) + 1 + c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off)) + 1 end do !Place the index in the correct position, growing if necessary @@ -102,7 +103,7 @@ module neighbors pointloop: do i = 1, n !First check to see if the point is a filler point, if so then skip it - if(r_list(1,i) < -Huge(-1.0_dp)+1) cycle + if(r_list(1,i) < box_bd(1)) cycle !c is the position of the cell that the point c = which_cell(:,i) diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 78f8029..c9236da 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -644,7 +644,7 @@ module opt_group if(group_ele_num > 0) then orig_atom_num = atom_num !Estimate number of atoms we are adding, this doesn't have to be exact - add_atom_num = group_ele_num*basisnum(lat_ele(element_index(1)))*size_ele(element_index(1))**3 + add_atom_num = group_ele_num*max_basisnum*(max_esize)**3 call grow_ele_arrays(0,add_atom_num) do i = 1, group_ele_num ie = element_index(i)