From cffb460e09cb8959204b5e7ae6e60ce81628c540 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Tue, 24 May 2022 06:52:30 -0400 Subject: [PATCH] Add periodic boundaries to metric calculation and energy to vtk out --- src/io.f90 | 24 +++++++++++++++++++++ src/mode_metric.f90 | 24 +++++++++++++++++---- src/neighbors.f90 | 51 +++++++++++++++++++++++++++++++++++---------- 3 files changed, 84 insertions(+), 15 deletions(-) diff --git a/src/io.f90 b/src/io.f90 index 6519862..3bc511a 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -436,6 +436,7 @@ module io 16 format(/'POINT_DATA', i16) 17 format('SCALARS weight float', / & 'LOOKUP_TABLE default') + 18 format('SCALARS atom_type float', / & 'LOOKUP_TABLE default') @@ -445,6 +446,9 @@ module io 21 format('SCALARS esize float', /& 'LOOKUP_TABLE default') +22 format('SCALARS energy float', /& + 'LOOKUP_TABLE default') + !First we write the vtk file containing the atoms open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind') @@ -467,6 +471,14 @@ module io do i = 1, atom_num write(11, '(i16)') type_atom(i) end do + + !Write the atom data + if(allocated(force_atom)) then + write(11,22) + do i =1, atom_num + write(11, "(f23.15)") energy_atom(i) + end do + end if close(11) !Now we write the vtk file for the elements @@ -499,6 +511,18 @@ module io do i = 1, ele_num write(11, '(i16)') size_ele(i) end do + + if(allocated(force_node)) then + write(11,16) node_num + write(11,22) + do i = 1, ele_num + do inod=1, ng_node(lat_ele(i)) + do ibasis=1, basisnum(lat_ele(i)) + write(11,'(f23.15)') energy_node(ibasis, inod, i) + end do + end do + end do + end if close(11) end subroutine diff --git a/src/mode_metric.f90 b/src/mode_metric.f90 index dab29f5..e6f6567 100644 --- a/src/mode_metric.f90 +++ b/src/mode_metric.f90 @@ -10,7 +10,7 @@ module mode_metric integer :: nfiles character(len=100) :: metric_type - real(kind=dp) :: rc_off + real(kind=dp) :: rc_off, old_box_len(3) !Save reference positions integer :: np, npreal, nmet @@ -57,6 +57,10 @@ module mode_metric print *, "Getting initial neighbor list" call calc_neighbor(rc_off, r_zero, np) + !Save box lengths for applying periodic boundaries + do i=1,3 + old_box_len(i) = box_bd(2*i) - box_bd(2*i-1) + end do !Reset element and box call reset_data call reset_box @@ -128,8 +132,11 @@ module mode_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) + U(3,3), R(3,3), Rskew(3,3), oldrij(3), box_len(3) + do i = 1, 3 + box_len(i) = box_bd(2*i) - box_bd(2*i-1) + end do !Loop over all points do ip = 1, np eta(:,:) = 0.0_dp @@ -138,8 +145,17 @@ module mode_metric 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) + rij = r_curr(:,nei) - r_curr(:,ip) +periodvec(:,jp,ip)*box_len + if (norm2(rij) > 5*rc_off) then + do j = 1, 3 + if (r_curr(j,nei) > r_curr(j,ip)) then + rij(j) = rij(j) - box_len(j) + else if(r_curr(j,nei) < r_curr(j,ip)) then + rij(j) = rij(j) + box_len(j) + end if + end do + end if + oldrij = r_zero(:,nei) - r_zero(:,ip) + periodvec(:,jp,ip)*old_box_len !Calculate eta and omega do i = 1,3 diff --git a/src/neighbors.f90 b/src/neighbors.f90 index 6e81c38..c415eae 100644 --- a/src/neighbors.f90 +++ b/src/neighbors.f90 @@ -4,10 +4,12 @@ module neighbors use elements use subroutines use functions + use box - integer, allocatable :: nei_list(:,:), nei_num(:), nn(:) + integer, allocatable :: nei_list(:,:), nei_num(:), nn(:), periodvec(:,:,:) 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) @@ -82,44 +84,68 @@ module neighbors integer, intent(in) :: n real(kind=dp), dimension(3,n) :: r_list - integer :: i, c(3), ci, cj, ck, num_nei, nei + integer :: i, j, c(3),cn(3), ci, cj, ck, num_nei, nei, v(3), period_dir(3) !Variables for cell list code integer, dimension(3) ::cell_num integer, allocatable :: num_in_cell(:,:,:), cell_list(:,:,:,:) integer :: which_cell(3,n) + real(kind=dp) :: r(3), box_len(3) + logical :: period_bd(3) !First reallocate the neighbor list codes if (allocated(nei_list)) then - deallocate(nei_list,nei_num) + deallocate(nei_list,nei_num, periodvec) end if - allocate(nei_list(100,n),nei_num(n)) + allocate(nei_list(100,n),nei_num(n), periodvec(3,100,n)) nei_list=0 + periodvec=0 nei_num=0 !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) + do i=1, 3 + if (box_bc(i:i) == 'p') then + period_bd(i) = .true. + else + period_bd(i)=.false. + end if + box_len(i) = box_bd(2*i) - box_bd(2*i-1) + end do !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) < box_bd(1)) cycle - !c is the position of the cell that the point + !c is the position of the cell that the point belongs to 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 + !First try to apply periodic boundaries + v=(/ ck, cj, ci /) + cn=0 + period_dir=0 + do j=1, 3 + if ((c(j) + v(j) == 0).and.period_bd(j)) then + cn(j) = cell_num(j) + period_dir(j) = 1 + else if ((c(j) + v(j) > cell_num(j)).and.period_bd(j)) then + cn(j) = 1 + period_dir(j) = -1 + + else if ((c(j)+v(j) >= 1) .and. (c(j)+v(j) <= cell_num(j))) then + cn(j) = c(j) + v(j) + end if + end do 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) + do num_nei = 1, num_in_cell(cn(1),cn(2),cn(3)) + nei = cell_list(num_nei,cn(1),cn(2),cn(3)) !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 @@ -127,10 +153,13 @@ module neighbors !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 + + r = r_list(:,nei) + period_dir*box_len + if (norm2(r-r_list(:,i)) < rc_off) then nei_num(i) = nei_num(i) + 1 nei_list(nei_num(i), i) = nei + periodvec(:,nei_num(i),i) = period_dir end if end if