From 6e085176975f908e1e5421cceb76892b8edc4fed Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Tue, 20 Oct 2020 01:03:31 -0400 Subject: [PATCH] Working continuum metric calculation code --- src/Makefile | 5 +- src/call_mode.f90 | 3 + src/elements.f90 | 10 +- src/functions.f90 | 120 ++++++++++++++++++++++ src/io.f90 | 4 +- src/main.f90 | 12 ++- src/mode_metric.f90 | 245 ++++++++++++++++++++++++++++++++++++++++++++ src/neighbors.f90 | 142 +++++++++++++++++++++++++ 8 files changed, 531 insertions(+), 10 deletions(-) create mode 100644 src/mode_metric.f90 create mode 100644 src/neighbors.f90 diff --git a/src/Makefile b/src/Makefile index c58e37c..6106104 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,9 +1,9 @@ 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 mode_metric.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 +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 @@ -42,3 +42,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/call_mode.f90 b/src/call_mode.f90 index 7eae059..8d73e25 100644 --- a/src/call_mode.f90 +++ b/src/call_mode.f90 @@ -5,6 +5,7 @@ subroutine call_mode(arg_pos) use mode_create use mode_convert use mode_merge + use mode_metric use parameters implicit none @@ -18,6 +19,8 @@ subroutine call_mode(arg_pos) 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/elements.f90 b/src/elements.f90 index ab1a60e..1254a43 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -747,24 +747,30 @@ module elements !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 becaus of:", allostat + 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 becaus of:", allostat + print *, "Error allocating atom data arrays in mode_metric because of:", allostat stop end if end if 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 556ce1f..9c02ba1 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -984,6 +984,7 @@ module io 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') @@ -1033,7 +1034,8 @@ module io !Read atom header read(11,*) textholder do ia = 1, in_atoms - read(11,*) tag, type, ra(:), ea, fa(:), va(:) + 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 diff --git a/src/main.f90 b/src/main.f90 index b305138..9f7a484 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -111,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..279e710 --- /dev/null +++ b/src/mode_metric.f90 @@ -0,0 +1,245 @@ +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), allocatable :: met(:,:) + + !Save reference positions + integer :: np, npreal, nmet + real(kind=dp), allocatable :: r_zero(:,:), r_curr(:,:) + + 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 + print *,np + 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 number of files to read and allocate the variables + call get_command_argument(3, 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(3+i, textholder, arglen) + call get_in_file(textholder) + end do + + arg_pos = 4+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 + print *, atom_num + max_ng_node*max_basisnum*ele_num + print *, rout(:,1) + + 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