diff --git a/src/Makefile b/src/Makefile index c84186e..e1a1145 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,6 +1,7 @@ + FC=ifort -#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -FFLAGS=-mcmodel=large -Ofast +FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin +#FFLAGS=-mcmodel=large -Ofast MODES=mode_create.o mode_merge.o mode_convert.o OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES) @@ -8,7 +9,7 @@ OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box .SUFFIXES: .c .f .f90 .F90 .o cacmb: $(OBJECTS) - $(FC) $(FFLAGS) $(OBJECTS) -o $@ + $(FC) $(FFLAGS) $(OBJECTS) parameters.o -o $@ .f90.o: $(FC) $(FFLAGS) -c $< @@ -30,4 +31,4 @@ main.o io.o build_subroutines.o: elements.o call_mode.o : $(MODES) $(MODES) io.o: atoms.o box.o $(MODES) main.o : io.o -testfuncs.o elements.o mode_create.o: subroutines.o +testfuncs.o elements.o mode_create.o: subroutines.o \ No newline at end of file diff --git a/src/elements.f90 b/src/elements.f90 index 22e6eff..d6545d6 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -101,10 +101,28 @@ module elements real(kind=dp), dimension(3,max_ng_node), intent(out) :: cell_mat + integer :: inod, i + real(kind=dp), dimension(3,max_ng_node) :: adjustVar + + adjustVar(:,:) = 0.0_dp + select case(trim(ele_type)) case('fcc') - cell_mat(:,1:8) = lapa * ((esize-1)*matmul(orient_mat, fcc_cell)) + if(lmpcac) then + do inod = 1, 8 + do i = 1,3 + if(is_equal(cubic_cell(i, inod),0.0_dp)) then + adjustVar(i,inod) = -0.5_dp + else + adjustVar(i, inod) = 0.5_dp + end if + end do + end do + adjustVar(:,1:8) = matmul(fcc_mat, adjustVar(:,1:8)) + end if + cell_mat(:, 1:8) = fcc_cell + adjustVar(:,1:8) + cell_mat(:,1:8) = lapa * ((esize-1)*matmul(orient_mat, cell_mat(:,1:8))) case default print *, "Element type ", trim(ele_type), " currently not accepted" stop diff --git a/src/io.f90 b/src/io.f90 index 6e914b2..9d5537a 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -9,6 +9,7 @@ module io integer :: outfilenum = 0, infilenum = 0 character(len=100) :: outfiles(10), infiles(10) + public contains @@ -62,9 +63,14 @@ module io outfilenum=outfilenum+1 outfiles(outfilenum) = temp_outfile exit + case('cac') + lmpcac = .true. + 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." + 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, cac." read(*,*) temp_outfile end select @@ -72,7 +78,6 @@ module io end subroutine get_out_file - subroutine write_out !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine @@ -94,6 +99,8 @@ module io call write_mb(outfiles(i)) case('restart') call write_pycac(outfiles(i)) + case('cac') + call write_lmpcac(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" @@ -193,6 +200,71 @@ module io end do end subroutine write_lmp + subroutine write_lmpcac(file) + !This subroutine writes out a .lmp style dump file + character(len=100), intent(in) :: file + integer :: write_num, i, inod, ibasis + real(kind=dp) :: mass + +1 format(i16, ' Eight_Node', 4i16) +2 format(i16, ' Atom', 4i16) +3 format(3i16,3f23.15) + + open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') + + !Comment line + write(11, '(a)') '# CAC input file made with cacmb' + write(11, '(a)') + !Calculate total atom number + write_num = atom_num + ele_num + + !Write total number of atoms + elements + write(11, '(i16, a)') write_num, ' cac elements' + !Write number of atom types + write(11, '(i16, a)') atom_types, ' atom types' + + write(11,'(a)') ' ' + !Write box bd + write(11, '(2f23.15, a)') box_bd(1:2), ' xlo xhi' + write(11, '(2f23.15, a)') box_bd(3:4), ' ylo yhi' + write(11, '(2f23.15, a)') box_bd(5:6), ' zlo zhi' + + !Masses + write(11, '(a)') 'Masses' + + write(11, '(a)') ' ' + do i =1, atom_types + call atommass(type_to_name(i),mass) + write(11, '(i16, f23.15, 2a)') i, mass, ' # ', type_to_name(i) + end do + write(11, '(a)') ' ' + + write(11, '(a)') 'CAC Elements' + write(11, '(a)') ' ' + + !Write element nodal positions + do i = 1, ele_num + select case(trim(adjustl(type_ele(i)))) + case('fcc') + !The first entry is the element specifier + write(11,1) i, basisnum(lat_ele(i)), size_ele(i), size_ele(i), size_ele(i) + do ibasis = 1, basisnum(lat_ele(i)) + do inod = 1, 8 + !Nodal information for every node + write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i) + end do + end do + end select + end do + + do i = 1, atom_num + !Element specifier dictating that it is an atom + write(11,2) ele_num+i, 1, 1, 1, 1 + !Write the atomic information + write(11,3) 1, 1, type_atom(i), r_atom(:,i) + end do + end subroutine write_lmpcac + subroutine write_vtk(file) !This subroutine writes out a vtk style dump file integer :: i, j, inod, ibasis @@ -249,7 +321,7 @@ module io 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)) + write(11, '(3f23.15)') sum(r_node(:,:,inod,i),2)/basisnum(lat_ele(i)) end do end do end do @@ -601,4 +673,4 @@ module io !just overwrite the arrays if(in_eles > 0) lattice_types = lattice_types + new_lattice_types end subroutine read_mb -end module io \ No newline at end of file +end module io diff --git a/src/main.f90 b/src/main.f90 index a681d33..f90b731 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -27,6 +27,17 @@ program main ! Command line parsing arg_num = command_argument_count() + !First check if we are writing out to lammpscac format by looping over all arguments + do i = 1, arg_num + call get_command_argument(i, argument) + select case(argument(scan(argument,'.',.true.)+1:)) + case('cac') + lmpcac = .true. + case default + continue + end select + end do + !Determine if a mode is being used and what it is. The first argument has to be the mode !if a mode is being used call get_command_argument(1, argument) diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 23b9100..a669f93 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -12,7 +12,7 @@ module mode_create 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), duplicate(3) + orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3), adjustVar(3,8) integer :: esize, 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 @@ -299,7 +299,7 @@ module mode_create !Internal variables integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,8), rzero(3), & vlat(3), temp_lat(3,8), m, n, o - real(kind=dp) :: v(3), temp_nodes(3,1,8) + real(kind=dp) :: v(3), temp_nodes(3,1,8), adjustVar(3,8) real(kind=dp), allocatable :: resize_lat_array(:,:) logical, allocatable :: lat_points(:,:,:) logical :: node_in_bd(8) @@ -307,6 +307,23 @@ module mode_create !Do some value initialization max_esize = esize + !If we are writing out to lammpscac format then we have to adjust the nodal positions + + if(lmpcac) then + do inod = 1, 8 + do i = 1,3 + if(is_equal(cubic_cell(i, inod),0.0_dp)) then + adjustVar(i,inod) = -0.5_dp + else + adjustVar(i, inod) = 0.5_dp + end if + end do + end do + + adjustVar(:,1:8) = matmul(orient,matmul(fcc_mat,adjustVar(:,1:8))) + else + adjustVar(:,:)=0.0_dp + end if !First find the bounding lattice points (min and max points for the box in each dimension) numlatpoints = 1 do i = 1, 3 @@ -324,7 +341,6 @@ module mode_create continue end select - !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 if (esize==2) then @@ -437,9 +453,11 @@ module mode_create end do if(all(node_in_bd)) then - lat_ele_num = lat_ele_num+1 - r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) - + lat_ele_num = lat_ele_num+1 + do inod = 1, 8 + r_lat(:,inod,lat_ele_num) = temp_nodes(:,1,inod) + adjustVar(:,inod) + end do + !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,:)) @@ -505,4 +523,4 @@ module mode_create end subroutine error_message -end module mode_create \ No newline at end of file +end module mode_create diff --git a/src/parameters.f90 b/src/parameters.f90 index 0443622..9e9d9d8 100644 --- a/src/parameters.f90 +++ b/src/parameters.f90 @@ -5,4 +5,6 @@ module parameters integer, parameter :: dp= selected_real_kind(15,307) real(kind=dp), parameter :: lim_zero = epsilon(1.0_dp), & lim_large = huge(1.0_dp) + logical, save :: lmpcac + end module parameters \ No newline at end of file