diff --git a/src/Makefile.dep b/src/Makefile.dep index 55f461d..54ee834 100644 --- a/src/Makefile.dep +++ b/src/Makefile.dep @@ -9,7 +9,6 @@ obj/main : \ obj/io.o \ obj/main.o \ obj/mode_calc.o \ - obj/mode_check.o \ obj/mode_convert.o \ obj/mode_create.o \ obj/mode_merge.o \ @@ -24,6 +23,7 @@ obj/main : \ obj/opt_slip_plane.o \ obj/parameters.o \ obj/sorts.o \ + obj/str.o \ obj/subroutines.o obj/atoms.o : \ @@ -37,7 +37,6 @@ obj/box.o : \ obj/caller.o : \ obj/box.o \ obj/mode_calc.o \ - obj/mode_check.o \ obj/mode_convert.o \ obj/mode_create.o \ obj/mode_merge.o \ @@ -65,7 +64,8 @@ obj/io.o : \ obj/atoms.o \ obj/box.o \ obj/elements.o \ - obj/parameters.o + obj/parameters.o \ + obj/str.o obj/main.o : \ obj/caller.o \ @@ -80,16 +80,6 @@ obj/mode_calc.o : \ obj/parameters.o \ obj/subroutines.o -obj/mode_check.o : \ - obj/atoms.o \ - obj/box.o \ - obj/elements.o \ - obj/functions.o \ - obj/io.o \ - obj/neighbors.o \ - obj/parameters.o \ - obj/subroutines.o - obj/mode_convert.o : \ obj/box.o \ obj/elements.o \ @@ -169,6 +159,8 @@ obj/parameters.o : obj/sorts.o : \ obj/parameters.o +obj/str.o : + obj/subroutines.o : \ obj/box.o \ obj/functions.o \ diff --git a/src/elements.f90 b/src/elements.f90 index fffa9cb..6d71e0a 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -14,10 +14,11 @@ module elements integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:), tag_ele(:) !Element size real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array !Element result data structures - real(kind=8), allocatable :: force_node(:,:,:,:), virial_node(:,:,:,:), energy_node(:,:,:) + real(kind=dp), allocatable :: force_node(:,:,:,:), virial_node(:,:,:,:), energy_node(:,:,:) integer, save :: ele_num !Number of elements integer, save :: node_num !Total number of nodes + integer, save :: node_atoms !Count of all basis atoms at nodes summed over all nodes !Data structure used to represent atoms integer, allocatable :: type_atom(:)!atom type @@ -120,6 +121,7 @@ module elements basisnum(:) = 0 ng_node(:) = 0 node_num = 0 + node_atoms = 0 ele_num = 0 atom_num = 0 end subroutine lattice_init @@ -288,6 +290,7 @@ module elements ele_num = ele_num + 1 node_num = node_num + ng_node(lat) + node_atoms = node_atoms + ng_node(lat)*basisnum(lat) if (tag==0) then newtag = ele_num !If we don't assign a tag then pass the tag as the ele_num @@ -386,7 +389,7 @@ module elements end if end subroutine - subroutine interpolate_atoms(type, esize, lat_type, r_in, type_interp, r_interp) + subroutine interpolate_atoms(type, esize, lat_type, r_in, type_interp, r_interp, eng, f, v, data_interp) !This subroutine returns the interpolated atoms from the elements. !Arguments @@ -396,6 +399,9 @@ module elements 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 + real(kind = dp), optional, intent(in) :: eng(max_basisnum, max_ng_node), f(3, max_basisnum, max_ng_node), & + v(6, max_basisnum, max_ng_node) + real(kind=dp), dimension(10, max_basisnum*max_esize**3), optional,intent(out) :: data_interp !Interpolated atomic positions !Internal variables integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp @@ -405,6 +411,7 @@ module elements !Initialize some variables r_interp(:,:) = 0.0_dp type_interp(:) = 0 + if(present(data_interp)) data_interp = 0.0_dp ia = 0 !Define bnum based on the input lattice type. If lat_type=0 then we are interpolating lattice points which means @@ -437,6 +444,12 @@ module elements 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 diff --git a/src/io.f90 b/src/io.f90 index f168a4f..6ee5e64 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -4,6 +4,7 @@ module io use parameters use atoms use box + use str implicit none @@ -59,7 +60,7 @@ module io cycle end if select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:)) - case('xyz', 'lmp', 'vtk', 'mb', 'restart') + case('xyz', 'lmp', 'vtk', 'mb', 'restart', 'dump') outfilenum=outfilenum+1 outfiles(outfilenum) = temp_outfile exit @@ -104,6 +105,8 @@ module io call write_pycac(outfiles(i)) case('cac') call write_lmpcac(outfiles(i)) + case('dump') + call write_ldump(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, lmp, vtk, mb, restart, cac and try again" @@ -124,7 +127,7 @@ module io open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') !Write total number of atoms + elements - write(11, '(i16)') node_num+atom_num + write(11, '(i16)') node_atoms+atom_num !Write comment line write(11, '(a)') "#Node + atom file created using cacmb" @@ -140,8 +143,8 @@ module io end do end do - if(outn /= node_num) then - print *, "outn", outn, " doesn't equal node_num ", node_num + if(outn /= node_atoms) then + print *, "outn", outn, " doesn't equal node_atoms ", node_atoms end if !Write atom positions @@ -150,8 +153,8 @@ module io outn = outn + 1 end do - if((outn-node_num) /= atom_num) then - print *, "outn", (outn-node_num), " doesn't equal atom_num ", atom_num + if((outn-node_atoms) /= atom_num) then + print *, "outn", (outn-node_atoms), " doesn't equal atom_num ", atom_num end if !Finish writing @@ -219,6 +222,87 @@ module io end do end subroutine write_lmp + subroutine write_ldump(file) + !This subroutine will only work if element data is defined + character(len = *), intent(in) :: file + integer :: write_num, i, iatom + logical :: write_dat + integer :: type_interp(max_basisnum*max_esize**3), interp_num + real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), data_interp(10, max_basisnum*max_esize**3) + + + +1 format('ITEM: TIMESTEP'/i16) +2 format('ITEM: NUMBER OF ATOMS' /i16) +3 format('ITEM: BOX BOUNDS ', 2a1, ' ', 2a1, ' ', 2a1 / & + 2f23.15 / 2f23.15 / 2f23.15) +4 format('ITEM: ATOMS id type x y z energy fx fy fz s11 s22 s33 s23 s13 s12') +5 format('ITEM: ATOMS id type x y z') + + open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') + + !Write header information + write(11,1) timestep + + !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 + if(type_ele(i) == 'bcc') write_num = write_num + size_ele(i)**3 + end do + !Write total number of atoms + write(11,2) write_num + !Write box information + write(11,3) box_bc(1:1), box_bc(1:1), box_bc(2:2), box_bc(2:2), box_bc(3:3), box_bc(3:3), box_bd(:) + + !Now pick if we are interpolating data or not + if(allocated(force_node).or.allocated(force_atom)) then + write(11, 4) + write_dat = .true. + else + write(11, 5) + write_dat = .false. + end if + + if (write_dat) then + do i = 1, atom_num + write(11, '(2i16, 13f23.15)') i, type_atom(i), r_atom(:,i), energy_atom(i), force_atom(:,i), virial_atom(:,i) + end do + else + do i = 1, atom_num + write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i) + end do + end if + + !Write refined element atomic positions + interp_num = 0 + do i = 1, ele_num + if(write_dat) then + call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp, & + energy_node(:,:,i), force_node(:,:,:,i), virial_node(:,:,:,i), data_interp) + else + 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') + if(write_dat) then + do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 + interp_num = interp_num+1 + call apply_periodic(r_interp(:,iatom)) + write(11, '(2i16, 13f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom), & + data_interp(:,iatom) + end do + else + do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 + interp_num = interp_num+1 + call apply_periodic(r_interp(:,iatom)) + write(11, '(2i16, 3f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom) + end do + end if + end select + end do + end subroutine write_ldump + subroutine write_lmpcac(file) !This subroutine writes out a .lmp style dump file character(len=100), intent(in) :: file @@ -392,54 +476,30 @@ module io character(len=100), intent(in) :: file integer :: interp_max, i, j, inod, ibasis, ip, unique_index(50), unique_size(50), unique_type(50), unique_num, & etype - real(kind=dp) :: box_vec(3) + real(kind=dp) :: box_vec(3), masses(10) 1 format('time' / i16, f23.15) 2 format('number of elements' / i16) 3 format('number of nodes' / i16) -4 format('element types' / i16) 5 format('number of atoms' / i16) -6 format('number of grains' / i16) 7 format('boundary ' / 3a1) 8 format('box bound' / 6f23.15) 9 format('box length' / 3f23.15) 10 format('box matrix') 11 format(3f23.15) 12 format('coarse-grained domain') -13 format('ie ele_type grain_ele lat_type_ele'/ 'ip ibasis x y z') -14 format('atomistic domain' / 'ia grain_atom type_atom x y z') -15 format('maximum lattice periodicity length' / 3f23.15) -16 format('Number of lattice types and atom types '/ 2i16) -17 format('lattice type IDs') -18 format('lattice types for grains') -19 format('max nodes per element' / i16) +13 format('ie basis_num ng_node esize'/ 'ip ibasis type x y z') +14 format('atomistic domain' / 'ia type_atom x y z') +19 format('max nodes per element and basis per nodes' / 2i16) 20 format('max interpo per element' / i16) 21 format('atom types to elements') open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') + !Below writes the header information for the restart file write(11,1) timestep, total_time write(11,2) ele_num - !Below writes the header information for the restart file - - - !First figure out all of the unique element types - unique_num = 0 - unique_index(:) = 0 - eleloop:do i = 1, ele_num - do j =1 , unique_num - if ( ( size_ele(i) == size_ele( unique_index(j) ) ).and. & - ( lat_ele(i) == lat_ele(unique_index(j)) ) ) then - cycle eleloop - end if - end do - unique_num = unique_num + 1 - unique_index(unique_num) = i - unique_size(unique_num) = size_ele(i) - unique_type(unique_num) = lat_ele(i) - end do eleloop - !Calculate the max number of atoms per element select case(max_ng_node) case(8) @@ -447,31 +507,19 @@ module io case default interp_max = 0 end select + write(11,20) interp_max + !Write write(11,3) node_num - write(11,19) max_ng_node - write(11,4) unique_num + write(11,19) max_ng_node, max_basisnum write(11,5) atom_num - write(11,6) 1 !Grain_num is ignored - write(11,16) lattice_types, atom_types - write(11,21) do i = 1, atom_types - write(11,*) i, type_to_name(i) + call atommass(type_to_name(i),masses(i)) end do + write(11,*) "masses: " + write(11, *) (masses(i), i = 1, atom_types) write(11,7) box_bc(1:1), box_bc(2:2), box_bc(3:3) - write(11,18) - write(11,'(2i16)') 1,1 !This is another throwaway line that is meaningless - write(11,17) - !This may have to be updated in the future but currently the only 8 node element is fcc - do i = 1, lattice_types - select case(ng_node(i)) - case(8) - write(11, *) i, 'fcc' - end select - end do - write(11,15) 1.0_dp, 1.0_dp, 1.0_dp !Another throwaway line that isn't needed write(11,8) box_bd - write(11,9) box_bd(2)-box_bd(1), box_bd(4) - box_bd(3), box_bd(6)-box_bd(5) write(11,10) !Current boxes are limited to being rectangular do i = 1,3 @@ -479,35 +527,18 @@ module io box_vec(i) = box_bd(2*i) - box_bd(2*i-1) write(11,11) box_vec end do - !We write this as box_mat ori and box_mat current - do i = 1,3 - box_vec(:) = 0.0_dp - box_vec(i) = box_bd(2*i) - box_bd(2*i-1) - write(11,11) box_vec - end do !write the element information if(ele_num > 0) then write(11,12) - - do i = 1, unique_num - write(11,'(3i16)') i, size_ele(unique_index(i))-1, basis_type(1,lat_ele(unique_index(i))) - end do ip = 0 write(11,13) do i = 1, ele_num - !Figure out the ele type - do j = 1, unique_num - if ( (unique_size(j) == size_ele(i)).and.(unique_type(j) == lat_ele(i))) then - etype = j - exit - endif - end do - write(11, '(4i16)') i, etype, 1, basis_type(1,lat_ele(i)) + write(11, '(4i16)') i, basisnum(lat_ele(i)), 2, (size_ele(i)-1) do inod = 1, ng_node(lat_ele(i)) do ibasis = 1, basisnum(lat_ele(i)) ip = ip + 1 - write(11, '(2i16, 3f23.15)') ip, ibasis, r_node(:, ibasis, inod, i) + write(11, '(3i16, 3f23.15)') ip, ibasis, basis_type(ibasis, lat_ele(i)), r_node(:, ibasis, inod, i) end do end do end do @@ -517,7 +548,7 @@ module io if(atom_num /= 0) then write(11,14) do i = 1, atom_num - write(11, '(3i16, 3f23.15)') i, 1, type_atom(i), r_atom(:,i) + write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i) end do end if @@ -803,7 +834,7 @@ module io integer :: i, inod, ibasis, j, k, l, in_eles, in_atoms, ele_types, in_lat_num, in_atom_types, & atom_type_map(100), etype_map(100), etype, lat_type, new_lattice_map(100), & atom_type, stat - real(kind=dp) :: newdisplace(3), r_in(3,1,8), r_in_atom(3) + real(kind=dp) :: newdisplace(3), r_in(3,1,8), r_in_atom(3), atomic_masses(10) character(len=100) :: textholder, in_lattype_map(10) character(len=2) :: atomic_element !First open the file @@ -817,31 +848,25 @@ module io read(11,*) textholder read(11,*) in_eles - !Discard info and read ng_max_node - do i = 1, 5 + !Discard info and read ng_max_node and max_basisnum + do i = 1, 5 read(11,*) textholder end do - read(11,*) max_ng_node - - !Read element types (only needed inside this subroutine) - read(11,*) textholder - read(11,*) ele_types + read(11,*) max_ng_node, max_basisnum !Read in atom num read(11,*) textholder read(11,*) in_atoms - !read in lattice_types and atom types - do i = 1,3 - read(11,*) textholder - end do - read(11,*) in_lat_num, in_atom_types - + !read in atom masses + read(11, *) textholder + read(11, '(a)') textholder + j = tok_count(textholder) + read(textholder, *) (atomic_masses(i), i=1, j) - !Read define atom_types by name - read(11,*) textholder + !Read define atom_types by mass do i = 1, in_atom_types - read(11,*) j, atomic_element + call atommassspecies(atomic_masses(i), atomic_element) call add_atom_type(atomic_element, atom_type_map(i)) end do @@ -849,22 +874,6 @@ module io read(11,*) textholder read(11,*) box_bc - !Disregard useless info - do i = 1, 3 - read(11,*) textholder - end do - - !Read in lat_type map - do i = 1, in_lat_num - read(11,*) j, in_lattype_map(i) - ng_node(lattice_types+i) = 8 !Only cubic elements are currently supported in pyCAC - end do - - !Disregard useless info - do i =1 , 3 - read(11,*) textholder - end do - !Read box boundaries and displace them if necessary read(11,*) temp_box_bd(:) do i = 1, 3 @@ -892,56 +901,13 @@ module io sub_box_bd(:, sub_box_num+1) = temp_box_bd !Read in more useless info - do i = 1, 10 + do i = 1, 9 read(11,*) textholder end do - !Now read the ele_type to size and lat map - do i = 1, ele_types - read(11,*) j, etype_map(i) - etype_map(i) = etype_map(i) + 1 - end do - - - !Now set up the lattice types. In this code it assumes that lattice_type 1 maps to - !atom type 1 because it only allows 1 atom per basis in pyCAC at the moment. - do i = 1, in_lat_num - basisnum(lattice_types+i) = 1 - basis_type(1,lattice_types+i) = atom_type_map(i) - end do - - !Figure out the lattice type maps in case we have repeated lattice_types - k = lattice_types + 1 - new_lattice_map(:) = 0 - new_loop:do i = 1, in_lat_num - old_loop:do j = 1, lattice_types - !First check all the lattice level variables - if ((basisnum(lattice_types+i) == basisnum(j)).and. & - (ng_node(lattice_types+i) == ng_node(j))) then - !Now check the basis level variables - do ibasis =1, basisnum(i) - if(basis_type(ibasis,lattice_types+i) /= basis_type(ibasis,j)) then - cycle old_loop - end if - end do - new_lattice_map(i) = j - cycle new_loop - end if - end do old_loop - new_lattice_map(i) = k - k = k+1 - end do new_loop - - !Read more useless data - read(11,*) textholder - - !set max values and allocate variables - max_basisnum = maxval(basisnum) - max_ng_node = maxval(ng_node) - call grow_ele_arrays(in_eles, in_atoms) - !Now start reading the elements if(in_eles > 0) then + read(11,*) textholder read(11,*) textholder do i = 1, in_eles read(11,*) j, etype, k, lat_type @@ -992,9 +958,9 @@ module io !Internal Variables integer :: i, in_eles, in_atoms, inbtypes(10), lat_type, ia, ie, inod, & - id, type_node, ilat, esize, tag, type + id, type_node, ilat, esize, tag, type, bnum, n, ibasis, ip real(kind=dp) :: newdisplace(3), ra(3), in_lapa, ea, fa(3), va(6), & - ee(1,8), fe(3,1,8), ve(3,1,8), re(3,1,8) + ee(10,8), fe(3,10,8), ve(6,10,8), re(3,10,8) character(len=100) :: textholder, fcc character(len=1000) :: line @@ -1002,7 +968,7 @@ module io open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') !Now initialize some important variables if they aren't defined - if (max_basisnum==0) max_basisnum = 1 + if (max_basisnum==0) max_basisnum = 10 if (max_ng_node==0) max_ng_node=8 fcc="fcc" @@ -1010,6 +976,9 @@ module io read(11, *) textholder read(11, *) textholder + !Read the timestep + read(11, *) textholder, timestep + !Read atom number and element number and grow element arrays by needed amount read(11,*) textholder, in_atoms, textholder, in_eles call grow_ele_arrays(in_eles, in_atoms) @@ -1055,30 +1024,28 @@ module io end if if(in_eles > 0) then - !Add the lattice_types based on the atom types - inbtypes=0 - do i = 1, maxval(type_atom) - inbtypes(1) = i - call lattice_map(1, inbtypes, 8 , 1.0_dp, ilat) !Please check documentation on pycac.out formats - end do !Read element and node headers read(11,*) textholder read(11,*) textholder !read element information, currently only 8 node elements with 1 basis do ie =1, in_eles - read(11,*) tag, lat_type, textholder, textholder, esize - do inod =1, 8 - read(11,*) textholder, textholder, textholder, re(:,1,inod), ee(1,inod), fe(:,1,inod), ve(:,1,inod) + read(11,*) tag, n, bnum, esize + inbtypes(:) = 0 + do inod =1, n*bnum + read(11,*) ip, ibasis, inbtypes(ibasis), re(:,ibasis,ip), ee(ibasis,ip), fe(:,ibasis,ip), ve(:,ibasis,ip) end do + call lattice_map(bnum, inbtypes, n, 1.0_dp, lat_type) call add_element(tag, fcc, esize+1, lat_type, sub_box_num, re) call add_element_data(ele_num, ee, fe, ve) end do end if call set_max_esize return - end subroutine + end subroutine + + subroutine read_lmpcac(file, displace, temp_box_bd) !This subroutine is used to read .cac files which are used with the lammpsCAC format character(len=100), intent(in) :: file diff --git a/src/str.f90 b/src/str.f90 new file mode 100644 index 0000000..2d0ed73 --- /dev/null +++ b/src/str.f90 @@ -0,0 +1,33 @@ +module str + + !this module has some string manipulation commands + public + contains + + pure function tok_count(text) + !counts number of tokens in a string + character(len = *), intent(in) :: text + integer :: tok_count + integer :: i, j + logical :: in_tok + + j = len(trim(adjustl(text))) + in_tok = .false. + tok_count = 0 + do i = 1, j + !This checks if it is a white space character which is the delimiter + if(trim(adjustl(text(i:i))) == ' ') then + !If previously we were in token and the current character is the delimiter + !Then we are no longer in the token + if(in_tok) in_tok = .false. + + !If the character isn't a white space character and we previously weren't in the token then set in_tok + !to true and increment token count + else if(.not.in_tok) then + in_tok = .true. + tok_count = tok_count + 1 + end if + end do + return + end function tok_count +end module str