diff --git a/README.md b/README.md index d038308..036981d 100644 --- a/README.md +++ b/README.md @@ -129,6 +129,11 @@ wrap This will wrap atomic positions back inside the box. Effectively as if periodic boundary conditions are applied so that atoms which exit from one side of the simulation cell enter back in through the other. +###Mode Metric +``` +--metric cmetric nfiles {filepaths} +``` + ## Options Options are additional portions of code which have additional functionality. Options are performed in the order that they appear in the argument list and can be added to any mode. If wanting to use strictly options use `--convert` to specify input and output files. diff --git a/src/Makefile b/src/Makefile index 4d6971d..22aaa6b 100644 --- a/src/Makefile +++ b/src/Makefile @@ -8,9 +8,9 @@ FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -w #FFLAGS=-mcmodel=large -O3 -g #FFLAGS=-mcmodel=large -O0 -g -fbacktrace -fcheck=all -MODES=mode_create.o mode_merge.o mode_convert.o +MODES=mode_create.o mode_merge.o mode_convert.o mode_metric.o OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o opt_deform.o opt_redef_box.o opt_slip_plane.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 +OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o neighbors.o $(MODES) $(OPTIONS) call_option.o sorts.o .SUFFIXES: .SUFFIXES: .c .f .f90 .F90 .o @@ -49,3 +49,4 @@ elements.o : sorts.o $(MODES) $(OPTIONS) subroutines.o io.o : atoms.o box.o $(MODES) main.o : io.o testfuncs.o elements.o mode_create.o $(OPTIONS) $(MODES): subroutines.o +opt_delete.o mode_metric.o : neighbors.o diff --git a/src/box.f90 b/src/box.f90 index faf268c..83d81a8 100644 --- a/src/box.f90 +++ b/src/box.f90 @@ -97,4 +97,10 @@ module box end do return end subroutine in_sub_box + + subroutine reset_box + !This subroutine just resets the box boundary and position + box_bc = "ppp" + box_bd(:) = 0 + end subroutine reset_box end module box diff --git a/src/call_mode.f90 b/src/call_mode.f90 index fb187c1..8d73e25 100644 --- a/src/call_mode.f90 +++ b/src/call_mode.f90 @@ -1,16 +1,16 @@ -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. use mode_create use mode_convert use mode_merge + use mode_metric use parameters implicit none integer, intent(out) :: arg_pos - character(len=100), intent(in) :: mode select case(mode) case('--create') @@ -19,6 +19,8 @@ subroutine call_mode(arg_pos,mode) call convert(arg_pos) case('--merge') call merge(arg_pos) + case('--metric') + call metric(arg_pos) case default print *, "Mode ", trim(adjustl(mode)), " currently not accepted. Please check documentation for ", & "accepted modes and rerun." diff --git a/src/call_option.f90 b/src/call_option.f90 index af63333..a195a0e 100644 --- a/src/call_option.f90 +++ b/src/call_option.f90 @@ -40,6 +40,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('-slip_plane') diff --git a/src/elements.f90 b/src/elements.f90 index 7e658c2..62b56f0 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 @@ -512,6 +516,7 @@ module elements !We go from largest index to smallest index just to make sure that we don't miss !accidentally overwrite values which need to be deleted do i = num, 1, -1 + node_num = node_num - ng_node(lat_ele(sorted_index(i))) if(sorted_index(i) == ele_num) then r_node(:,:,:,sorted_index(i)) = 0.0_dp type_ele(sorted_index(i)) ='' @@ -520,7 +525,6 @@ module elements sbox_ele(sorted_index(i)) = 0 tag_ele(sorted_index(i)) = 0 else - node_num = node_num - ng_node(lat_ele(sorted_index(i))) r_node(:,:,:,sorted_index(i)) = r_node(:,:,:,ele_num) type_ele(sorted_index(i)) = type_ele(ele_num) size_ele(sorted_index(i)) = size_ele(ele_num) @@ -549,8 +553,7 @@ module elements max_bd(:) = -huge(1.0_dp) min_bd(:) = huge(1.0_dp) - - do i = 1, atom_num +do i = 1, atom_num do j = 1, 3 if (r_atom(j,i) > max_bd(j)) max_bd(j) = r_atom(j,i) + tol if (r_atom(j,i) < min_bd(j)) min_bd(j) = r_atom(j,i) - tol @@ -775,4 +778,71 @@ module elements end subroutine + 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 + if (allocated(force_node)) then + deallocate(force_node, virial_node, energy_node) + end if + 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 because of:", allostat + stop + end if + + end if + + if (m > 0) then + if (allocated(force_atom)) then + deallocate(force_atom, virial_atom, energy_atom) + end if + 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 because 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 + + subroutine reset_data + !This function resets all of the data arrays for the elements and atoms + atom_num = 0 + ele_num = 0 + node_num = 0 + end subroutine reset_data + end module elements diff --git a/src/functions.f90 b/src/functions.f90 index 6fed2b7..0a8bc49 100644 --- a/src/functions.f90 +++ b/src/functions.f90 @@ -271,4 +271,124 @@ END FUNCTION StrDnCase norm_dis(1:3) = (rk - rl) norm_dis(4) = norm2(rk-rl) end function + + pure function matinv3(A) result(B) + !! Performs a direct calculation of the inverse of a 3×3 matrix. + real(kind=dp), intent(in) :: A(3,3) !! Matrix + real(kind=dp) :: B(3,3) !! Inverse matrix + real(kind=dp) :: detinv + + if(abs(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)& + - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)& + + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)) < lim_zero) then + B(:,:) = 0 + return + else + ! Calculate the inverse determinant of the matrix + + detinv = 1/(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)& + - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)& + + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)) + + ! Calculate the inverse of the matrix + B(1,1) = +detinv * (A(2,2)*A(3,3) - A(2,3)*A(3,2)) + B(2,1) = -detinv * (A(2,1)*A(3,3) - A(2,3)*A(3,1)) + B(3,1) = +detinv * (A(2,1)*A(3,2) - A(2,2)*A(3,1)) + B(1,2) = -detinv * (A(1,2)*A(3,3) - A(1,3)*A(3,2)) + B(2,2) = +detinv * (A(1,1)*A(3,3) - A(1,3)*A(3,1)) + B(3,2) = -detinv * (A(1,1)*A(3,2) - A(1,2)*A(3,1)) + B(1,3) = +detinv * (A(1,2)*A(2,3) - A(1,3)*A(2,2)) + B(2,3) = -detinv * (A(1,1)*A(2,3) - A(1,3)*A(2,1)) + B(3,3) = +detinv * (A(1,1)*A(2,2) - A(1,2)*A(2,1)) + end if + end function + + pure function transpose3(A) result(B) + !!Transposes matrix A + real(kind=dp), intent(in) :: A(3,3) + real(kind=dp) :: B(3,3) + + integer :: i, j + forall(i =1:3,j=1:3) B(i,j) = A(j,i) + + end function transpose3 + + pure function sqrt3(A) result(B) + !This calculates the square of matrix A fulfilling the equation B*B = A according to + !the algorithm by Franca1989 + + real(kind=dp), intent(in) :: A(3,3) + real(kind=dp) :: B(3,3) + + real(kind=dp) :: Ione, Itwo, Ithree, l, k, phi, Asq(3,3), lambda, Bone, Btwo, Bthree, p + + !Step 1 is calculating the invariants of C + Ione = A(1,1) + A(2,2) + A(3,3) + Asq = matmul(A,A) + Itwo = 0.5_dp *(Ione*Ione - (Asq(1,1) + Asq(2,2) + Asq(3,3))) + Ithree = (A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)& + - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)& + + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)) + + if (Ithree < 0) then + B(:,:)=0.0_dp + return + end if + !Check for an isotropic matrix + k = Ione*Ione - 3*Itwo + if (k < lim_zero) then + lambda = sqrt(Ione/3.0_dp) + B = lambda*identity_mat(3) + else + l = Ione*(Ione*Ione - 9.0_dp/2.0_dp * Itwo) + 27.0_dp/2.0_dp * Ithree + p = l/(k**(1.5_dp)) + + if (p > 1.0 ) then + B(:,:) = 0.0_dp + return + end if + + if ((p< -1).or.(p>1)) then + B(:,:) = 0.0_dp + return + end if + phi = acos(p) + lambda = sqrt(1.0_dp/3.0_dp * (Ione + 2*sqrt(k)*cos(phi/3))) + + !Now calculate invariantes of B + Bthree = sqrt(Ithree) + if((-lambda*lambda + Ione + 2*Ithree/lambda) < 0) then + B(:,:) = 0.0_dp + return + end if + Bone = lambda + sqrt(- lambda*lambda + Ione + 2*Ithree/lambda) + Btwo = (Bone*Bone - Ione)/2.0_dp + + !Now calculate B + if(abs(Bone*Btwo -Bthree) < lim_zero) then + B(:,:) = 0.0_dp + return + end if + B = 1/(Bone*Btwo - Bthree) *(Bone*Bthree*identity_mat(3) + (Bone*Bone - Btwo)*A - matmul(A,A)) + end if + end function sqrt3 + + pure function permutation(i,j,k) result(e) + !Calculates the permutation tensor + integer, intent(in) :: i,j,k + integer :: e + + if ( ((i==1).and.(j==2).and.(k==3)).or. & + ((i==2).and.(j==3).and.(k==1)).or. & + ((i==3).and.(j==1).and.(k==2))) then + e=1 + else if( ((i==2).and.(j==1).and.(k==3)).or. & + ((i==1).and.(j==3).and.(k==2)).or. & + ((i==3).and.(j==2).and.(k==1))) then + e=-1 + else + e=0 + end if + end function permutation + end module functions diff --git a/src/io.f90 b/src/io.f90 index 2787b73..3c9e096 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 @@ -583,34 +583,41 @@ module io 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 does then ask user if they would like to overwrite the file + inquire(file=trim(temp_infile), exist=file_exists) + if (.not.file_exists) then + print *, "The file ", trim(adjustl(filename)), " does not exist. Please input an existing file to read." + stop 3 + end if - !Check to see if file exists, if it does then ask user if they would like to overwrite the file - inquire(file=trim(temp_infile), exist=file_exists) - if (.not.file_exists) then - print *, "The file ", trim(adjustl(filename)), " does not exist. Please input a filename that exists" - read(*,*) temp_infile - cycle + select case(temp_infile(scan(temp_infile,'.',.true.)+1:)) + case('restart', 'mb', 'cac') + infilenum=infilenum+1 + infiles(infilenum) = temp_infile + 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(temp_infile(scan(temp_infile,'.',.true.)+1:)) - case('restart', 'mb', '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: mb, restart." - read(*,*) temp_infile - + select case(trim(adjustl(mode))) + case('--convert','--metric') + infilenum = infilenum+1 + infiles(infilenum) = temp_infile + case default + print *, "Files of type .out cannot be used with mode ", trim(adjustl(mode)) + stop 3 end select - end do + + 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, cac, or out." + stop 3 + end select end subroutine get_in_file subroutine read_in(i, displace, temp_box_bd) - !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine + !This subroutine reads in file i integer, intent(in) :: i real(kind=dp), dimension(3), intent(in) :: displace @@ -624,9 +631,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 +790,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 +969,105 @@ 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 + character(len=1000) :: line + + + 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,'(a)') line(:) + read(line,*) 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) + node_num = node_num + 8 + 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 +1088,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 +1217,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..9f7a484 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 @@ -108,11 +111,13 @@ program main if(bound_called) call def_new_box !Check to make sure a file was passed to be written out and then write out - ! Before building do a check on the file - if (outfilenum == 0) then - argument = 'none' - call get_out_file(argument) + ! Before building do a check on the file + if (trim(adjustl(mode)) /= "--metric") then + if ((outfilenum == 0)) then + argument = 'none' + call get_out_file(argument) + end if + call write_out end if - call write_out end program main diff --git a/src/mode_metric.f90 b/src/mode_metric.f90 new file mode 100644 index 0000000..10838c6 --- /dev/null +++ b/src/mode_metric.f90 @@ -0,0 +1,247 @@ +module mode_metric + !This mode is used to calculate continuum metrics for the j + + use parameters + use io + use elements + use neighbors + + implicit none + + integer :: nfiles + character(len=100) :: metric_type + real(kind=dp) :: rc_off + + !Save reference positions + integer :: np, npreal, nmet + real(kind=dp), allocatable :: r_zero(:,:), r_curr(:,:), met(:,:) + + public + contains + subroutine metric(arg_pos) + !This is the main calling subroutine for the metric code + integer, intent(out) :: arg_pos + character(len=100) :: infile, outfile + + integer :: i, ibasis, inod, np_temp, ppos + real(kind=dp), dimension(6) :: temp_box_bd + + !These are the variables containing the cell list information + integer, dimension(3) :: cell_num + integer, allocatable :: num_in_cell(:,:,:), which_cell(:,:) + integer, allocatable :: cell_list(:,:,:,:) + + !Parse the command arguments + call parse_command(arg_pos) + + !Now read the first file + call read_in(1, (/ 0.0_dp, 0.0_dp, 0.0_dp /), temp_box_bd) + np = atom_num + max_basisnum*max_ng_node*ele_num + allocate(r_zero(3,atom_num+max_basisnum*max_ng_node*ele_num), & + r_curr(3,atom_num+max_basisnum*max_ng_node*ele_num)) + r_zero(:,:) = -huge(1.0_dp) + + !Set up the met variable for the user desired metric + select case(trim(adjustl(metric_type))) + case('def_grad') + allocate(met(9, np)) + case('microrotation') + allocate(met(3,np)) + end select + + !Now set the reference positions + call convert_positions(r_zero, npreal) + + !Now calculate the neighbor list for the reference configuration + call calc_neighbor(5.0_dp, r_zero, np) + + !Reset element and box + call reset_data + call reset_box + + !Now loop over new files + do i = 2, nfiles + call read_in(i, (/ 0.0_dp, 0.0_dp, 0.0_dp /), temp_box_bd) + call convert_positions(r_curr, np_temp) + if (npreal /= np_temp) then + print *, "Error in mode_metric where number of points in ", i, "th file is ", np_temp, " and number of points in" & + // "reference file is", npreal + end if + call calc_metric + !Now create the output file num and write out to xyz format + ppos = scan(trim(infiles(i)),".", BACK= .true.) + if ( ppos > 0 ) then + outfile = infiles(i)(1:ppos)//'xyz' + else + outfile = infiles(i)//'.xyz' + end if + call write_metric_xyz(outfile) + call reset_data + call reset_box + end do + end subroutine metric + + subroutine parse_command(arg_pos) + !This subroutine parses the arguments for mode metric + integer, intent(out) :: arg_pos + + integer :: i, arglen + character(len=100) :: textholder + logical :: file_exists + + !First read the metric to be used + call get_command_argument(2,metric_type,arglen) + if (arglen == 0) stop "Incomplete mode metric command, check documentation" + select case(trim(adjustl(metric_type))) + case("microrotation") + continue + case default + print *, "Mode metric does not accept metric ", metric_type, ". Please select from: microrotation" + stop 3 + end select + + !Now read the cutoff radius + call get_command_argument(3,textholder,arglen) + if (arglen == 0) stop "Incomplete mode metric command, check documentation" + read(textholder, *) rc_off + + !Now read the number of files to read and allocate the variables + call get_command_argument(4, textholder, arglen) + if (arglen == 0) stop "Incomplete mode metric command, check documentation" + read(textholder, *) nfiles + + !Now read the files to be read + do i = 1, nfiles + call get_command_argument(4+i, textholder, arglen) + call get_in_file(textholder) + end do + + arg_pos = 5+nfiles + return + end subroutine parse_command + + subroutine calc_metric + !This subroutine calculates the continuum metric that we require + integer :: i, j, k, nei, ip, jp + real(kind=dp) :: def_grad(3,3), omega(3,3), eta(3,3), rij(3), eta_inv(3,3), ftf(3,3), & + U(3,3), R(3,3), Rskew(3,3), oldrij(3) + + !Loop over all points + do ip = 1, np + eta(:,:) = 0.0_dp + omega(:,:) = 0.0_dp + def_grad(:,:) = 0.0_dp + do jp = 1, nei_num(ip) + !Calculate the neighbor vec in current configuration + nei = nei_list(jp, ip) + rij = r_curr(:,nei) - r_curr(:,ip) + oldrij = r_zero(:,nei) - r_zero(:,ip) + + !Calculate eta and omega + do i = 1,3 + do j = 1,3 + omega(i,j) = omega(i,j) + rij(i) * oldrij(j) + eta(i,j) = eta(i,j) + oldrij(i) * oldrij(j) + end do + end do + end do + + eta_inv=matinv3(eta) + def_grad=matmul(omega,eta_inv) + + select case(trim(adjustl(metric_type))) + case('def_grad') + k = 1 + do i = 1,3 + do j = 1, 3 + met(k, ip) = def_grad(i,j) + end do + end do + case('microrotation') + met(:,ip) = 0.0_dp + if(.not.all(def_grad == 0)) then + !Now calculate microrotation + ftf = matmul(transpose3(def_grad), def_grad) + U = sqrt3(ftf) + if(.not.all(abs(U) < lim_zero)) then + R = matmul(def_grad, matinv3(U)) + Rskew = 0.5_dp * ( R - transpose3(R)) + do k =1,3 + do j = 1,3 + do i = 1,3 + met(k,ip) = met(k,ip) -0.5*permutation(i,j,k)*Rskew(i,j) + end do + end do + end do + end if + end if + end select + end do + return + end subroutine + + subroutine convert_positions(rout, npoints) + !This subroutine just converts current atom and element arrays to a single point based form + real(kind=dp), dimension(3,atom_num+max_ng_node*max_basisnum*ele_num), intent(inout) :: rout + integer, intent(out) :: npoints + + integer :: i, inod, ibasis + + npoints=0 + + if(atom_num > 0) then + do i = 1, atom_num + rout(:,tag_atom(i)) = r_atom(:,i) + npoints = npoints + 1 + end do + end if + + if (ele_num > 0) then + do i = 1, ele_num + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + rout(:, atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis) & + = r_node(:,ibasis,inod,i) + npoints = npoints + 1 + end do + end do + end do + end if + + end subroutine convert_positions + + subroutine write_metric_xyz(outfile) + character(len=100), intent(in) :: outfile + + integer :: i, inod, ibasis + real(kind = dp) :: r(3), eng + open (unit=11, file=trim(adjustl(outfile)), action='write', position='rewind', status = 'replace') + !Write the header + write(11,*) npreal + + select case(metric_type) + case('def_grad') + write(11,*) "type x y z F11 F12 F13 F21 F22 F23 F31 F32 F33" + case('microrotation') + write(11,*) "type x y z micro1 micro2 micro3" + end select + + if(atom_num > 0) then + do i = 1, atom_num + write(11,*) type_atom(i), r_atom(:,i), met(:,tag_atom(i)) + end do + end if + + if (ele_num > 0) then + do i = 1, ele_num + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + write(11,*) basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i), & + met(:, atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis) + end do + end do + end do + end if + end subroutine write_metric_xyz + +end module mode_metric diff --git a/src/neighbors.f90 b/src/neighbors.f90 new file mode 100644 index 0000000..0faebca --- /dev/null +++ b/src/neighbors.f90 @@ -0,0 +1,142 @@ +module neighbors + + use parameters + use elements + use subroutines + use functions + + integer, allocatable :: nei_list(:,:), nei_num(:) + real(kind=dp), allocatable :: init_vec(:,:,:), output(:), microrotation(:,:) + public + contains + + subroutine build_cell_list(numinlist, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell) + !This subroutine builds a cell list based on rc_off + + !----------------------------------------Input/output variables------------------------------------------- + + integer, intent(in) :: numinlist !The number of points within r_list + + real(kind=dp), dimension(3,numinlist), intent(in) :: r_list !List of points to be used for the construction of + !the cell list. + + real(kind=dp), intent(in) :: rc_off ! Cutoff radius which dictates the size of the cells + + integer, dimension(3), intent(inout) :: cell_num !Number of cells in each dimension. + + integer, allocatable, intent(inout) :: num_in_cell(:,:,:) !Number of points within each cell + + integer, allocatable, intent(inout) :: cell_list(:,:,:,:) !Index of points from r_list within each cell. + + integer, dimension(3,numinlist), intent(out) :: which_cell !The cell index for each point in r_list + + !----------------------------------------Begin Subroutine ------------------------------------------- + + integer :: i, j, cell_lim, c(3) + real(kind=dp) :: box_len(3) + integer, allocatable :: resize_cell_list(:,:,:,:) + + !First calculate the number of cells that we need in each dimension + do i = 1,3 + box_len(i) = box_bd(2*i) - box_bd(2*i-1) + cell_num(i) = int(box_len(i)/(rc_off/2))+1 + end do + + !Initialize/allocate variables + cell_lim = 10 + allocate(num_in_cell(cell_num(1),cell_num(2),cell_num(3)), cell_list(cell_lim, cell_num(1), cell_num(2), cell_num(3))) + + !Now place points within cell + do i = 1, numinlist + !Check to see if the current point is a filler point and if so just skip it + if(r_list(1,i) < -huge(1.0_dp)+1) cycle + + !c is the position of the cell that the point belongs to + do j = 1, 3 + c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off/2)) + 1 + end do + + !Place the index in the correct position, growing if necessary + num_in_cell(c(1),c(2),c(3)) = num_in_cell(c(1),c(2),c(3)) + 1 + if (num_in_cell(c(1),c(2),c(3)) > cell_lim) then + allocate(resize_cell_list(cell_lim+10,cell_num(1),cell_num(2),cell_num(3))) + resize_cell_list(1:cell_lim, :, :, :) = cell_list + resize_cell_list(cell_lim+1:, :, :, :) = 0 + call move_alloc(resize_cell_list, cell_list) + end if + + cell_list(num_in_cell(c(1),c(2),c(3)),c(1),c(2),c(3)) = i + which_cell(:,i) = c + end do + + return + end subroutine build_cell_list + + subroutine calc_neighbor(rc_off, r_list, n) + !This function populates the neighbor list in this module + + real(kind=dp), intent(in) :: rc_off + integer, intent(in) :: n + real(kind=dp), dimension(3,n) :: r_list + + integer :: i, c(3), ci, cj, ck, num_nei, nei + !Variables for cell list code + integer, dimension(3) ::cell_num + integer, allocatable :: num_in_cell(:,:,:), cell_list(:,:,:,:) + integer :: which_cell(3,n) + + !First reallocate the neighbor list codes + if (allocated(nei_list)) then + deallocate(nei_list,nei_num) + end if + + allocate(nei_list(100,n),nei_num(n)) + + !Now first pass the position list and and point num to the cell algorithm + call build_cell_list(n, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell) + + !Now loop over every point and find it's neighbors + pointloop: do i = 1, n + + !First check to see if the point is a filler point, if so then skip it + if(r_list(1,i) < -Huge(-1.0_dp)+1) cycle + + !c is the position of the cell that the point + c = which_cell(:,i) + + !Loop over all neighboring cells + do ci = -1, 1, 1 + do cj = -1, 1, 1 + do ck = -1, 1, 1 + !First check to make sure that the neighboring cell exists + if(any((c + (/ ck, cj, ci /)) == 0)) cycle + if( (c(1) + ck > cell_num(1)).or.(c(2) + cj > cell_num(2)).or. & + (c(3) + ci > cell_num(3))) cycle + + do num_nei = 1, num_in_cell(c(1) + ck, c(2) + cj, c(3) + ci) + nei = cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci) + + !Check to make sure the atom isn't the same index as the atom we are checking + !and that the neighbor hasn't already been deleted + if((nei /= i)) then + + !Now check to see if it is in the cutoff radius, if it is add it to the neighbor list for that + !atom and calculate the initial neighbor vector + if (norm2(r_list(:,nei)-r_list(:,i)) < rc_off) then + + nei_num(i) = nei_num(i) + 1 + nei_list(nei_num(i), i) = nei + + end if + end if + end do + end do + end do + end do + + end do pointloop + + return + end subroutine calc_neighbor + +end module neighbors diff --git a/src/opt_delete.f90 b/src/opt_delete.f90 index 09cefb1..c3bcc9e 100644 --- a/src/opt_delete.f90 +++ b/src/opt_delete.f90 @@ -3,6 +3,7 @@ module opt_delete use parameters use subroutines use elements + use neighbors implicit none diff --git a/src/opt_disl.f90 b/src/opt_disl.f90 index 19f8ed8..dc0e72c 100644 --- a/src/opt_disl.f90 +++ b/src/opt_disl.f90 @@ -561,7 +561,6 @@ module opt_disl !Now reset the list for the scanning algorithm delete_num = 0 - !Now scan over all atoms again and find the closest vloop_size number of atoms to the initial atom !that reside on the same plane. If loop_radius is > 0 then we build a circular vacancy cluster if(loop_radius > 0) then do i = 1, atom_num diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 0956677..6fe2896 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -87,7 +87,7 @@ module opt_group continue case default print *, "Select_type ", trim(adjustl(type)), " is not an accept group selection criteria. ", & - "Please select from atoms, nodes, or both." + "Please select from atoms, elements, or both." end select arg_pos = arg_pos + 1 diff --git a/src/parameters.f90 b/src/parameters.f90 index 8677381..c589b4b 100644 --- a/src/parameters.f90 +++ b/src/parameters.f90 @@ -13,5 +13,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 diff --git a/src/subroutines.f90 b/src/subroutines.f90 index 1015968..b272c16 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -198,66 +198,6 @@ module subroutines end do end subroutine - subroutine build_cell_list(numinlist, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell) - !This subroutine builds a cell list based on rc_off - - !----------------------------------------Input/output variables------------------------------------------- - - integer, intent(in) :: numinlist !The number of points within r_list - - real(kind=dp), dimension(3,numinlist), intent(in) :: r_list !List of points to be used for the construction of - !the cell list. - - real(kind=dp), intent(in) :: rc_off ! Cutoff radius which dictates the size of the cells - - integer, dimension(3), intent(inout) :: cell_num !Number of cells in each dimension. - - integer, allocatable, intent(inout) :: num_in_cell(:,:,:) !Number of points within each cell - - integer, allocatable, intent(inout) :: cell_list(:,:,:,:) !Index of points from r_list within each cell. - - integer, dimension(3,numinlist), intent(out) :: which_cell !The cell index for each point in r_list - - !----------------------------------------Begin Subroutine ------------------------------------------- - - integer :: i, j, cell_lim, c(3) - real(kind=dp) :: box_len(3) - integer, allocatable :: resize_cell_list(:,:,:,:) - - !First calculate the number of cells that we need in each dimension - do i = 1,3 - box_len(i) = box_bd(2*i) - box_bd(2*i-1) - cell_num(i) = int(box_len(i)/(rc_off/2))+1 - end do - - !Initialize/allocate variables - cell_lim = 10 - allocate(num_in_cell(cell_num(1),cell_num(2),cell_num(3)), cell_list(cell_lim, cell_num(1), cell_num(2), cell_num(3))) - - !Now place points within cell - do i = 1, numinlist - !c is the position of the cell that the point belongs to - do j = 1, 3 - c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off/2)) + 1 - end do - - !Place the index in the correct position, growing if necessary - num_in_cell(c(1),c(2),c(3)) = num_in_cell(c(1),c(2),c(3)) + 1 - if (num_in_cell(c(1),c(2),c(3)) > cell_lim) then - allocate(resize_cell_list(cell_lim+10,cell_num(1),cell_num(2),cell_num(3))) - resize_cell_list(1:cell_lim, :, :, :) = cell_list - resize_cell_list(cell_lim+1:, :, :, :) = 0 - call move_alloc(resize_cell_list, cell_list) - end if - - cell_list(num_in_cell(c(1),c(2),c(3)),c(1),c(2),c(3)) = i - which_cell(:,i) = c - end do - - return - end subroutine build_cell_list - - subroutine check_right_ortho(ori, isortho, isrighthanded) !This subroutine checks whether provided orientations in the form: ! | x1 x2 x3 | @@ -301,4 +241,5 @@ module subroutines return end subroutine check_right_ortho + end module subroutines