diff --git a/src/Makefile b/src/Makefile index feb933f..c58e37c 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,7 +1,7 @@ FC=ifort FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays #FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays -MODES=mode_create.o mode_merge.o mode_convert.o +MODES=mode_create.o mode_merge.o mode_convert.o OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o opt_deform.o opt_redef_box.o OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES) $(OPTIONS) call_option.o sorts.o diff --git a/src/call_mode.f90 b/src/call_mode.f90 index fb187c1..7eae059 100644 --- a/src/call_mode.f90 +++ b/src/call_mode.f90 @@ -1,4 +1,4 @@ -subroutine call_mode(arg_pos,mode) +subroutine call_mode(arg_pos) !This code is used to parse the command line argument for the mode information and calls the required !mode module. @@ -10,7 +10,6 @@ subroutine call_mode(arg_pos,mode) implicit none integer, intent(out) :: arg_pos - character(len=100), intent(in) :: mode select case(mode) case('--create') diff --git a/src/call_option.f90 b/src/call_option.f90 index d55e1ee..aeb17b5 100644 --- a/src/call_option.f90 +++ b/src/call_option.f90 @@ -39,6 +39,8 @@ subroutine call_option(option, arg_pos) call run_delete(arg_pos) case('-set_cac') arg_pos=arg_pos +3 + case('-set_types') + arg_pos = arg_pos + 3 + atom_types case('-redef_box') call redef_box(arg_pos) case default diff --git a/src/elements.f90 b/src/elements.f90 index 6df859b..a12d263 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -13,6 +13,8 @@ module elements character(len=100), allocatable :: type_ele(:) !Element type 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(:,:,:) integer, save :: ele_num !Number of elements integer, save :: node_num !Total number of nodes @@ -22,6 +24,8 @@ module elements integer, allocatable :: sbox_atom(:), tag_atom(:) real(kind =dp),allocatable :: r_atom(:,:) !atom position integer :: atom_num=0 !Number of atoms + !Atom result data structures information + real(kind=8), allocatable :: force_atom(:,:), virial_atom(:,:), energy_atom(:) !Mapping atom type to provided name character(len=2), dimension(10) :: type_to_name @@ -669,17 +673,17 @@ module elements esize = size_ele(ie) select case(iface) case(1) - pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp /) + pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, -10.0_dp**(-2.0_dp) /) case(2) - pos = (/ real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /) + pos = (/ real(esize-1,dp)/2.0_dp, -10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp /) case(3) - pos = (/ (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) + pos = (/ (esize-1)+10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) case(4) - pos = (/ real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /) + pos = (/ real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp /) case(5) - pos = (/ -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) + pos = (/ -10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) case(6) - pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp /) + pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**(-2.0_dp) /) end select !Now transform it to real space and adjust it to the position of the element in the first node. @@ -736,4 +740,57 @@ module elements end subroutine lattice_map + subroutine alloc_dat_arrays(n,m) + !This subroutine used to provide initial allocation for the atom and element data arrays + integer, intent(in) :: n,m !n-size of element arrays, m-size of atom arrays + integer :: allostat + + !Allocate element arrays + if (n > 0) then + allocate(force_node(3,max_basisnum, max_ng_node, n), & + virial_node(6,max_basisnum, max_ng_node, n), & + energy_node(max_basisnum,max_ng_node,n), & + stat=allostat) + if(allostat > 0) then + print *, "Error allocating element data arrays in mode_metric becaus of:", allostat + stop + end if + + end if + + if (m > 0) then + allocate(force_atom(3, m), & + virial_atom(6, m), & + energy_atom(m), & + stat=allostat) + if(allostat > 0) then + print *, "Error allocating atom data arrays in mode_metric becaus of:", allostat + stop + end if + end if + + end subroutine + + subroutine add_atom_data(ia, eng, force, virial) + !Function which sets the atom data arrays + integer, intent(in) :: ia + real(kind=dp), intent(in) :: eng, force(3), virial(6) + + energy_atom(ia) = eng + force_atom(:,ia) = force(:) + virial_atom(:,ia) = virial(:) + return + end subroutine add_atom_data + + subroutine add_element_data(ie, eng, force, virial) + !Function which sets the element data arrays + integer, intent(in) :: ie + real(kind=dp), intent(in) :: eng(max_basisnum, max_ng_node), & + force(3,max_basisnum, max_ng_node), & + virial(6,max_basisnum,max_ng_node) + energy_node(:,:,ie) = eng + force_node(:,:,:,ie) = force + virial_node(:,:,:,ie) = virial + return + end subroutine add_element_data end module elements diff --git a/src/io.f90 b/src/io.f90 index 7bbf3c2..e1ec449 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -106,7 +106,7 @@ module io 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" + " is not accepted for writing. Please select from: xyz, lmp, vtk, mb, restart, cac and try again" stop end select @@ -599,9 +599,24 @@ module io infilenum=infilenum+1 infiles(infilenum) = temp_infile exit + case('out') + if(atom_types == 0) then + print *, "Please run -set_types command prior to running code requiring reading in pycac_*.out files" + stop 3 + end if + select case(trim(adjustl(mode))) + case('--convert','--metric') + infilenum = infilenum+1 + infiles(infilenum) = temp_infile + exit + case default + print *, "Files of type .out cannot be used with mode ", trim(adjustl(mode)) + stop 3 + end select + case default print *, "File type: ", trim(temp_infile(scan(temp_infile,'.',.true.):)), "not currently accepted. ", & - "please input a filename with extension from following list: mb, restart." + "please input a filename with extension from following list: mb, restart, cac, or out." read(*,*) temp_infile end select @@ -624,9 +639,11 @@ module io call read_pycac(infiles(i), displace, temp_box_bd) case('cac') call read_lmpcac(infiles(i), displace, temp_box_bd) + case('out') + call read_pycac_out(infiles(i), displace, temp_box_bd) case default print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), & - " is not accepted for writing. Please select from: mb and try again" + " is not accepted for reading. Please select from: mb,restart,cac,out and try again" stop end select @@ -781,7 +798,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 + atom_type, stat real(kind=dp) :: newdisplace(3), r_in(3,1,8), r_in_atom(3) character(len=100) :: textholder, in_lattype_map(10) character(len=2) :: atomic_element @@ -960,6 +977,102 @@ module io end if end subroutine read_pycac + subroutine read_pycac_out(file, displace, temp_box_bd) + !This subroutine reads in the pyCAC dump file + + + !Arguments + character(len=100), intent(in) :: file + real(kind=dp), dimension(3), intent(in) :: displace + real(kind=dp), dimension(6), intent(out) :: temp_box_bd + + !Internal Variables + integer :: i, in_eles, in_atoms, inbtypes(10), lat_type, ia, ie, inod, & + id, type_node, ilat, esize, tag, type + 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) + character(len=100) :: textholder, fcc + + + 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_ng_node==0) max_ng_node=8 + fcc="fcc" + + !Skip header comment lines + read(11, *) textholder + read(11, *) textholder + + !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) + call alloc_dat_arrays(in_eles, in_atoms) + + !Read boundary information + read(11,*) textholder, box_bc(1:1), box_bc(2:2), box_bc(3:3), temp_box_bd(:) + do i = 1, 3 + if (abs(displace(i)) > lim_zero) then + newdisplace(i) = displace(i) - temp_box_bd(2*i-1) + else + newdisplace(i)=displace(i) + end if + temp_box_bd(2*i-1) = temp_box_bd(2*i-1) + newdisplace(i) + temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i) + end do + + call grow_box(temp_box_bd) + + !Allocate sub_box boundaries + if (sub_box_num == 0) then + call alloc_sub_box(1) + else + call grow_sub_box(1) + end if + + !Because orientations and other needed sub_box information isn't really + !present within the .cac file we just default a lot of it. + sub_box_ori(:,:,sub_box_num+1) = identity_mat(3) + sub_box_bd(:, sub_box_num+1) = temp_box_bd + sub_box_num = sub_box_num + 1 + + if(in_atoms > 0 ) then + !Read atom header + read(11,*) textholder + do ia = 1, in_atoms + read(11,*) tag, type, ra(:), ea, fa(:), va(:) + call add_atom(tag, type, sub_box_num, ra) + call add_atom_data(atom_num, ea, fa, va) + end do + + 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) + end do + 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 + 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 @@ -980,7 +1093,7 @@ module io !Open the file open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') - !Now initialiaze some important variables if they aren't defined + !Now initialize some important variables if they aren't defined if (max_basisnum==0) max_basisnum = 10 if (max_ng_node==0) max_ng_node=8 @@ -1109,4 +1222,25 @@ module io print *, in_lattice_type end subroutine set_cac + + subroutine set_types(apos) + !This code + integer, intent(in) :: apos + integer :: i, j,arglen, arg_pos, ntypes + + character(len=100) :: textholder + + arg_pos = apos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) stop "Missing numtypes in io" + read(textholder,*) ntypes + + do i=1,ntypes + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + call add_atom_type(textholder, j) + end do + + return + end subroutine set_types end module io diff --git a/src/main.f90 b/src/main.f90 index 53fa4fa..b305138 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -63,6 +63,8 @@ program main case('-set_cac') call set_cac(i) + case('-set_types') + call set_types(i) end select end do !Determine if a mode is being used and what it is. The first argument has to be the mode @@ -71,7 +73,8 @@ program main argument = trim(adjustl(argument)) if (argument(1:2) == '--') then - call call_mode(end_mode_arg, argument) + mode = argument + call call_mode(end_mode_arg) end if !Check to make sure a mode was called diff --git a/src/parameters.f90 b/src/parameters.f90 index f261552..825245d 100644 --- a/src/parameters.f90 +++ b/src/parameters.f90 @@ -12,5 +12,7 @@ module parameters !Numeric constants real(kind=dp), parameter :: pi = 3.14159265358979323846_dp - + + !Mode type that is being called + character(len=100) :: mode end module parameters