diff --git a/src/elements.f90 b/src/elements.f90 index 21de971..3b66afb 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -30,19 +30,15 @@ module elements !Below are lattice type arrays which provide information on the general form of the elements. !We currently have a limit of 10 lattice types for simplicities sake but this can be easily increased. - + integer :: lattice_types = 0 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), 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 - real(kind=dp) :: box_bd(6) - + integer :: max_basisnum, basisnum(10) !Max basis atom number, number of basis atoms in each lattice type + integer :: basis_type(10,10) public contains @@ -89,7 +85,6 @@ module elements max_basisnum = 0 basisnum(:) = 0 - basis_pos(:,:,:) = 0.0_dp ng_node(:) = 0 end subroutine lattice_init @@ -304,4 +299,5 @@ module elements return end subroutine rhombshape + end module elements \ No newline at end of file diff --git a/src/io.f90 b/src/io.f90 index 7850ee2..4801de1 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -3,6 +3,7 @@ module io use elements use parameters use atoms + use box implicit none @@ -56,18 +57,10 @@ module io cycle end if select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:)) - case('xyz') + case('xyz', 'lmp', 'vtk', 'mb') outfilenum=outfilenum+1 outfiles(outfilenum) = temp_outfile exit - case('lmp') - 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, lmp, vtk." @@ -96,6 +89,8 @@ module io call write_lmp(outfiles(i)) case('vtk') call write_vtk(outfiles(i)) + case('mb') + call write_mb(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" @@ -276,4 +271,53 @@ module io end do close(11) end subroutine + + subroutine write_mb(file) + + !This subroutine writes the cacmb formatted file which provides necessary information for building models + character(len=100), intent(in) :: file + + integer :: i, j, inod, ibasis + + !Open the .mb file for writing + open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') + + !First write the box boundary information + !Write the global box boundaries + write(11,*) box_bd(:) + !Write the number of sub_boxes in the system + write(11,*) sub_box_num + !For every subbox write the orientation and sub box boundary + do i = 1, sub_box_num + write(11,*) sub_box_ori(:,:,i) + write(11,*) sub_box_bd(:,i) + end do + + !Write the number of atom types in the current model and all of their names + write(11,*) atom_types, (type_to_name(i), i=1, atom_types) + !Write the number of lattice_types, basisnum and number of nodes for each lattice type + write(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types) + !Now for every lattice type write the basis atom types + write(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = 1, lattice_types) + + !Now write the numbers of elements and atoms + write(11,*) atom_num, ele_num + + !Write out atoms first + do i = 1, atom_num + write(11,*) type_atom(i), r_atom(:,i) + end do + + !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, *) i, lat_ele(i), size_ele(i), 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) + end do + end do + end do + + end subroutine write_mb end module io \ No newline at end of file diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 52342e2..2d41ae8 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -6,13 +6,15 @@ module mode_create use io use subroutines use elements + use box 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), 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) + integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), & + basis_pos(3,10) logical :: dup_flag, dim_flag real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:) @@ -98,7 +100,7 @@ module mode_create !Add the basis atoms to the unit cell do inod = 1, max_ng_node do ibasis = 1, basisnum(1) - r_node_temp(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis,1) + origin(:) + r_node_temp(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis) + origin(:) end do end do do i = 1,3 @@ -115,7 +117,7 @@ module mode_create if(lat_atom_num > 0) then do i = 1, lat_atom_num do ibasis = 1, basisnum(1) - call add_atom(basis_type(ibasis,1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis,1)) + call add_atom(basis_type(ibasis, 1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis)) end do end do deallocate(r_atom_lat) @@ -125,7 +127,7 @@ module mode_create 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) + r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis) end do end do call add_element(element_type, esize, 1, r_node_temp) @@ -258,13 +260,17 @@ module mode_create !Now normalize the orientation matrix orient = matrix_normal(orient,3) + !Set lattice_num to 1 + lattice_types = 1 + !If we haven't defined a basis then define the basis and add the default basis atom type and position if (basisnum(1) == 0) then max_basisnum = 1 basisnum(1) = 1 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 + basis_pos(:,1) = 0.0_dp end if + ! end subroutine subroutine build_with_rhomb(box_in_lat, transform_matrix) @@ -432,12 +438,11 @@ module mode_create 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 iz = 1, bd_in_array(3) do iy = 1, bd_in_array(2) - do iz = 1, bd_in_array(1) + do ix = 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) /) @@ -453,7 +458,6 @@ module mode_create end do end do - print *, lat_atom_num end if end subroutine build_with_rhomb