From 54aa50b605f5e8735f944cbf71e9108505417bee Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Fri, 29 Nov 2019 13:36:19 -0500 Subject: [PATCH 1/5] Working writing to lammpsCAC format --- src/io.f90 | 161 +++++++++++++++++++++++++++++++++++++++----- src/mode_create.f90 | 2 - 2 files changed, 146 insertions(+), 17 deletions(-) diff --git a/src/io.f90 b/src/io.f90 index 7850ee2..948fd52 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -6,8 +6,8 @@ module io implicit none - integer :: outfilenum = 0 - character(len=100) :: outfiles(10) + integer :: outfilenum = 0, infilenum = 0 + character(len=100) :: outfiles(10), infiles(10) public contains @@ -41,6 +41,7 @@ module io if((scan(overwrite, "n") > 0).or.(scan(overwrite, "N") > 0)) then print *, "Please specify a new filename with extension:" read(*,*) temp_outfile + cycle else if((scan(overwrite, "y") > 0).or.(scan(overwrite, "Y") > 0)) then continue else @@ -56,21 +57,13 @@ module io cycle end if select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:)) - case('xyz') + case('xyz','lmp','vtk','cac') 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." + 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 @@ -78,6 +71,52 @@ module io end subroutine get_out_file + subroutine get_in_file(filename) + + implicit none + + character(len=100), intent(in) :: filename + character(len=100) :: temp_infile + logical :: file_exists + + !If no filename is provided then this function is called with none and prompts user input + if (filename=='none') then + print *, "Please specify a filename with extension to read in:" + read(*,*) temp_infile + else + temp_infile = filename + end if + + !Infinite loop which only exists if user provides valid filetype + do while(.true.) + + !Check to see if file exists, if it doesn't then ask the user for another input + inquire(file=trim(temp_infile), exist=file_exists) + if (.not.file_exists) then + print *, "The file ", temp_infile, "does not exist, please input an existing file to read in." + read(*,*) temp_infile + cycle + end if + + if (scan(temp_outfile,'.',.true.) == 0) then + print *, "No extension included on filename, please type a full filename that includes an extension." + read(*,*) temp_infile + cycle + end if + select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:)) + case('cac') + infilenum=infilenum+1 + infiles(infilenum) = temp_infile + exit + case default + print *, "File type: ", trim(temp_infile(scan(temp_outfile,'.',.true.):)), " not currently accepted. ", & + "please input a filename with extension from following list: cac." + read(*,*) temp_infile + + end select + end do + + end subroutine get_in_file subroutine write_out !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine @@ -96,6 +135,8 @@ module io call write_lmp(outfiles(i)) case('vtk') call write_vtk(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" @@ -201,6 +242,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 @@ -257,7 +363,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 @@ -276,4 +382,29 @@ module io end do close(11) end subroutine -end module io \ No newline at end of file + + +! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! READ SUBROUTINES !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! subroutine read_lmpcac(file, box_bd) +! !This subroutine reads in a lmpcac file which can be used with different options and modes + +! !Arguments +! character(len=100), intent(in) :: file +! real(kind=wp), dimension(6), intent(out) :: box_bd + +! !Internal variables +! character(len=1000) :: line +! integer :: read_num, atom_lim, ele_lim + +! !Open the lmpcac file +! open(unit=11, file=file, action='read', position='rewind') + +! !Skip header lines +! read(11,*) line +! read(11,*) line + +! !Read total number of elements + + +! end subroutine read_lmpcac +end module io diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 52342e2..e822798 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -432,7 +432,6 @@ 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) @@ -453,7 +452,6 @@ module mode_create end do end do - print *, lat_atom_num end if end subroutine build_with_rhomb From c1150ae4b13b1418d9419452b455d41d028310fb Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Fri, 29 Nov 2019 22:23:22 -0500 Subject: [PATCH 2/5] First attempts at figuring out the element format for lammpscac --- src/io.f90 | 14 ++++++++++---- src/mode_create.f90 | 46 ++++++++++++++++++++++++++++++++++++++------- 2 files changed, 49 insertions(+), 11 deletions(-) diff --git a/src/io.f90 b/src/io.f90 index 948fd52..af66da3 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -8,6 +8,7 @@ module io integer :: outfilenum = 0, infilenum = 0 character(len=100) :: outfiles(10), infiles(10) + logical lmpcac public contains @@ -57,7 +58,12 @@ module io cycle end if select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:)) - case('xyz','lmp','vtk','cac') + case('xyz','lmp','vtk') + outfilenum=outfilenum+1 + outfiles(outfilenum) = temp_outfile + exit + case('cac') + lmpcac = .true. outfilenum=outfilenum+1 outfiles(outfilenum) = temp_outfile exit @@ -98,18 +104,18 @@ module io cycle end if - if (scan(temp_outfile,'.',.true.) == 0) then + if (scan(temp_infile,'.',.true.) == 0) then print *, "No extension included on filename, please type a full filename that includes an extension." read(*,*) temp_infile cycle end if - select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:)) + select case(temp_infile(scan(temp_infile,'.',.true.)+1:)) case('cac') infilenum=infilenum+1 infiles(infilenum) = temp_infile exit case default - print *, "File type: ", trim(temp_infile(scan(temp_outfile,'.',.true.):)), " not currently accepted. ", & + print *, "File type: ", trim(temp_infile(scan(temp_infile,'.',.true.):)), " not currently accepted. ", & "please input a filename with extension from following list: cac." read(*,*) temp_infile diff --git a/src/mode_create.f90 b/src/mode_create.f90 index e822798..53880c9 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -11,7 +11,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) + orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), adjustVar(3,8) 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 @@ -91,14 +91,30 @@ module mode_create continue else - call cell_init(lattice_parameter, esize, element_type, orient, cell_mat) + 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*esize + else + adjustVar(i, inod) = 0.5_dp*esize + end if + end do + end do + else + adjustVar(:,:)=0.0_dp + end if + + ! call cell_init(lattice_parameter, esize, element_type, orient, cell_mat) !If the user doesn't pass any build instructions than we just put the cell mat into the element_array call alloc_ele_arrays(1,0) !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) = lattice_parameter*matmul(orient, & + matmul(fcc_mat, esize*(cubic_cell(:,inod)+adjustVar(:,inod)))) & + + basis_pos(:,ibasis,1) end do end do do i = 1,3 @@ -276,7 +292,7 @@ module mode_create !Internal variables 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) :: v(3), temp_nodes(3,1,8), ele_atoms(3,esize**3), r_interp(3,basisnum(1)*esize**3), adjustVar(3,8) real(kind=dp), allocatable :: resize_lat_array(:,:) logical, allocatable :: lat_points(:,:,:) logical :: node_in_bd(8) @@ -284,6 +300,21 @@ 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.25_dp + else + adjustVar(i, inod) = 0.25_dp + end if + end do + end do + 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 @@ -301,7 +332,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 @@ -415,8 +445,10 @@ module mode_create if(all(node_in_bd)) then lat_ele_num = lat_ele_num+1 - r_lat(:,:,lat_ele_num) = temp_nodes(:,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,:)) From 591762d2591ee26b79313f1ea661685b0b7edf4c Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Wed, 4 Dec 2019 21:03:13 -0500 Subject: [PATCH 3/5] Latest attempts to get the lammpscac data format working --- src/mode_create.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 53880c9..6c465a3 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -95,9 +95,9 @@ module mode_create 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*esize + adjustVar(i,inod) = -0.5_dp else - adjustVar(i, inod) = 0.5_dp*esize + adjustVar(i, inod) = 0.5_dp end if end do end do @@ -113,7 +113,7 @@ module mode_create do inod = 1, max_ng_node do ibasis = 1, basisnum(1) r_node_temp(:,ibasis,inod) = lattice_parameter*matmul(orient, & - matmul(fcc_mat, esize*(cubic_cell(:,inod)+adjustVar(:,inod)))) & + matmul(fcc_mat, (esize+1)*cubic_cell(:,inod)+adjustVar(:,inod))) & + basis_pos(:,ibasis,1) end do end do From ba28734fd7ccf34b2a816910528c57d225cabcc4 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Mon, 13 Jan 2020 16:54:52 -0500 Subject: [PATCH 4/5] Working changes to write lammps cac --- src/mode_create.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 6c465a3..b65d62e 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -306,9 +306,9 @@ module mode_create do inod = 1, 8 do i = 1,3 if(is_equal(cubic_cell(i, inod),0.0_dp)) then - adjustVar(i,inod) = -0.25_dp + adjustVar(i,inod) = -0.5_dp else - adjustVar(i, inod) = 0.25_dp + adjustVar(i, inod) = 0.5_dp end if end do end do From fbae3a15946d52fa55513ee2759ba18389913bc4 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Mon, 13 Jan 2020 20:12:46 -0500 Subject: [PATCH 5/5] Working version of CACmb capable of building simple models for lammpsCAC --- src/Makefile | 9 +++++---- src/elements.f90 | 20 +++++++++++++++++- src/io.f90 | 49 +-------------------------------------------- src/main.f90 | 11 ++++++++++ src/mode_create.f90 | 18 +++-------------- src/parameters.f90 | 2 ++ 6 files changed, 41 insertions(+), 68 deletions(-) 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 dd49096..9d5537a 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -9,7 +9,7 @@ module io integer :: outfilenum = 0, infilenum = 0 character(len=100) :: outfiles(10), infiles(10) - logical lmpcac + public contains @@ -78,53 +78,6 @@ module io end subroutine get_out_file - subroutine get_in_file(filename) - - implicit none - - character(len=100), intent(in) :: filename - character(len=100) :: temp_infile - logical :: file_exists - - !If no filename is provided then this function is called with none and prompts user input - if (filename=='none') then - print *, "Please specify a filename with extension to read in:" - read(*,*) temp_infile - else - temp_infile = filename - end if - - !Infinite loop which only exists if user provides valid filetype - do while(.true.) - - !Check to see if file exists, if it doesn't then ask the user for another input - inquire(file=trim(temp_infile), exist=file_exists) - if (.not.file_exists) then - print *, "The file ", temp_infile, "does not exist, please input an existing file to read in." - read(*,*) temp_infile - cycle - end if - - if (scan(temp_infile,'.',.true.) == 0) then - print *, "No extension included on filename, please type a full filename that includes an extension." - read(*,*) temp_infile - cycle - end if - select case(temp_infile(scan(temp_infile,'.',.true.)+1:)) - case('cac') - infilenum=infilenum+1 - infiles(infilenum) = temp_infile - exit - case default - print *, "File type: ", trim(temp_infile(scan(temp_infile,'.',.true.):)), " not currently accepted. ", & - "please input a filename with extension from following list: cac." - read(*,*) temp_infile - - end select - end do - - end subroutine get_in_file - subroutine write_out !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine 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 d5995e8..a669f93 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -91,20 +91,6 @@ module mode_create end do else - 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 - else - adjustVar(:,:)=0.0_dp - end if - call cell_init(lattice_parameter, esize, element_type, orient, cell_mat) !If the user doesn't pass any build instructions than we just put the cell mat into the element_array call alloc_ele_arrays(1,0) @@ -333,6 +319,8 @@ module mode_create end if end do end do + + adjustVar(:,1:8) = matmul(orient,matmul(fcc_mat,adjustVar(:,1:8))) else adjustVar(:,:)=0.0_dp end if @@ -465,7 +453,7 @@ module mode_create end do if(all(node_in_bd)) then - lat_ele_num = lat_ele_num+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 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