aselimov/CACmb
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add ability for even numbered elements

Browse Source
This commit is contained in:
Alex Selimov 2022-03-09 14:06:09 -05:00
parent a712db3274
commit f9e7f68676
1 changed files with 22 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
src/elements.f90
View File

@ -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
@ -104,6 +106,14 @@ module elements
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, &
0.5_dp, 0.5_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