diff --git a/src/Makefile b/src/Makefile index c83df06..58507b4 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,24 +1,24 @@ FC=ifort -#FFLAGS=-c -mcmodel=large -debug -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -FFLAGS=-c -mcmodel=large -Ofast +FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone +#FFLAGS=-c -mcmodel=large -Ofast MODES=mode_create.o -OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o build_subroutines.o $(MODES) +OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o $(MODES) .SUFFIXES: .SUFFIXES: .c .f .f90 .F90 .o cacmb: $(OBJECTS) - $(FC) $(OBJECTS) -o $@ + $(FC) $(FFLAGS) $(OBJECTS) -o $@ .f90.o: - $(FC) $(FFLAGS) $< + $(FC) $(FFLAGS) -c $< .PHONY: clean clean: $(RM) cacmb *.o testfuncs: testfuncs.o functions.o subroutines.o - $(FC) testfuncs.o functions.o build_subroutines.o subroutines.o elements.o -o $@ + $(FC) testfuncs.o functions.o subroutines.o elements.o -o $@ .PHONY: cleantest cleantest: @@ -31,4 +31,3 @@ call_mode.o : $(MODES) $(MODES) io.o: atoms.o $(MODES) main.o : io.o testfuncs.o elements.o mode_create.o: subroutines.o -testfuncs.o : build_subroutines.o diff --git a/src/call_mode.f90 b/src/call_mode.f90 index 8fcdb0a..1ad4f0f 100644 --- a/src/call_mode.f90 +++ b/src/call_mode.f90 @@ -12,7 +12,7 @@ subroutine call_mode(arg_num,mode) select case(mode) case('--create') - call create(arg_num) + call create() case default print *, "Mode ", mode, " currently not accepted. Please check documentation for ", & diff --git a/src/elements.f90 b/src/elements.f90 index 53b827d..21de971 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -13,13 +13,17 @@ module elements real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array integer :: ele_num=0 !Number of elements + integer :: node_num=0 !Total number of nodes !Data structure used to represent atoms - integer, allocatable :: tag_atom(:) !atom id - character(len=100), allocatable:: type_atom(:) !atom type + integer, allocatable :: tag_atom(:), type_atom(:)!atom id real(kind =dp),allocatable :: r_atom(:,:) !atom position integer :: atom_num=0 !Number of atoms + !Mapping atom type to provided name + character(len=2), dimension(10) :: type_to_name + 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) integer :: cubic_faces(4,6) @@ -28,12 +32,12 @@ module elements !We currently have a limit of 10 lattice types for simplicities sake but this can be easily increased. integer :: max_ng_node, ng_node(10) !Max number of nodes per element and number of nodes per element for each lattice type + integer :: max_esize=0 !Maximum number of atoms per side of element !These variables contain information on the basis, for simplicities sake we limit !the user to the definition of 10 lattice types with 10 basis atoms at each lattice point. !This can be easily increased with no change to efficiency - integer :: max_basisnum, basisnum(10)!Max basis atom number, number of basis atoms in each lattice type - character(len=2) :: basis_type(10,10) !Atom type of each basis + integer :: max_basisnum, basisnum(10), basis_type(10,10)!Max basis atom number, number of basis atoms in each lattice type real(kind=dp) :: basis_pos(3,10,10) !Basis atom positions !Simulation cell parameters @@ -150,13 +154,14 @@ module elements size_ele(ele_num) = size lat_ele(ele_num) = lat r_node(:,:,:,ele_num) = r(:,:,:) + node_num = node_num + ng_node(lat) end subroutine add_element subroutine add_atom(type, r) !Subroutine which adds an atom to the atom arrays - character(len=2), intent(in) :: type + integer, intent(in) :: type real(kind=dp), intent(in), dimension(3) :: r atom_num = atom_num+1 @@ -166,6 +171,32 @@ module elements end subroutine add_atom + subroutine add_atom_type(type, inttype) + !This subroutine adds a new atom type to the module list of atom types + character(len=2), intent(in) :: type + integer, intent(out) :: inttype + + integer :: i + logical :: exists + + exists = .false. + do i=1, 10 + if(type == type_to_name(i)) exists = .true. + inttype = i + end do + + if (exists.eqv..false.) then + atom_types = atom_types+1 + if(atom_types > 10) then + print *, "Defined atom types are greater than 10 which is currently not supported." + stop 3 + end if + type_to_name(atom_types) = type + inttype = atom_types + end if + return + end subroutine add_atom_type + subroutine def_ng_node(n, element_types) !This subroutine defines the maximum node number among n element types integer, intent(in) :: n !Number of element types @@ -184,4 +215,93 @@ module elements if(ng_node(i) > max_ng_node) max_ng_node = ng_node(i) end do end subroutine def_ng_node + + subroutine set_max_esize + !This subroutine sets the maximum esize + max_esize=maxval(size_ele) + end subroutine + + subroutine interpolate_atoms(type, esize, lat_type, r_in, type_interp, r_interp) + !This subroutine returns the interpolated atoms from the elements. + + !Arguments + character(len=100), 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(3,max_basisnum, max_ng_node), intent(in) :: r_in !Nodal positions + integer, dimension(max_basisnum*max_esize**3), intent(out) :: type_interp !Interpolated atomic positions + real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: r_interp !Interpolated atomic positions + + !Internal variables + integer :: i, it, is, ir, ibasis, inod, ia, bnum, lat_type_temp + real(kind=dp), allocatable :: a_shape(:) + real(kind=dp) :: t, s, r + + !Initialize some variables + r_interp(:,:) = 0.0_dp + type_interp(:) = 0 + 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(trim(adjustl(type))) + case('fcc') + allocate(a_shape(8)) + !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 rhombshape(r,s,t,a_shape) + + do ibasis = 1, bnum + ia = ia + 1 + do inod = 1, 8 + type_interp(ia) = basis_type(ibasis,lat_type_temp) + r_interp(:,ia) = r_interp(:,ia) + a_shape(inod) * r_in(:,ibasis,inod) + + end do + end do + + end do + end do + end do + end select + + if (ia /= esize**3) then + print *, "Incorrect interpolation" + stop 3 + end if + return + end subroutine interpolate_atoms + + subroutine rhombshape(r,s,t, shape_fun) + !Shape function for rhombohedral elements + real(kind=8), intent(in) :: r, s, t + real(kind=8), intent(out) :: shape_fun(8) + + shape_fun(1) = (1.0-r)*(1.0-s)*(1.0-t)/8.0 + shape_fun(2) = (1.0+r)*(1.0-s)*(1.0-t)/8.0 + shape_fun(3) = (1.0+r)*(1.0+s)*(1.0-t)/8.0 + shape_fun(4) = (1.0-r)*(1.0+s)*(1.0-t)/8.0 + shape_fun(5) = (1.0-r)*(1.0-s)*(1.0+t)/8.0 + shape_fun(6) = (1.0+r)*(1.0-s)*(1.0+t)/8.0 + shape_fun(7) = (1.0+r)*(1.0+s)*(1.0+t)/8.0 + shape_fun(8) = (1.0-r)*(1.0+s)*(1.0+t)/8.0 + + + return + end subroutine rhombshape end module elements \ No newline at end of file diff --git a/src/functions.f90 b/src/functions.f90 index e7c7a05..eb8e813 100644 --- a/src/functions.f90 +++ b/src/functions.f90 @@ -225,7 +225,7 @@ END FUNCTION StrDnCase real(kind=dp), intent(in) :: A, B logical :: is_equal - if((A>(B - 10.0_dp**-10)).and.(A < (B+10.0_dp**-10))) then + if((A>(B - 10.0_dp**(-10))).and.(A < (B+10.0_dp**(-10)))) then is_equal = .true. else is_equal = .false. diff --git a/src/io.f90 b/src/io.f90 index 1bfa649..7850ee2 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -64,10 +64,13 @@ module io outfilenum=outfilenum+1 outfiles(outfilenum) = temp_outfile exit - + case('vtk') + outfilenum=outfilenum+1 + outfiles(outfilenum)=temp_outfile + exit case default print *, "File type: ", trim(temp_outfile(scan(temp_outfile,'.',.true.):)), "not currently accepted. ", & - "please input a filename with extension from following list: xyz." + "please input a filename with extension from following list: xyz, lmp, vtk." read(*,*) temp_outfile end select @@ -81,6 +84,9 @@ module io integer :: i + !Find max esize which will be needed later + call set_max_esize + do i = 1, outfilenum !Pull out the extension of the file and call the correct write subroutine select case(trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:)))) @@ -88,6 +94,8 @@ module io call write_xyz(outfiles(i)) case('lmp') call write_lmp(outfiles(i)) + case('vtk') + call write_vtk(outfiles(i)) case default print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), & " is not accepted for writing. Please select from: xyz and try again" @@ -138,10 +146,10 @@ module io end subroutine write_xyz subroutine write_lmp(file) - - integer :: write_num, i - character(len=100), intent(in) :: file !This subroutine writes out a .lmp style dump file + character(len=100), intent(in) :: file + integer :: write_num, i, iatom, type_interp(max_basisnum*max_esize**3) + real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), mass open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') @@ -150,10 +158,13 @@ module io write(11, '(a)') !Calculate total atom number write_num = atom_num + do i = 1,ele_num + if(type_ele(i) == 'fcc') write_num = write_num + size_ele(i)**3 + end do !Write total number of atoms + elements write(11, '(i16, a)') write_num, ' atoms' !Write number of atom types - write(11, '(i16, a)') 1, ' atom types' + write(11, '(i16, a)') atom_types, ' atom types' write(11,'(a)') ' ' !Write box bd @@ -163,15 +174,106 @@ module io !Masses write(11, '(a)') 'Masses' + write(11, '(a)') ' ' - write(11, '(i16, f23.15)') 1, 63.546 + do i =1, atom_types + call atommass(type_to_name(i),mass) + write(11, '(i16, f23.15)') i, mass + end do write(11, '(a)') ' ' !Write atom positions write(11, '(a)') 'Atoms' write(11, '(a)') ' ' do i = 1, atom_num - write(11, '(2i16, 3f23.15)') i, 1, r_atom(:,i) + write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i) + end do + + !Write refined element atomic positions + 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') + do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 + write(11, '(2i16, 3f23.15)') atom_num+iatom, type_interp(iatom), r_interp(:,iatom) + end do + end select end do end subroutine write_lmp + + subroutine write_vtk(file) + !This subroutine writes out a vtk style dump file + integer :: i, j, inod, ibasis + character(len=100), intent(in) :: file + +1 format('# vtk DataFile Version 4.0.1', / & + 'CAC output -- cg', / & + 'ASCII') +11 format('# vtk DataFile Version 4.0.1', / & + 'CACmb output -- atoms', / & + 'ASCII') +2 format('DATASET UNSTRUCTURED_GRID') +3 format('POINTS', i16, ' float') +4 format(/'CELLS', 2i16) +5 format(/'CELL_TYPES', i16) +12 format(/'CELL_DATA', i16) +16 format(/'POINT_DATA', i16) +17 format('SCALARS weight float', / & + 'LOOKUP_TABLE default') +18 format('SCALARS atom_type float', / & + 'LOOKUP_TABLE default') + +20 format('SCALARS lattice_type float', /& + 'LOOKUP_TABLE default') + + !First we write the vtk file containing the atoms + open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind') + + write(11, 11) + write(11, 2) + write(11, 3) atom_num + do i = 1, atom_num + write(11, '(3f23.15)') r_atom(:,i) + end do + write(11,4) atom_num, atom_num*2 + do i = 1, atom_num + write(11, '(2i16)') 1, i-1 + end do + write(11, 5) atom_num + do i = 1, atom_num + write(11, '(i16)') 1 + end do + write(11, 16) atom_num + write(11, 18) + do i = 1, atom_num + write(11, '(i16)') type_atom(i) + end do + close(11) + + open(unit=11, file='cg_'//trim(adjustl(file)), action='write', status='replace',position='rewind') + write(11,1) + write(11,2) + write(11,3) node_num + do i = 1, ele_num + do inod=1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + write(11, '(3f23.1)') sum(r_node(:,:,inod,i),2)/basisnum(lat_ele(i)) + end do + end do + end do + write(11, 4) ele_num, ele_num + node_num + do i =1, ele_num + write(11, '(9i16)') ng_node(lat_ele(i)), (j, j = (i-1)*ng_node(lat_ele(i)), i*ng_node(lat_ele(i))-1) + end do + write(11,5) ele_num + do i = 1, ele_num + if(trim(adjustl(type_ele(i))) == 'fcc') write(11, '(i16)') 12 + end do + write(11,12) ele_num + write(11,20) + do i = 1, ele_num + write(11, '(i16)') lat_ele(i) + end do + close(11) + end subroutine end module io \ No newline at end of file diff --git a/src/mode_create.f90 b/src/mode_create.f90 index e6bb0be..52342e2 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -5,28 +5,27 @@ module mode_create use atoms use io use subroutines + use elements implicit none 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), ox_bd(6), maxbd(3), lattice_space(3) - integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_num, lat_atom_num, bd_in_lat(6) + orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3) + integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6) logical :: dup_flag, dim_flag real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:) public contains - subroutine create(arg_num) + subroutine create() ! Main subroutine which controls execution - integer, intent(in) :: arg_num character(len=100) :: textholder integer :: i, ibasis, inod - real(kind=dp) :: r(3), periodvone(3), periodvtwo(3) - real(kind=dp), allocatable :: r_node(:,:,:) + real(kind=dp), allocatable :: r_node_temp(:,:,:) !Initialize default parameters orient = reshape((/ 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp /), shape(orient)) @@ -40,11 +39,11 @@ module mode_create dup_flag = .false. dim_flag = .false. basisnum = 0 - lat_num = 0 + lat_ele_num = 0 lat_atom_num = 0 !First we parse the command - call parse_command(arg_num) + call parse_command() ! Before building do a check on the file if (outfilenum == 0) then @@ -55,7 +54,7 @@ module mode_create !Now we setup the unit element and call other init subroutines call def_ng_node(1, element_type) - allocate(r_node(3,max_basisnum,max_ng_node)) + allocate(r_node_temp(3,max_basisnum,max_ng_node)) if(dup_flag) then @@ -64,7 +63,7 @@ module mode_create do i = 1, 8 - box_vert(:,i) = duplicate(:)*lattice_space(:)*cubic_cell(:,i) + origin(:) + box_vert(:,i) = duplicate(:)*esize*lattice_space(:)*cubic_cell(:,i) + origin(:) end do call matrix_inverse(orient,3,orient_inv) !Now get the rotated box vertex positions in lattice space. Should be integer units @@ -72,31 +71,14 @@ module mode_create !Find the new maxlen maxbd = maxval(matmul(orient,matmul(fcc_mat,box_lat_vert)),2) do i = 1, 3 - box_bd(2*i) = maxval(box_vert(i,:)) - 0.1_dp*lattice_space(i) - box_bd(2*i-1) = origin(i) + box_bd(2*i) = maxval(box_vert(i,:)) - 0.25_dp*lattice_space(i) + box_bd(2*i-1) = origin(i)-0.25_dp*lattice_space(i) end do - !and then call the build function with the correct transformation matrix select case(trim(adjustl(element_type))) case('fcc') - ! periodvone(:) = matmul(orient, reshape((/ 1, 1, 0 /),(/ 3 /))) - ! periodvtwo(:) = matmul(orient, reshape((/ 1, 1, 2 /),(/ 3 /))) - - ! do i = 1, 3 - ! if (periodvone(i) < lim_zero) then - ! ! box_bd(2*i) =floor(box_bd(2*i)/periodvtwo(i))*periodvtwo(i) - ! box_bd(2*i) = box_bd(2*i) - 0.5*periodvtwo(i) - ! else if(periodvtwo(i) < lim_zero) then - ! ! box_bd(2*i) =floor(box_bd(2*i)/periodvone(i))*periodvone(i) - ! box_bd(2*i) = box_bd(2*i) - 0.5*periodvone(i) - ! else - ! ! box_bd(2*i) = floor(box_bd(2*i)/lcm(periodvone(i),periodvtwo(i)))*lcm(periodvone(i),periodvtwo(i)) - ! box_bd(2*i) = box_bd(2*i) - 0.5*lcm(periodvone(i),periodvtwo(i)) - - ! end if - ! end do - - call lattice_in_box(box_lat_vert, fcc_mat) + + call build_with_rhomb(box_lat_vert, fcc_mat) case default print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", & "element type" @@ -116,33 +98,45 @@ module mode_create !Add the basis atoms to the unit cell do inod = 1, max_ng_node do ibasis = 1, basisnum(1) - r_node(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis,1) + origin(:) + r_node_temp(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis,1) + origin(:) end do end do - - call add_element(element_type, esize, 1, r_node) + do i = 1,3 + box_bd(2*i) = maxval(r_node_temp(i,:,:)) + box_bd(2*i-1) = origin(i) + end do + call add_element(element_type, esize, 1, r_node_temp) end if !If we passed the dup_flag or dim_flag then we have to convert the lattice points and add them to the atom/element arrays if(dup_flag.or.dim_flag) then !Allocate variables - call alloc_ele_arrays(lat_num, lat_atom_num*basisnum(1)) + call alloc_ele_arrays(lat_ele_num, lat_atom_num*basisnum(1)) if(lat_atom_num > 0) then - !Check for periodicity do i = 1, lat_atom_num - do ibasis = 1, basisnum(1) + do ibasis = 1, basisnum(1) call add_atom(basis_type(ibasis,1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis,1)) end do end do deallocate(r_atom_lat) end if + + if(lat_ele_num > 0) then + do i = 1, lat_ele_num + do inod= 1, ng_node(1) + do ibasis = 1, basisnum(1) + r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis,1) + end do + end do + call add_element(element_type, esize, 1, r_node_temp) + end do + end if end if end subroutine create !This subroutine parses the command and pulls out information needed for mode_create - subroutine parse_command(arg_num) + subroutine parse_command() - integer, intent(in) :: arg_num integer :: arg_pos, ori_pos, i, j, arglen, stat character(len=100) :: textholder @@ -217,7 +211,6 @@ module mode_create read(textholder, *) origin(i) arg_pos = arg_pos + 1 end do - print *, origin !If a filetype is passed then we add name.ext to the outfiles list case('xyz') textholder = trim(adjustl(name)) //'.xyz' @@ -244,13 +237,11 @@ module mode_create select case(trim(adjustl(element_type))) case('fcc') do i = 1,3 - print *, orient(i,:) !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.& (is_equal(abs(orient(i,2)), abs(orient(i,3))).and.(is_equal(orient(i,1),0.0_dp))).or.& (is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(orient(i,2),0.0_dp)))) then - print *, '110', i lattice_space(i) = 0.5_dp * lattice_space(i) !Check if one direction is 112 @@ -259,7 +250,6 @@ module mode_create (is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(abs(orient(i,2)),2.0_dp*abs(orient(i,3))))))& then - print *, '112 ', i lattice_space(i) = 0.5_dp * lattice_space(i) end if @@ -272,21 +262,27 @@ module mode_create if (basisnum(1) == 0) then max_basisnum = 1 basisnum(1) = 1 - basis_type(1,1) = name !If basis command not defined then we use name as the atom_name + call add_atom_type(name, basis_type(1,1)) !If basis command not defined then we use name as the atom_name basis_pos(:,1,1) = 0.0_dp end if end subroutine - subroutine lattice_in_box(box_in_lat, transform_matrix) + subroutine build_with_rhomb(box_in_lat, transform_matrix) !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 !Internal variables - integer :: i, j, bd_in_lat(6), ix, iy, iz, numlatpoints - real(kind=dp) :: box_face_center(3,6), face_normals(3,6), Cx(2), Cy, Cz, A(2), v(3) + integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, templatpoints, ele(3,8), rzero(3), ilat, & + type_interp(basisnum(1)*esize**3), vlat(3), temp_lat(3,8), m, n, o + real(kind=dp) :: v(3), temp_nodes(3,1,8), ele_atoms(3,esize**3), r_interp(3,basisnum(1)*esize**3) real(kind=dp), allocatable :: resize_lat_array(:,:) + logical, allocatable :: lat_points(:,:,:) + logical :: node_in_bd(8) + + !Do some value initialization + max_esize = esize !First find the bounding lattice points (min and max points for the box in each dimension) numlatpoints = 1 @@ -304,82 +300,10 @@ module mode_create case default continue end select - !Calculate the box_face centroids and box face normals. This is used in the centroid code. - box_face_center(:,:) = 0.0_dp - face_normals = reshape((/ -1.0_dp, 0.0_dp, 0.0_dp, & - 1.0_dp, 0.0_dp, 0.0_dp, & - 0.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, & - 0.0_dp, 0.0_dp, 1.0_dp /),& - shape(face_normals)) - !Face normals - select case(trim(adjustl(element_type))) - case('fcc') - do i = 1, 6 - !Box face normal - face_normals(:,i) = matmul(fcc_inv, matmul(orient_inv, face_normals(:,i))) - end do - end select - !Face centroids - do i =1, 6 - - !Initialize variables - A(:) = 0.0_dp - Cx(:) = 0.0_dp - Cy = 0.0_dp - Cz = 0.0_dp - - !Calculate all terms using a projection onto the xy and xz planes and then using the 2d equation - !for centroid of a plane. This is why we calculate the x centroid twice. - do j = 1, 4 - ! A(1) = A(1) + 0.5*(box_in_lat(1,cubic_faces(j,i))*box_in_lat(2,cubic_faces(j+1,i)) & - ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(2,cubic_faces(j,i))) - - ! !Centroid in x from xy plane - ! Cx(1) = Cx(1) + (box_in_lat(1,cubic_faces(j,i))+box_in_lat(1,cubic_faces(j+1,i)))* & - ! (box_in_lat(1,cubic_faces(j,i))*box_in_lat(2,cubic_faces(j+1,i)) & - ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(2,cubic_faces(j,i))) - - ! !Centroid in y from xy plane - ! Cy = Cy + (box_in_lat(2,cubic_faces(j,i))+box_in_lat(2,cubic_faces(j+1,i)))* & - ! (box_in_lat(1,cubic_faces(j,i))*box_in_lat(2,cubic_faces(j+1,i)) & - ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(2,cubic_faces(j,i))) - - ! A(2) = A(2) + 0.5*(box_in_lat(1,cubic_faces(j,i))*box_in_lat(3,cubic_faces(j+1,i)) & - ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(3,cubic_faces(j,i))) - - ! !Centroid in x from xz plane - ! Cx(2) = Cx(2) + (box_in_lat(1,cubic_faces(j,i))+box_in_lat(1,cubic_faces(j+1,i)))* & - ! (box_in_lat(1,cubic_faces(j,i))*box_in_lat(3,cubic_faces(j+1,i)) & - ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(3,cubic_faces(j,i))) - - ! !Centroid in z from xz plane - ! Cz = Cz + (box_in_lat(3,cubic_faces(j,i))+box_in_lat(3,cubic_faces(j+1,i)))* & - ! (box_in_lat(1,cubic_faces(j,i))*box_in_lat(3,cubic_faces(j+1,i)) & - ! - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(3,cubic_faces(j,i))) - - ! print *, j, i, Cx, Cy, Cz, A - Cx(1) = Cx(1) + box_in_lat(1,cubic_faces(j,i))*0.25 - Cy = Cy + box_in_lat(2,cubic_faces(j,i))*0.25 - Cz = Cz + box_in_lat(3,cubic_faces(j,i))*0.25 - - end do - - ! Cx = Cx * 1/(6*A) - ! if(Cx(1) /= Cx(2)) then - ! call error_message(7) - ! end if - ! Cy = Cy* 1/(6*A(1)) - ! Cz = Cz*1/(6*A(2)) - - box_face_center(:,i) = (/ Cx(1), Cy, Cz /) - end do !Loop over all of lattice points within the boundary, we choose between two loops. One for the atomistic case !and one for the regular case - print *, box_bd if (esize==2) then !atomistics do iz = bd_in_lat(5)-5, bd_in_lat(6)+5 @@ -398,10 +322,141 @@ module mode_create end do end do else - continue + !If we are working with elements we have to use more complex code + + !Initialize finite element + ele(:,:) = (esize-1) * cubic_cell(:,:) + + !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 + 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. + do iz = bd_in_lat(5)-5, bd_in_lat(6)+5 + do iy = bd_in_lat(3)-5, bd_in_lat(4)+5 + do ix = bd_in_lat(1)-5, bd_in_lat(2)+5 + v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /) + + !Transform point back to real space for easier checking + v = matmul(orient, matmul(transform_matrix,v)) + !If within the boundaries + if(in_block_bd(v, box_bd)) then + lat_points(ix-bd_in_lat(1)+5,iy-bd_in_lat(3)+5,iz-bd_in_lat(5) + 5) = .true. + end if + end do + end do + end do + + !Now we redefine bd_in_lat The first 3 indices contains limits for the lat_points array + 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 + !Figure out where the starting point is. This is the first piont which fully contains the finite element + outerloop: do iz = 1, bd_in_array(3) + do iy = 1, bd_in_array(2) + do ix = 1, bd_in_array(1) + node_in_bd(8) = .false. + do inod = 1, 8 + vlat = ele(:,inod) + (/ ix, iy, iz /) + + !Check to see if the node resides at a position containing a lattice point within the box + if(any(vlat > shape(lat_points))) then + continue + else if(lat_points(vlat(1),vlat(2),vlat(3))) then + node_in_bd(inod) = .true. + end if + end do + + if(all(node_in_bd)) then + rzero = (/ ix, iy, iz /) + exit outerloop + end if + end do + end do + end do outerloop + + !Now build the finite element region + lat_ele_num = 0 + lat_atom_num = 0 + allocate(r_lat(3,8,numlatpoints/esize)) + + !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 + ele(:,inod) = ele(:,inod) + rzero + end do + do iz = -bd_in_array(6), bd_in_array(6) + do iy = -bd_in_array(5), bd_in_array(5) + do ix = -bd_in_array(4), bd_in_array(4) + node_in_bd(:) = .false. + temp_nodes(:,:,:) = 0.0_dp + temp_lat(:,:) = 0 + do inod = 1, 8 + 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)) + do i = 1,3 + v(i) = real(vlat(i) + bd_in_lat(2*i-1) - 5) + end do + temp_nodes(:,1, inod) = matmul(orient, matmul(transform_matrix, v)) + temp_lat(:,inod) = vlat + + !Check to see if the lattice point values are greater then the array limits + if(any(vlat > shape(lat_points)).or.any(vlat < 1)) then + continue + !If within array boundaries check to see if it is a lattice point + else if(lat_points(vlat(1),vlat(2),vlat(3))) then + node_in_bd(inod) = .true. + end if + end do + + if(all(node_in_bd)) then + lat_ele_num = lat_ele_num+1 + r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) + + !Now set all the lattice points contained within an element to false + do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:)) + do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:)) + do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:)) + lat_points(m,n,o) = .false. + end do + end do + end do + + end if + end do + end do + end do + + !Now figure out how many lattice points could not be contained in elements + print *, count(lat_points) + allocate(r_atom_lat(3,count(lat_points))) + lat_atom_num = 0 + do ix = 1, bd_in_array(3) + do iy = 1, bd_in_array(2) + do iz = 1, bd_in_array(1) + !If this point is a lattice point then save the lattice point as an atom + if (lat_points(ix,iy,iz)) then + v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /) + do i = 1,3 + v(i) = v(i) + real(bd_in_lat(2*i-1) - 5) + end do + !Transform point back to real space + v = matmul(orient, matmul(transform_matrix,v)) + lat_atom_num = lat_atom_num + 1 + r_atom_lat(:,lat_atom_num) = v + end if + end do + end do + end do + + print *, lat_atom_num end if - end subroutine lattice_in_box + end subroutine build_with_rhomb subroutine error_message(errorid)