Add periodic boundaries to metric calculation and energy to vtk out

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
 

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)
+                rij = r_curr(:,nei) - r_curr(:,ip)  +periodvec(:,jp,ip)*box_len 
-                oldrij = r_zero(:,nei) - r_zero(:,ip)
+                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

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)
+                        do num_nei = 1, num_in_cell(cn(1),cn(2),cn(3))
-                            nei = cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci)
+                            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