From f9e7f686768e1944234c55763e16f3f77e670d34 Mon Sep 17 00:00:00 2001
From: Alex Selimov <alex@alexselimov.xyz>
Date: Wed, 9 Mar 2022 14:06:09 -0500
Subject: [PATCH] Add ability for even numbered elements

---
 src/elements.f90 | 34 ++++++++++++++++++++++------------
 1 file changed, 22 insertions(+), 12 deletions(-)

diff --git a/src/elements.f90 b/src/elements.f90
index 655f4e5..8a14b86 100644
--- a/src/elements.f90
+++ b/src/elements.f90
@@ -2,6 +2,7 @@ module elements
     
     !This module contains the elements data structures, structures needed for building regions
     !and operations that are done on elements
+    use atoms
     use parameters
     use functions
     use subroutines
@@ -31,20 +32,21 @@ module elements
     !Mapping atom type to provided name
     character(len=2), dimension(10) :: type_to_name
     integer :: atom_types = 0
+    real(kind=dp) :: masses(10)
 
     !Variables for creating elements based on primitive cells
     real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3), bcc_cell(3,8), bcc_mat(3,3), bcc_inv(3,3), &
                      cube20(3,20)   
-    integer :: cubic_faces(4,6)
+    integer :: cubic_faces(4,6), oneonetwopairs(2,6)
 
     !Below are lattice type arrays which provide information on the general form of the elements.
     !We currently have a limit of 10 lattice types for simplicities sake but this can be easily increased.
     integer :: lattice_types = 0
     integer :: max_ng_node, ng_node(10) !Max number of nodes per element and number of nodes per element for each lattice type
     integer :: max_esize=0 !Maximum number of atoms per side of element
-    real(kind=dp) :: lapa(10), masses(10)
+    real(kind=dp) :: lapa(10)
 
-    !These variables contain information on the basis, for simplicities sake we limit 
+    !These variables contain informatione on the basis, for simplicities sake we limit 
     !the user to the definition of 10 lattice types with 10 basis atoms at each lattice point.
     !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 
@@ -103,6 +105,14 @@ module elements
         cubic_faces(:,4) = (/ 1, 4, 8, 5 /)
         cubic_faces(:,5) = (/ 4, 3, 7, 8 /)
         cubic_faces(:,6) = (/ 5, 6, 7, 8 /)
+        
+        !Now mark which node pairs represent the 112 directions. This is used for the dislocation analysis.
+        oneonetwopairs(:,1) = (/ 1, 3 /)
+        oneonetwopairs(:,2) = (/ 1, 6 /)
+        oneonetwopairs(:,3) = (/ 2, 7 /)
+        oneonetwopairs(:,4) = (/ 1, 8 /)
+        oneonetwopairs(:,5) = (/ 4, 7 /)
+        oneonetwopairs(:,6) = (/ 5, 7 /)
 
         !!Now initialize the fcc primitive cell
         fcc_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, &
@@ -345,9 +355,9 @@ module elements
 
     end subroutine add_atom
 
-    subroutine add_atom_type(type, inttype, force_new)
+    subroutine add_atom_type(mass, inttype, force_new)
         !This subroutine adds a new atom type to the module list of atom types
-        character(len=2), intent(in) :: type
+        real(kind=dp), intent(in) :: mass
         integer, intent(out) :: inttype
 
         logical, optional, intent(in) :: force_new
@@ -364,7 +374,7 @@ module elements
         exists = .false.
         if(.not.force) then 
             do i=1, 10
-                if(type == type_to_name(i)) then 
+                if(is_equal(mass, masses(i))) then 
                     exists = .true.
                     inttype = i
                     exit
@@ -378,7 +388,8 @@ module elements
                 print *, "Defined atom types are greater than 10 which is currently not supported."
                 stop 3
             end if
-            type_to_name(atom_types) = type
+            call atommassspecies(mass, type_to_name(atom_types))
+            masses(atom_types) = mass
             inttype = atom_types
         end if
         return
@@ -429,7 +440,7 @@ module elements
         real(kind=dp), dimension(10, max_basisnum*max_esize**3), optional,intent(out) :: data_interp !Interpolated atomic positions
 
         !Internal variables
-        integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp
+        integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp, adjust
         real(kind=dp) :: a_shape(max_ng_node)
         real(kind=dp) :: t, s, r
         
@@ -451,16 +462,15 @@ module elements
             lat_type_temp = lat_type
         end select
 
-
         select case(type)
         case('fcc','bcc')
             !Now loop over all the possible sites
             do it = 1, esize
-                t = (1.0_dp*(it-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
+                t=-1.0_dp +(it-1)*(2.0_dp/(esize-1))   
                 do is =1, esize
-                    s = (1.0_dp*(is-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
+                    s=-1.0_dp +(is-1)*(2.0_dp/(esize-1))   
                     do ir = 1, esize
-                        r = (1.0_dp*(ir-1) - (esize-1)/2)/(1.0_dp*(esize-1)/2)
+                        r=-1.0_dp +(ir-1)*(2.0_dp/(esize-1))   
                         call rhombshape(r,s,t,a_shape)
 
                         do ibasis = 1, bnum