Changes to how the adjustment to nodal positions is performed for lammpscac output

src/elements.f90

@@ -9,7 +9,7 @@ module elements
 
     !Data structures used to represent the CAC elements. Each index represents an element
     character(len=100), allocatable :: type_ele(:) !Element type 
-    integer, allocatable :: size_ele(:), lat_ele(:) !Element siz
+    integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:) !Element size
     real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array
 
     integer, save :: ele_num !Number of elements
@@ -39,6 +39,7 @@ module elements
     !This can be easily increased with no change to efficiency
     integer :: max_basisnum, basisnum(10) !Max basis atom number, number of basis atoms in each lattice type 
     integer :: basis_type(10,10)
+    real(kind=dp) :: lapa(10)
 
     public 
     contains
@@ -140,7 +141,7 @@ module elements
         
         !Allocate element arrays
         if(n > 0) then 
-            allocate(type_ele(n), size_ele(n), lat_ele(n), r_node(3,max_basisnum, max_ng_node,n), &
+            allocate(type_ele(n), size_ele(n), lat_ele(n), sbox_ele(n), r_node(3,max_basisnum, max_ng_node,n), &
                        stat=allostat)
             if(allostat > 0) then 
                 print *, "Error allocating element arrays in elements.f90 because of: ", allostat
@@ -179,12 +180,17 @@ module elements
             allocate(temp_int(n+ele_num+buffer_size))
             temp_int(1:ele_size) = lat_ele
             temp_int(ele_size+1:) = 0
-            call move_alloc(temp_int(1:ele_size), lat_ele)
+            call move_alloc(temp_int, lat_ele)
 
             allocate(temp_int(n+ele_num+buffer_size))
             temp_int(1:ele_size) = size_ele
             temp_int(ele_size+1:) = 0
-            call move_alloc(temp_int(1:ele_size), size_ele)
+            call move_alloc(temp_int, size_ele)
+
+            allocate(temp_int(n+ele_num+buffer_size))
+            temp_int(1:ele_size) = lat_ele
+            temp_int(ele_size+1:) = 0
+            call move_alloc(temp_int, sbox_ele)
 
             allocate(char_temp(n+ele_num+buffer_size))
             char_temp(1:ele_size) = type_ele
@@ -210,9 +216,9 @@ module elements
         end if
     end subroutine
 
-    subroutine add_element(type, size, lat, r)
+    subroutine add_element(type, size, lat, sbox, r)
         !Subroutine which adds an element to the element arrays
-        integer, intent(in) :: size, lat
+        integer, intent(in) :: size, lat, sbox
         character(len=100), intent(in) :: type
         real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node)
 
@@ -220,6 +226,7 @@ module elements
         type_ele(ele_num) = type
         size_ele(ele_num) = size
         lat_ele(ele_num) = lat
+        sbox_ele(ele_num) = sbox
        r_node(:,:,:,ele_num) = r(:,:,:)
        node_num = node_num + ng_node(lat)
 

src/io.f90

@@ -204,7 +204,7 @@ module io
         !This subroutine writes out a .lmp style dump file
         character(len=100), intent(in) :: file
         integer :: write_num, i, inod, ibasis
-        real(kind=dp) :: mass
+        real(kind=dp) :: mass, fcc_adjust(3,8), local_adjust(3,8), rout(3)
 
 1       format(i16, '   Eight_Node', 4i16)
 2       format(i16, '   Atom', 4i16)
@@ -242,16 +242,32 @@ module io
         write(11, '(a)') 'CAC Elements'
         write(11, '(a)') ' '
 
+        !Set up the nodal adjustment variables for all the different element types. This adjusts the node centers
+        !from the center of the unit cell (as formulated in this code) to the corners of the unit cells
+         do inod = 1, 8
+            do i = 1,3
+                if(is_equal(cubic_cell(i, inod),0.0_dp)) then 
+                    fcc_adjust(i,inod) = -0.5_dp
+                else 
+                    fcc_adjust(i, inod) = 0.5_dp
+                end if
+            end do
+        end do
+        fcc_adjust = matmul(fcc_mat, fcc_adjust)
+
         !Write element nodal positions
         do i = 1, ele_num
             select case(trim(adjustl(type_ele(i))))
             case('fcc')
+                !Now orient the current adjustment vector to the correct orientation
+                local_adjust = matmul(sub_box_ori(:,:,sbox_ele(i)), fcc_adjust) * lapa(lat_ele(i))
                 !The first entry is the element specifier 
                 write(11,1) i, basisnum(lat_ele(i)), size_ele(i), size_ele(i), size_ele(i)
                 do ibasis = 1, basisnum(lat_ele(i))
                     do inod = 1, 8
                         !Nodal information for every node
-                        write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i)
+                        rout = r_node(:,ibasis,inod,i) + local_adjust(:,inod)
+                        write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), rout
                     end do
                 end do
             end select
@@ -486,6 +502,8 @@ module io
         write(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types)
         !Now for every lattice type write the basis atom types
         write(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = 1, lattice_types)
+        !Now for every lattice type write the lattice parameters
+        write(11,*) (lapa(i), i = 1, lattice_types)
 
         !Now write the numbers of elements and atoms
         write(11,*) atom_num, ele_num
@@ -498,7 +516,7 @@ module io
         !Write out the elements, this is written in two stages, one line for the element and then 1 line for 
         !every basis at every node
         do i = 1, ele_num
-            write(11, *) i, lat_ele(i), size_ele(i), type_ele(i)
+            write(11, *) i, lat_ele(i), size_ele(i), sbox_ele(i), type_ele(i)
             do inod = 1, ng_node(lat_ele(i))
                 do ibasis =1, basisnum(lat_ele(i))
                     write(11,*) inod, ibasis, r_node(:, ibasis, inod, i)
@@ -582,7 +600,7 @@ module io
         real(kind = dp), dimension(6), intent(out) :: temp_box_bd
 
         integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles, new_atom_types, &
-                   new_type_to_type(10), new_lattice_types
+                   new_type_to_type(10), new_lattice_types, sbox
         character(len=100) :: etype
         real(kind=dp) ::  r(3), newdisplace(3)
         real(kind=dp), allocatable :: r_innode(:,:,:)
@@ -623,8 +641,6 @@ module io
         sub_box_array_bd(:,1,sub_box_num+1:) = sub_box_array_bd(:,1,sub_box_num+1:) + atom_num
         sub_box_array_bd(:,2,sub_box_num+1:) = sub_box_array_bd(:,2,sub_box_num+1:) + ele_num
 
-        sub_box_num = sub_box_num + n
-
         !Read in the number of atom types and all their names
         read(11, *) new_atom_types, (new_type_to_name(i), i = 1, new_atom_types)
         !Now fit these into the global list of atom types, after this new_type_to_type is the actual global 
@@ -646,6 +662,8 @@ module io
                 basis_type(i,j) = new_type_to_type(basis_type(i,j))
             end do 
         end do 
+        !Read the lattice parameters for every lattice type
+        read(11,*) (lapa(i), i = lattice_types+1, lattice_types+new_lattice_types)
         !Read number of elements and atoms and allocate arrays
         read(11, *) in_atoms, in_eles
         call grow_ele_arrays(in_eles, in_atoms)
@@ -659,7 +677,7 @@ module io
 
         !Read the elements
         do i = 1, in_eles
-            read(11, *) n, type, size, etype
+            read(11, *) n, type, size, sbox, etype
             do inod = 1, ng_node(type)
                 do ibasis =1, basisnum(type)
                     read(11,*) j, k, r_innode(:, ibasis, inod)
@@ -667,7 +685,7 @@ module io
                 end do
             end do
             type = type + lattice_types
-            call add_element(etype, size, type, r_innode)
+            call add_element(etype, size, type, sbox+n, r_innode)
         end do
 
         !Close the file being read
@@ -676,5 +694,9 @@ module io
         !Only increment the lattice types if there are elements, if there are no elements then we 
         !just overwrite the arrays
         if(in_eles > 0) lattice_types = lattice_types + new_lattice_types
+
+
+        sub_box_num = sub_box_num + n
+
     end subroutine read_mb
 end module io

src/main.f90

@@ -28,17 +28,6 @@ program main
     ! Command line parsing
     arg_num = command_argument_count()
 
-    !First check if we are writing out to lammpscac format by looping over all arguments
-    do i = 1, arg_num
-        call get_command_argument(i, argument)
-        select case(argument(scan(argument,'.',.true.)+1:))
-        case('cac')
-            lmpcac = .true.
-        case default
-            continue
-        end select
-    end do 
-
     !Determine if a mode is being used and what it is. The first argument has to be the mode
     !if a mode is being used
     call get_command_argument(1, argument)

src/mode_create.f90

@@ -12,7 +12,7 @@ module mode_create
 
     character(len=100) :: name, element_type
     real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), &
-                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3), adjustVar(3,8)
+                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3)
     integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), &
                basis_pos(3,10)
     logical :: dup_flag, dim_flag
@@ -105,7 +105,7 @@ module mode_create
                 box_bd(2*i) = maxval(r_node_temp(i,:,:))
                 box_bd(2*i-1) = origin(i)
             end do
-            call add_element(element_type, esize, 1, r_node_temp)
+            call add_element(element_type, esize, 1, 1, r_node_temp)
         end if
         
         !If we passed the dup_flag or dim_flag then we have to convert the lattice points and add them to the atom/element arrays
@@ -141,7 +141,7 @@ module mode_create
                             r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis)
                         end do
                     end do 
-                    call add_element(element_type, esize, 1, r_node_temp)
+                    call add_element(element_type, esize, 1, 1, r_node_temp)
                 end do 
             end if
         end if
@@ -199,7 +199,20 @@ module mode_create
                     call get_command_argument(arg_pos, orient_string, arglen)
                     if (arglen==0) STOP "Missing orientation in orient command of mode create"
                     arg_pos = arg_pos+1
-                    call parse_ori_vec(orient_string, orient(i,:))
+                    ori_pos=2
+                    do j = 1,3
+                        if (orient_string(ori_pos:ori_pos) == '-') then 
+                            ori_pos = ori_pos + 1
+                            read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
+                            if (stat>0) STOP "Error reading orient value"
+                            orient(i,j) = -orient(i,j)
+                            ori_pos = ori_pos + 1
+                        else
+                            read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
+                            if(stat>0) STOP "Error reading orient value"
+                            ori_pos=ori_pos + 1
+                        end if
+                    end do
                 end do
 
 
@@ -234,6 +247,7 @@ module mode_create
                 exit
             end select
         end do
+
         !Calculate the lattice periodicity length in lattice units
         do i = 1, 3
             lattice_space(i) = norm2(orient(i,:))
@@ -264,8 +278,9 @@ module mode_create
         !Now normalize the orientation matrix
         orient =  matrix_normal(orient,3)
 
-        !Set lattice_num to 1 
+        !Set lattice_num to 1 and add the lattice_parameter to the elements module lattice paramter variable
         lattice_types = 1
+        lapa(1) = lattice_parameter
 
         !If we haven't defined a basis then define the basis and add the default basis atom type and position
         if (basisnum(1) == 0) then
@@ -286,7 +301,7 @@ module mode_create
         !Internal variables
         integer :: i,  inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz,  numlatpoints, ele(3,8), rzero(3), &
                    vlat(3), temp_lat(3,8), m, n, o
-        real(kind=dp) :: v(3), temp_nodes(3,1,8), adjustVar(3,8)
+        real(kind=dp) :: v(3), temp_nodes(3,1,8)
         real(kind=dp), allocatable :: resize_lat_array(:,:)
         logical, allocatable :: lat_points(:,:,:)
         logical :: node_in_bd(8)
@@ -294,23 +309,6 @@ module mode_create
         !Do some value initialization
         max_esize = esize
 
-        !If we are writing out to lammpscac format then we have to adjust the nodal positions
-
-        if(lmpcac) then 
-             do inod = 1, 8
-                do i = 1,3
-                    if(is_equal(cubic_cell(i, inod),0.0_dp)) then 
-                        adjustVar(i,inod) = -0.5_dp
-                    else 
-                        adjustVar(i, inod) = 0.5_dp
-                    end if
-                end do
-            end do
-
-            adjustVar(:,1:8) = matmul(orient,matmul(fcc_mat,adjustVar(:,1:8)))
-        else
-            adjustVar(:,:)=0.0_dp
-        end if
         !First find the bounding lattice points (min and max points for the box in each dimension)
         numlatpoints = 1
         do i = 1, 3
@@ -328,6 +326,7 @@ module mode_create
             continue
         end select
 
+
         !Loop over all of lattice points within the boundary, we choose between two loops. One for the atomistic case 
         !and one for the regular case
         if (esize==2) then  
@@ -441,9 +440,7 @@ module mode_create
 
                         if(all(node_in_bd)) then 
                             lat_ele_num = lat_ele_num+1
-                            do inod = 1, 8
+                            r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)
-                                r_lat(:,inod,lat_ele_num) = temp_nodes(:,1,inod) + adjustVar(:,inod)
-                            end do
 
                             !Now set all the lattice points contained within an element to false
                             do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))