Merge pull request #3 from aselimov/development

v0.4.1
master
aselimov 5 years ago committed by GitHub
commit 520babe6f6
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23

@ -1,2 +1,114 @@
# CAC_Model_Builder # CAC_Model_Builder
This is a tool for building models in CAC This is a tool for building models in CAC. Commands and usage options are below. This code is intended to follow the atomsk code fairly closely.
## Modes
The modes follow similarly to the modes you find when using atomsk. The modes will be listed below alongside their syntax and other usage instructions. As a note, if a mode is being used then it has to come first.
### Mode Create
```
cacmb --create name element_type lattice_parameter esize
```
Mode create has the following parameters:
`name` - User defined name that either defines the atom type or the lattice type if using the basis option
`element_type` - Specifies which element type to use, this dictates the crystal being build. Current acceptable options for element_type are:
* FCC - Uses the Rhombohedral primitive fcc unit cell as the finite element.
`lattice_parameter` - The lattice parameter for the crystal structure.
`esize` - Number of atoms per edge of the finite element. A value of 2 signifies full atomistic resolution and is the lowest number acceptable.
**Example**
```
cacmb --create Cu fcc 3.615 11
```
Creates a copper element with a lattice parameter of 3.615 with 11 atoms per side
#### Optional keywords
**Orient**
``` orient [hkl] [hkl] [hkl]
orient [hkl] [hkl] [hkl]
```
Default orientation is `[100] [010] [001]`. If this keyword is present then the user must provide the orientation matrix in form `[hkl] [hkl] [hkl]`.
*Example:* `orient [-112] [110] [-11-1]`
**Basis**
```
basis num atom_name x y z
```
Default basis has `atom_name = name` with position (0,0,0). If used then the `atom_name x y z` must be include `num` times.
*Example:* `basis 2 Mg 0 0 0 Mg 0.5 0.288675 0.81647`
**Duplicate**
```
duplicate numx numy numz
```
Default duplicate is `1 1 1`. This is used to replicate the element along each dimensions. This cannot be used if the keyword dimensions is included. By default jagged edges along boundaries are filled if duplicate is greater than `1 1 1`.
*Example:* `duplicate 10 10 10`
**Dimensions**
```
dimensions dimx dimy dimz
```
There is no default dimensions as duplicate is the default option. This command assigns a box with user-assigned dimensions and fills it with the desired element. By default atoms fill in the jagged edges at the boundaries if the dimensions command is included.
Example: `dimensions 100 100 100`
**ZigZag**
```
zigzag boolx booly boolz
```
Default zigzag is `f f f`. This command specifies whether a boundary should be left jagged (i.e. in essence not filled in). If `boolx` is `t` than the x dimension is left jagged and if it is `f` then the x dimension is filled.
*Example:* `zigzag t f t` gives a box with jagged edges in the x and z and filled edges in the y.
**Origin**
```
origin x y z
```
Default origin is `0 0 0`. This command just sets the origin for where the simulation cell starts building.
*Example:* `origin 10 0 1`
### Mode Convert
```
cacmb --convert infile outfile
```
This mode converts a file `infile` to a file of `outfile`. The extensions determine the conversion process.
### Mode Merge
```
cacmb --merge dim N infiles outfile
```
This mode merges multiple data files and creates one big simulation cell. The parameters are:
`N` - The number of files which are being read
`dim` - the dimension they are to be stacked along, can be either `x`, `y`, or `z`. If the argument `none` is passed then the cells are just overlaid. Future options will include a delete overlap command.

@ -0,0 +1,33 @@
FC=ifort
FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin
#FFLAGS=-c -mcmodel=large -Ofast
MODES=mode_create.o mode_merge.o mode_convert.o
OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES)
.SUFFIXES:
.SUFFIXES: .c .f .f90 .F90 .o
cacmb: $(OBJECTS)
$(FC) $(FFLAGS) $(OBJECTS) -o $@
.f90.o:
$(FC) $(FFLAGS) -c $<
.PHONY: clean
clean:
$(RM) cacmb *.o
testfuncs: testfuncs.o functions.o subroutines.o
$(FC) testfuncs.o functions.o subroutines.o elements.o -o $@
.PHONY: cleantest
cleantest:
$(RM) testfuncs testfuncs.o
$(OBJECTS) : parameters.o
atoms.o subroutines.o testfuncs.o : functions.o
main.o io.o build_subroutines.o: elements.o
call_mode.o : $(MODES)
$(MODES) io.o: atoms.o box.o
$(MODES) main.o : io.o
testfuncs.o elements.o mode_create.o: subroutines.o

File diff suppressed because it is too large Load Diff

@ -0,0 +1,75 @@
module box
!This module contains information on the properties of the current box.
use parameters
implicit none
real(kind=dp) :: box_bd(6) !Global box boundaries
!The subbox variables contain values for each subbox, being the boxes read in through some
!command. Currently only mode_merge will require sub_boxes, for mode_create it will always
!allocate to only 1 sub_box
integer :: sub_box_num = 0
integer, allocatable :: sub_box_array_bd(:,:,:)!Boundaries in the atom and element arrays for each sub_box
real(kind=dp), allocatable :: sub_box_ori(:,:,:)!Orientations for each of the subboxes
real(kind=dp), allocatable :: sub_box_bd(:,:)!Boundaries for each of the sub_boxes
public
contains
subroutine box_init
!Initialize some box functions
box_bd(:) = 0.0_dp
end subroutine box_init
subroutine alloc_sub_box(n)
!Allocate the sub_box variables
integer, intent(in) :: n
allocate(sub_box_ori(3,3,n), sub_box_bd(6,n), sub_box_array_bd(2,2,n))
end subroutine alloc_sub_box
subroutine grow_sub_box(n)
!Grows sub box arrays, this is only called when a new file is read in
integer, intent(in) :: n
integer, allocatable :: temp_array_bd(:,:,:), temp_file(:)
real(kind=dp), allocatable :: temp_ori(:,:,:), temp_bd(:,:)
!Allocate temporary arrays
allocate(temp_ori(3,3,sub_box_num+n),temp_bd(6,sub_box_num+n), &
temp_array_bd(2,2,sub_box_num+n), temp_file(sub_box_num+n))
!Move allocation for all sub_box_arrays
temp_ori(:,:,1:sub_box_num) = sub_box_ori
temp_ori(:,:,sub_box_num+1:) = 0.0_dp
call move_alloc(temp_ori, sub_box_ori)
temp_bd(:, 1:sub_box_num) = sub_box_bd
temp_bd(:, sub_box_num+1:) = 0.0_dp
call move_alloc(temp_bd, sub_box_bd)
temp_array_bd(:,:,1:sub_box_num) = sub_box_array_bd
temp_array_bd(:,:,sub_box_num+1:) = 0.0_dp
call move_alloc(temp_array_bd, sub_box_array_bd)
return
end subroutine grow_sub_box
subroutine grow_box(temp_box_bd)
!This function takes in a temporary box boundary and adjusts the overall box boundaries
!to include it
real(kind=dp), dimension(6), intent(in) :: temp_box_bd
integer :: i
do i = 1, 3
if(temp_box_bd(2*i-1) < box_bd(2*i-1)) box_bd(2*i-1) = temp_box_bd(2*i-1)
if(temp_box_bd(2*i) > box_bd(2*i)) box_bd(2*i) = temp_box_bd(2*i)
end do
return
end subroutine grow_box
end module box

@ -0,0 +1,29 @@
subroutine call_mode(arg_num,mode)
!This code is used to parse the command line argument for the mode information and calls the required
!mode module.
use mode_create
use mode_convert
use mode_merge
use parameters
implicit none
integer, intent(in) :: arg_num
character(len=100), intent(in) :: mode
select case(mode)
case('--create')
call create
case('--convert')
call convert
case('--merge')
call merge
case default
print *, "Mode ", trim(adjustl(mode)), " currently not accepted. Please check documentation for ", &
"accepted modes and rerun."
stop 3
end select
end subroutine call_mode

@ -0,0 +1,352 @@
module elements
!This module contains the elements data structures, structures needed for building regions
!and operations that are done on elements
use parameters
use subroutines
implicit none
!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
real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array
integer :: ele_num=0 !Number of elements
integer :: node_num=0 !Total number of nodes
!Data structure used to represent atoms
integer, allocatable :: type_atom(:)!atom type
real(kind =dp),allocatable :: r_atom(:,:) !atom position
integer :: atom_num=0 !Number of atoms
!Mapping atom type to provided name
character(len=2), dimension(10) :: type_to_name
integer :: atom_types = 0
!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)
integer :: cubic_faces(4,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
!These variables contain information 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
integer :: basis_type(10,10)
public
contains
subroutine lattice_init
!This subroutine just intializes variables needed for building the different finite
!element types
!First initialize the cubic cell
cubic_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, &
1.0_dp, 0.0_dp, 0.0_dp, &
1.0_dp, 1.0_dp, 0.0_dp, &
0.0_dp, 1.0_dp, 0.0_dp, &
0.0_dp, 0.0_dp, 1.0_dp, &
1.0_dp, 0.0_dp, 1.0_dp, &
1.0_dp, 1.0_dp, 1.0_dp, &
0.0_dp, 1.0_dp, 1.0_dp /), &
shape(fcc_cell))
!Now we create a list containing the list of vertices needed to describe the 6 cube faces
cubic_faces(:,1) = (/ 1, 4, 8, 5 /)
cubic_faces(:,2) = (/ 2, 3, 7, 6 /)
cubic_faces(:,3) = (/ 1, 2, 6, 5 /)
cubic_faces(:,4) = (/ 3, 4, 8, 7 /)
cubic_faces(:,5) = (/ 1, 2, 3, 4 /)
cubic_faces(:,6) = (/ 5, 6, 7, 8 /)
!!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, &
0.5_dp, 1.0_dp, 0.5_dp, &
0.0_dp, 0.5_dp, 0.5_dp, &
0.5_dp, 0.0_dp, 0.5_dp, &
1.0_dp, 0.5_dp, 0.5_dp, &
1.0_dp, 1.0_dp, 1.0_dp, &
0.5_dp, 0.5_dp, 1.0_dp /), &
shape(fcc_cell))
fcc_mat = reshape((/ 0.5_dp, 0.5_dp, 0.0_dp, &
0.0_dp, 0.5_dp, 0.5_dp, &
0.5_dp, 0.0_dp, 0.5_dp /), &
shape(fcc_mat))
call matrix_inverse(fcc_mat,3,fcc_inv)
max_basisnum = 0
basisnum(:) = 0
ng_node(:) = 0
end subroutine lattice_init
subroutine cell_init(lapa,esize,ele_type, orient_mat, cell_mat)
!This subroutine uses the user provided information to transform the finite element cell to the correct
!size, orientation, and dimensions using the esize, lattice parameter, element_type, and orientation provided
!by the user
integer, intent(in) :: esize
real(kind=dp), intent(in) :: lapa, orient_mat(3,3)
character(len=100), intent(in) :: ele_type
real(kind=dp), dimension(3,max_ng_node), intent(out) :: cell_mat
select case(trim(ele_type))
case('fcc')
cell_mat(:,1:8) = lapa * ((esize-1)*matmul(orient_mat, fcc_cell))
case default
print *, "Element type ", trim(ele_type), " currently not accepted"
stop
end select
end subroutine cell_init
subroutine alloc_ele_arrays(n,m)
!This subroutine used to provide initial allocation for the atom and element arrays
integer, intent(in) :: n, m !n-size of element arrays, m-size of atom arrays
integer :: allostat
!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), &
stat=allostat)
if(allostat > 0) then
print *, "Error allocating element arrays in elements.f90 because of: ", allostat
stop
end if
end if
if(m > 0) then
!Allocate atom arrays
allocate(type_atom(m), r_atom(3,m), stat=allostat)
if(allostat > 0) then
print *, "Error allocating atom arrays in elements.f90 because of: ", allostat
stop
end if
end if
end subroutine
subroutine grow_ele_arrays(n, m)
integer, intent(in) :: n, m
integer :: ele_size, atom_size, buffer_size
integer, allocatable :: temp_int(:)
real(kind=dp), allocatable :: temp_ele_real(:,:,:,:), temp_real(:,:)
character(len=100), allocatable :: char_temp(:)
!The default size we grow the
buffer_size = 1024
!Figure out the size of the atom and element arrays
ele_size = size(size_ele)
atom_size = size(type_atom)
!Check if we need to grow the ele_size, if so grow all the variables
if ( n+ele_num > size(size_ele)) then
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)
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)
allocate(char_temp(n+ele_num+buffer_size))
char_temp(1:ele_size) = type_ele
call move_alloc(char_temp, type_ele)
allocate(temp_ele_real(3, max_basisnum, max_ng_node, n+ele_num+buffer_size))
temp_ele_real(:,:,:,1:ele_size) = r_node
temp_ele_real(:,:,:,ele_size+1:) = 0.0_dp
call move_alloc(temp_ele_real, r_node)
end if
!Now grow atom arrays if needed
if (m+atom_num > atom_size) then
allocate(temp_int(m+atom_num+buffer_size))
temp_int(1:atom_size) = type_atom
temp_int(atom_size+1:) = 0
call move_alloc(temp_int, type_atom)
allocate(temp_real(3,m+atom_num+buffer_size))
temp_real(:,1:atom_size) = r_atom
temp_real(:, atom_size+1:) = 0.0_dp
call move_alloc(temp_real, r_atom)
end if
end subroutine
subroutine add_element(type, size, lat, r)
!Subroutine which adds an element to the element arrays
integer, intent(in) :: size, lat
character(len=100), intent(in) :: type
real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node)
ele_num = ele_num + 1
type_ele(ele_num) = type
size_ele(ele_num) = size
lat_ele(ele_num) = lat
r_node(:,:,:,ele_num) = r(:,:,:)
node_num = node_num + ng_node(lat)
end subroutine add_element
subroutine add_atom(type, r)
!Subroutine which adds an atom to the atom arrays
integer, intent(in) :: type
real(kind=dp), intent(in), dimension(3) :: r
atom_num = atom_num+1
type_atom(atom_num) = type
r_atom(:,atom_num) = r(:)
end subroutine add_atom
subroutine add_atom_type(type, inttype)
!This subroutine adds a new atom type to the module list of atom types
character(len=2), intent(in) :: type
integer, intent(out) :: inttype
integer :: i
logical :: exists
exists = .false.
do i=1, 10
if(type == type_to_name(i)) exists = .true.
inttype = i
end do
if (exists.eqv..false.) then
atom_types = atom_types+1
if(atom_types > 10) then
print *, "Defined atom types are greater than 10 which is currently not supported."
stop 3
end if
type_to_name(atom_types) = type
inttype = atom_types
end if
return
end subroutine add_atom_type
subroutine def_ng_node(n, element_types)
!This subroutine defines the maximum node number among n element types
integer, intent(in) :: n !Number of element types
character(len=100), dimension(n) :: element_types !Array of element type strings
integer :: i
max_ng_node = 0
do i=1, n
select case(trim(adjustl(element_types(i))))
case('fcc')
ng_node(i) = 8
end select
if(ng_node(i) > max_ng_node) max_ng_node = ng_node(i)
end do
end subroutine def_ng_node
subroutine set_max_esize
!This subroutine sets the maximum esize
max_esize=maxval(size_ele)
end subroutine
subroutine interpolate_atoms(type, esize, lat_type, r_in, type_interp, r_interp)
!This subroutine returns the interpolated atoms from the elements.
!Arguments
character(len=100), intent(in) :: type !The type of element that it is
integer, intent(in) :: esize !The number of atoms per side
integer, intent(in) :: lat_type !The integer lattice type of the element
real(kind=dp), dimension(3,max_basisnum, max_ng_node), intent(in) :: r_in !Nodal positions
integer, dimension(max_basisnum*max_esize**3), intent(out) :: type_interp !Interpolated atomic positions
real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: r_interp !Interpolated atomic positions
!Internal variables
integer :: i, it, is, ir, ibasis, inod, ia, bnum, lat_type_temp
real(kind=dp), allocatable :: a_shape(:)
real(kind=dp) :: t, s, r
!Initialize some variables
r_interp(:,:) = 0.0_dp
type_interp(:) = 0
ia = 0
!Define bnum based on the input lattice type. If lat_type=0 then we are interpolating lattice points which means
!the basis is 0,0,0, and the type doesn't matter
select case(lat_type)
case(0)
bnum=1
lat_type_temp = 1
case default
bnum = basisnum(lat_type)
lat_type_temp = lat_type
end select
select case(trim(adjustl(type)))
case('fcc')
allocate(a_shape(8))
!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)
do is =1, esize
s = (1.0_dp*(is-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
do ir = 1, esize
r = (1.0_dp*(ir-1) - (esize-1)/2)/(1.0_dp*(esize-1)/2)
call rhombshape(r,s,t,a_shape)
do ibasis = 1, bnum
ia = ia + 1
do inod = 1, 8
type_interp(ia) = basis_type(ibasis,lat_type_temp)
r_interp(:,ia) = r_interp(:,ia) + a_shape(inod) * r_in(:,ibasis,inod)
end do
end do
end do
end do
end do
end select
if (ia /= esize**3) then
print *, "Incorrect interpolation"
stop 3
end if
return
end subroutine interpolate_atoms
subroutine rhombshape(r,s,t, shape_fun)
!Shape function for rhombohedral elements
real(kind=8), intent(in) :: r, s, t
real(kind=8), intent(out) :: shape_fun(8)
shape_fun(1) = (1.0-r)*(1.0-s)*(1.0-t)/8.0
shape_fun(2) = (1.0+r)*(1.0-s)*(1.0-t)/8.0
shape_fun(3) = (1.0+r)*(1.0+s)*(1.0-t)/8.0
shape_fun(4) = (1.0-r)*(1.0+s)*(1.0-t)/8.0
shape_fun(5) = (1.0-r)*(1.0-s)*(1.0+t)/8.0
shape_fun(6) = (1.0+r)*(1.0-s)*(1.0+t)/8.0
shape_fun(7) = (1.0+r)*(1.0+s)*(1.0+t)/8.0
shape_fun(8) = (1.0-r)*(1.0+s)*(1.0+t)/8.0
return
end subroutine rhombshape
end module elements

@ -0,0 +1,235 @@
module functions
use parameters
implicit none
public
contains
! Functions below this comment are taken from the functions module of atomsk
!********************************************************
! STRUPCASE
! This function reads a string of any length
! and capitalizes all letters.
!********************************************************
FUNCTION StrUpCase (input_string) RESULT (UC_string)
!
IMPLICIT NONE
CHARACTER(*),PARAMETER:: lower_case = 'abcdefghijklmnopqrstuvwxyz'
CHARACTER(*),PARAMETER:: upper_case = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
CHARACTER(*),INTENT(IN):: input_string
CHARACTER(LEN(Input_String)):: UC_string !Upper-Case string that is produced
INTEGER:: i, n
!
IF(LEN(input_string)==0) RETURN
UC_string = input_string
! Loop over string elements
DO i=1,LEN(UC_string)
!Find location of letter in lower case constant string
n = INDEX( lower_case, UC_string(i:i) )
!If current substring is a lower case letter, make it upper case
IF(n>0) THEN
UC_string(i:i) = upper_case(n:n)
ENDIF
END DO
!
END FUNCTION StrUpCase
!********************************************************
! STRDNCASE
! This function reads a string of any length
! and transforms all letters to lower case.
!********************************************************
FUNCTION StrDnCase (input_string) RESULT (lc_string)
!
IMPLICIT NONE
CHARACTER(*),PARAMETER:: lower_case = 'abcdefghijklmnopqrstuvwxyz'
CHARACTER(*),PARAMETER:: upper_case = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
CHARACTER(*),INTENT(IN):: input_string
CHARACTER(LEN(Input_String)):: lc_string !Lower-Case string that is produced
INTEGER:: i, n
!
IF(LEN(input_string)==0) RETURN
lc_string = input_string
! Loop over string elements
DO i=1,LEN(lc_string)
!Find location of letter in upper case constant string
n = INDEX( upper_case, lc_string(i:i) )
!If current substring is an upper case letter, make it lower case
IF(n>0) THEN
lc_string(i:i) = lower_case(n:n)
ENDIF
END DO
!
END FUNCTION StrDnCase
pure function matrix_normal(a, n)
integer :: i
integer, intent(in) :: n
real(kind = dp), dimension(n) :: v
real(kind = dp), dimension(n, n),intent(in) :: a
real(kind = dp), dimension(n,n) :: matrix_normal
matrix_normal(:, :) = a(:, :)
do i = 1, n
v(:) = a(i,:)
matrix_normal(i, :) = v(:) / norm2(v)
end do
return
end function matrix_normal
pure function cross_product(a, b)
!Function which finds the cross product of two vectors
real(kind = dp), dimension(3), intent(in) :: a, b
real(kind = dp), dimension(3) :: cross_product
cross_product(1) = a(2) * b(3) - a(3) * b(2)
cross_product(2) = a(3) * b(1) - a(1) * b(3)
cross_product(3) = a(1) * b(2) - a(2) * b(1)
return
end function cross_product
pure function identity_mat(n)
!Returns the nxn identity matrix
integer :: i
integer, intent(in) :: n
real(kind = dp), dimension(n, n) :: identity_mat
identity_mat(:, :) = 0.0_dp
do i = 1, n
identity_mat(i, i) = 1.0_dp
end do
return
end function identity_mat
pure function triple_product(a, b, c)
!triple product between three 3*1 vectors
real(kind = dp) :: triple_product
real(kind = dp), dimension(3), intent(in) :: a, b, c
triple_product = dot_product(a, cross_product(b, c))
return
end function triple_product
function in_bd_lat(v, box_faces, box_norms)
!This function returns whether the point is within the transformed box boundaries. The transformed
!box being the transformed simulation cell in the lattice basis
!Input/output variables
real(kind=dp), dimension(3), intent(in) :: v !integer lattice position
real(kind=dp), dimension(3,6), intent(in) :: box_faces !Centroid of all the box faces
real(kind=dp), dimension(3,6), intent(in) :: box_norms !Box face normals
logical :: in_bd_lat
!Other variables
integer :: i
real(kind=dp) :: pt_to_face(3)
in_bd_lat = .true.
!Check if point is in box bounds, this works by comparing the dot product of the face normal and the
!vector from the point to the face. If the dot product is greater than 0 then the point is behind the face
!if it is equal to zero then the point is on the face, if is less than 0 the the point is in front of the face.
do i = 1, 6
pt_to_face(:) = box_faces(:, i) - v
if(dot_product(pt_to_face, box_norms(:,i)) <= 0) then
in_bd_lat = .false.
exit
end if
end do
return
end function in_bd_lat
function in_block_bd(v, box_bd)
!This function returns whether a point is within a block in 3d
!Input/output
real(kind=dp), dimension(3), intent(in) :: v
real(kind=dp), dimension(6), intent(in) :: box_bd
logical :: in_block_bd
!Other variables
integer :: i
in_block_bd = .true.
do i =1 ,3
!Check upper bound
if(v(i) > (box_bd(2*i)+10.0_dp**(-10)) ) then
in_block_bd =.false.
exit
!Check lower bound
else if (v(i) < (box_bd(2*i-1)-10.0_dp**(-10))) then
in_block_bd = .false.
exit
end if
end do
end function in_block_bd
function lcm(a,b)
!This function returns the smallest least common multiple of two numbers
real(kind=dp), intent(in) :: a, b
real(kind=dp) :: lcm
integer :: aint, bint, gcd, remainder, placeholder
!Cast the vector positions to ints. There will be some error associated with this calculation
aint = a*10**2
bint = b*10**2
!Calculate greated common divisor
gcd = aint
placeholder = bint
do while(placeholder /= 0)
remainder = modulo(gcd, placeholder)
gcd = placeholder
placeholder=remainder
end do
lcm = real((aint*bint),dp)/(real(gcd,dp))* 10.0_dp**(-2.0_dp)
end function lcm
function is_neighbor(rl, rk, r_in, r_out)
!This function checks to see if two atoms are within a shell with an inner radius r_in and outer radius
!r_out
real(kind=dp), intent(in) :: r_in, r_out
real(kind=dp), dimension(3), intent(in) :: rl, rk
logical :: is_neighbor
!Internal variable
real(kind=dp) :: rlk
rlk = norm2(rk - rl)
is_neighbor=.true.
if((rlk>r_out).or.(rlk < r_in)) is_neighbor = .false.
return
end function is_neighbor
function is_equal(A, B)
!This function checks if too numbers are equal within a tolerance
real(kind=dp), intent(in) :: A, B
logical :: is_equal
if((A>(B - 10.0_dp**(-10))).and.(A < (B+10.0_dp**(-10)))) then
is_equal = .true.
else
is_equal = .false.
end if
return
end function is_equal
end module functions

@ -0,0 +1,472 @@
module io
use elements
use parameters
use atoms
use box
implicit none
integer :: outfilenum = 0, infilenum = 0
character(len=100) :: outfiles(10), infiles(10)
public
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Subroutines for writing out data files !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine get_out_file(filename)
implicit none
character(len=100), intent(in) :: filename
character(len=100) :: temp_outfile
character(len=1) :: overwrite
logical :: file_exists
!If no filename is provided then this function is called with none and prompts user input
if (filename=='none') then
print *, "Please specify a filename or extension to output to:"
read(*,*) temp_outfile
else
temp_outfile = filename
end if
!Infinite loop which only exists if user provides valid filetype
overwrite = 'r'
do while(.true.)
!Check to see if file exists, if it does then ask user if they would like to overwrite the file
inquire(file=trim(temp_outfile), exist=file_exists)
if (file_exists) then
if (overwrite == 'r') print *, "File ", trim(temp_outfile), " already exists. Would you like to overwrite? (Y/N)"
read(*,*) overwrite
if((scan(overwrite, "n") > 0).or.(scan(overwrite, "N") > 0)) then
print *, "Please specify a new filename with extension:"
read(*,*) temp_outfile
else if((scan(overwrite, "y") > 0).or.(scan(overwrite, "Y") > 0)) then
continue
else
print *, "Please pick either y or n"
read(*,*) overwrite
end if
end if
if (scan(temp_outfile,'.',.true.) == 0) then
print *, "No extension included on filename, please type a full filename that includes an extension."
read(*,*) temp_outfile
cycle
end if
select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:))
case('xyz', 'lmp', 'vtk', 'mb')
outfilenum=outfilenum+1
outfiles(outfilenum) = temp_outfile
exit
case default
print *, "File type: ", trim(temp_outfile(scan(temp_outfile,'.',.true.):)), "not currently accepted. ", &
"please input a filename with extension from following list: xyz, lmp, vtk."
read(*,*) temp_outfile
end select
end do
end subroutine get_out_file
subroutine write_out
!This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine
integer :: i
!Find max esize which will be needed later
call set_max_esize
do i = 1, outfilenum
!Pull out the extension of the file and call the correct write subroutine
select case(trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))))
case('xyz')
call write_xyz(outfiles(i))
case('lmp')
call write_lmp(outfiles(i))
case('vtk')
call write_vtk(outfiles(i))
case('mb')
call write_mb(outfiles(i))
case default
print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
" is not accepted for writing. Please select from: xyz and try again"
stop
end select
end do
end subroutine write_out
subroutine write_xyz(file)
!This is the simplest visualization subroutine, it writes out all nodal positions and atom positions to an xyz file
character(len=100), intent(in) :: file
integer :: node_num, i, inod, ibasis
open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
!Calculate total node number
node_num=0
do i = 1, ele_num
node_num = node_num + basisnum(lat_ele(i))*ng_node(lat_ele(i))
end do
!Write total number of atoms + elements
write(11, '(i16)') node_num+atom_num
!Write comment line
write(11, '(a)') "#Node + atom file created using cacmb"
!Write nodal positions
do i = 1, ele_num
do inod = 1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
write(11, '(a, 3f23.15)') basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i)
end do
end do
end do
!Write atom positions
do i = 1, atom_num
write(11, '(a, 3f23.15)') type_atom(i), r_atom(:,i)
end do
!Finish writing
close(11)
end subroutine write_xyz
subroutine write_lmp(file)
!This subroutine writes out a .lmp style dump file
character(len=100), intent(in) :: file
integer :: write_num, i, iatom, type_interp(max_basisnum*max_esize**3)
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), mass
open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
!Comment line
write(11, '(a)') '# lmp file made with cacmb'
write(11, '(a)')
!Calculate total atom number
write_num = atom_num
do i = 1,ele_num
if(type_ele(i) == 'fcc') write_num = write_num + size_ele(i)**3
end do
!Write total number of atoms + elements
write(11, '(i16, a)') write_num, ' atoms'
!Write number of atom types
write(11, '(i16, a)') atom_types, ' atom types'
write(11,'(a)') ' '
!Write box bd
write(11, '(2f23.15, a)') box_bd(1:2), ' xlo xhi'
write(11, '(2f23.15, a)') box_bd(3:4), ' ylo yhi'
write(11, '(2f23.15, a)') box_bd(5:6), ' zlo zhi'
!Masses
write(11, '(a)') 'Masses'
write(11, '(a)') ' '
do i =1, atom_types
call atommass(type_to_name(i),mass)
write(11, '(i16, f23.15)') i, mass
end do
write(11, '(a)') ' '
!Write atom positions
write(11, '(a)') 'Atoms'
write(11, '(a)') ' '
do i = 1, atom_num
write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i)
end do
!Write refined element atomic positions
do i = 1, ele_num
call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp)
select case(trim(adjustl(type_ele(i))))
case('fcc')
do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
write(11, '(2i16, 3f23.15)') atom_num+iatom, type_interp(iatom), r_interp(:,iatom)
end do
end select
end do
end subroutine write_lmp
subroutine write_vtk(file)
!This subroutine writes out a vtk style dump file
integer :: i, j, inod, ibasis
character(len=100), intent(in) :: file
1 format('# vtk DataFile Version 4.0.1', / &
'CAC output -- cg', / &
'ASCII')
11 format('# vtk DataFile Version 4.0.1', / &
'CACmb output -- atoms', / &
'ASCII')
2 format('DATASET UNSTRUCTURED_GRID')
3 format('POINTS', i16, ' float')
4 format(/'CELLS', 2i16)
5 format(/'CELL_TYPES', i16)
12 format(/'CELL_DATA', i16)
16 format(/'POINT_DATA', i16)
17 format('SCALARS weight float', / &
'LOOKUP_TABLE default')
18 format('SCALARS atom_type float', / &
'LOOKUP_TABLE default')
20 format('SCALARS lattice_type 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')
write(11, 11)
write(11, 2)
write(11, 3) atom_num
do i = 1, atom_num
write(11, '(3f23.15)') r_atom(:,i)
end do
write(11,4) atom_num, atom_num*2
do i = 1, atom_num
write(11, '(2i16)') 1, i-1
end do
write(11, 5) atom_num
do i = 1, atom_num
write(11, '(i16)') 1
end do
write(11, 16) atom_num
write(11, 18)
do i = 1, atom_num
write(11, '(i16)') type_atom(i)
end do
close(11)
open(unit=11, file='cg_'//trim(adjustl(file)), action='write', status='replace',position='rewind')
write(11,1)
write(11,2)
write(11,3) node_num
do i = 1, ele_num
do inod=1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
write(11, '(3f23.1)') sum(r_node(:,:,inod,i),2)/basisnum(lat_ele(i))
end do
end do
end do
write(11, 4) ele_num, ele_num + node_num
do i =1, ele_num
write(11, '(9i16)') ng_node(lat_ele(i)), (j, j = (i-1)*ng_node(lat_ele(i)), i*ng_node(lat_ele(i))-1)
end do
write(11,5) ele_num
do i = 1, ele_num
if(trim(adjustl(type_ele(i))) == 'fcc') write(11, '(i16)') 12
end do
write(11,12) ele_num
write(11,20)
do i = 1, ele_num
write(11, '(i16)') lat_ele(i)
end do
close(11)
end subroutine
subroutine write_mb(file)
!This subroutine writes the cacmb formatted file which provides necessary information for building models
character(len=100), intent(in) :: file
integer :: i, j, k, inod, ibasis
!Open the .mb file for writing
open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
!First write the box boundary information
!Write the global box boundaries
write(11,*) box_bd(:)
!Write the number of sub_boxes in the system
write(11,*) sub_box_num
!For every subbox write the orientation, sub box boundary, and sub_box_array_bds
do i = 1, sub_box_num
write(11,*) sub_box_ori(:,:,i)
write(11,*) sub_box_bd(:,i)
write(11,*) ((sub_box_array_bd(j,k,i), j = 1, 2), k = 1, 2)
end do
!Write the number of atom types in the current model and all of their names
write(11,*) atom_types, (type_to_name(i), i=1, atom_types)
!Write the number of lattice_types, basisnum and number of nodes for each lattice type
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 write the numbers of elements and atoms
write(11,*) atom_num, ele_num
!Write out atoms first
do i = 1, atom_num
write(11,*) i, type_atom(i), r_atom(:,i)
end do
!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)
do inod = 1, ng_node(lat_ele(i))
do ibasis =1, basisnum(lat_ele(i))
write(11,*) inod, ibasis, r_node(:, ibasis, inod, i)
end do
end do
end do
end subroutine write_mb
!!!!!!!!!!!!! Below are subroutines for reading files !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine get_in_file(filename)
implicit none
character(len=100), intent(in) :: filename
character(len=100) :: temp_infile
logical :: file_exists
!If no filename is provided then this function is called with none and prompts user input
if (filename=='none') then
print *, "Please specify a filename or extension to output to:"
read(*,*) temp_infile
else
temp_infile = filename
end if
!Infinite loop which only exists if user provides valid filetype
do while(.true.)
!Check to see if file exists, if it does then ask user if they would like to overwrite the file
inquire(file=trim(temp_infile), exist=file_exists)
if (.not.file_exists) then
print *, "The file ", trim(adjustl(filename)), " does not exist. Please input a filename that exists"
read(*,*) temp_infile
cycle
end if
select case(temp_infile(scan(temp_infile,'.',.true.)+1:))
case('xyz', 'lmp', 'vtk', 'mb')
infilenum=infilenum+1
infiles(infilenum) = temp_infile
exit
case default
print *, "File type: ", trim(temp_infile(scan(temp_infile,'.',.true.):)), "not currently accepted. ", &
"please input a filename with extension from following list: mb."
read(*,*) temp_infile
end select
end do
end subroutine get_in_file
subroutine read_in(i, displace, temp_box_bd)
!This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine
integer, intent(in) :: i
real(kind=dp), dimension(3), intent(in) :: displace
real(kind=dp), dimension(6), intent(out) :: temp_box_bd
!Pull out the extension of the file and call the correct write subroutine
select case(trim(adjustl(infiles(i)(scan(infiles(i),'.',.true.)+1:))))
case('mb')
call read_mb(infiles(i), displace, temp_box_bd)
case default
print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
" is not accepted for writing. Please select from: mb and try again"
stop
end select
end subroutine read_in
subroutine read_mb(file, displace, temp_box_bd)
!This subroutine reads in an mb file for operation
character(len=100), intent(in) :: file
real(kind=dp), dimension(3), intent(in) :: displace
real(kind = dp), dimension(6), intent(out) :: temp_box_bd
integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles
character(len=100) :: etype
real(kind=dp) :: r(3), newdisplace(3)
real(kind=dp), allocatable :: r_innode(:,:,:)
!First open the file
open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
!Read in the box boundary and grow the current active box bd
read(11, *) temp_box_bd(:)
do i = 1, 3
newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
temp_box_bd(2*i-1) = temp_box_bd(2*i-1) + newdisplace(i)
temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i)
end do
!Read in the number of sub_boxes and allocate the variables
read(11, *) n
if (sub_box_num == 0) then
call alloc_sub_box(n)
else
call grow_sub_box(n)
end if
!Read in subbox orientations and boundaries
do i = 1, n
!Read in orientation with column major ordering
read(11,*) ((sub_box_ori(j, k, sub_box_num+i), j = 1, 3), k = 1, 3)
!Read in subbox boundaries
read(11,*) sub_box_bd(:,sub_box_num+i)
sub_box_bd(:,sub_box_num+i) = sub_box_bd(:, sub_box_num+i) + displace(:)
!Read in sub_box_array_bd
read(11,*) ((sub_box_ori(j, k, sub_box_num+i), j = 1, 2), k = 1, 2)
end do
sub_box_num = sub_box_num + n
!Read in the number of atom types and all their names
read(11, *) atom_types, (type_to_name(i), i = 1, atom_types)
!Read the number of lattice types, basisnum, and number of nodes for each lattice type
read(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types)
!Define max_ng_node and max_basis_num
max_basisnum = maxval(basisnum)
max_ng_node = maxval(ng_node)
!Read the basis atom types for every lattice
read(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = 1, 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)
allocate(r_innode(3,max_basisnum, max_ng_node))
!Read the atoms
do i = 1, in_atoms
read(11,*) j, type, r(:)
call add_atom(type, r+newdisplace)
end do
!Read the elements
do i = 1, in_eles
read(11, *) n, type, size, etype
do inod = 1, ng_node(type)
do ibasis =1, basisnum(type)
read(11,*) j, k, r_innode(:, ibasis, inod)
r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace
end do
end do
call add_element(etype, size, type, r_innode)
end do
!Close the file being read
close(11)
end subroutine read_mb
end module io

@ -0,0 +1,41 @@
program main
!**************************** CACmb *******************************
!* CAC model building toolkit *
! ____________ *
! / / *
! / / *
! /___________/ *
! _|_ _|_ _|____________ *
! / / *
! / / *
! /___________/ *
! *
!*******************************************************************
use parameters
use elements
use io
integer :: arg_num
character(len=100) :: mode
!Call initialization functions
call lattice_init
call box_init
! Command line parsing
arg_num = command_argument_count()
!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, mode)
mode = trim(adjustl(mode))
if (mode(1:2) == '--') then
call call_mode(arg_num, mode)
end if
!Finish by writing the files
call write_out
end program main

@ -0,0 +1,28 @@
module mode_convert
use parameters
use box
use elements
use io
public
contains
subroutine convert
!This subroutine converts a single input file from one format to another
character(len=100) :: infile, outfile
real(kind = dp) :: temp_box_bd(6)
!We have to allocate the element and atom arrays with a size of 1 for the read in code to work
call alloc_ele_arrays(1,1)
!First read in the file
call get_command_argument(2, infile)
call get_in_file(infile)
call read_in(1, (/0.0_dp,0.0_dp,0.0_dp/), temp_box_bd)
call grow_box(temp_box_bd)
!Now get the outfile, writing is done after all the codes complete
call get_command_argument(3, outfile)
call get_out_file(outfile)
end subroutine convert
end module mode_convert

@ -0,0 +1,499 @@
module mode_create
!This mode is intended for creating element/atom regions and writing them to specific files.
use parameters
use atoms
use io
use subroutines
use elements
use box
implicit none
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)
integer :: esize, duplicate(3), 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
real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:)
public
contains
subroutine create()
! Main subroutine which controls execution
character(len=100) :: textholder
integer :: i, ibasis, inod
real(kind=dp), allocatable :: r_node_temp(:,:,:)
!Initialize default parameters
orient = reshape((/ 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp /), shape(orient))
cell_mat(:,:)=0.0_dp
name =''
element_type = ''
esize=0
lattice_parameter=0.0_dp
duplicate(:) = 0
box_len(:) = 0.0_dp
dup_flag = .false.
dim_flag = .false.
basisnum = 0
lat_ele_num = 0
lat_atom_num = 0
!First we parse the command
call parse_command()
! Before building do a check on the file
if (outfilenum == 0) then
textholder = 'none'
call get_out_file(textholder)
end if
!Now we setup the unit element and call other init subroutines
call def_ng_node(1, element_type)
allocate(r_node_temp(3,max_basisnum,max_ng_node))
if(dup_flag) then
!We initialize the cell with a lattice_parameter of 1 because we will add the lattice parameter later
call cell_init(1.0_dp, esize, element_type, orient, cell_mat)
do i = 1, 8
box_vert(:,i) = duplicate(:)*esize*lattice_space(:)*cubic_cell(:,i) + origin(:)
end do
call matrix_inverse(orient,3,orient_inv)
!Now get the rotated box vertex positions in lattice space. Should be integer units
box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
!Find the new maxlen
maxbd = maxval(matmul(orient,matmul(fcc_mat,box_lat_vert)),2)
do i = 1, 3
box_bd(2*i) = maxval(box_vert(i,:)) - 0.25_dp*lattice_space(i)
box_bd(2*i-1) = origin(i)-0.25_dp*lattice_space(i)
end do
!and then call the build function with the correct transformation matrix
select case(trim(adjustl(element_type)))
case('fcc')
call build_with_rhomb(box_lat_vert, fcc_mat)
case default
print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", &
"element type"
stop 3
end select
!Now that it is multiply by the lattice parameter
box_bd = box_bd*lattice_parameter
else if(dim_flag) then
continue
else
call cell_init(lattice_parameter, esize, element_type, orient, cell_mat)
!If the user doesn't pass any build instructions than we just put the cell mat into the element_array
call alloc_ele_arrays(1,0)
!Add the basis atoms to the unit cell
do inod = 1, max_ng_node
do ibasis = 1, basisnum(1)
r_node_temp(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis) + origin(:)
end do
end do
do i = 1,3
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)
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
if(dup_flag.or.dim_flag) then
!Allocate variables
call alloc_ele_arrays(lat_ele_num, lat_atom_num*basisnum(1))
if(lat_atom_num > 0) then
do i = 1, lat_atom_num
do ibasis = 1, basisnum(1)
call add_atom(basis_type(ibasis, 1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis))
end do
end do
deallocate(r_atom_lat)
end if
if(lat_ele_num > 0) then
do i = 1, lat_ele_num
do inod= 1, ng_node(1)
do ibasis = 1, basisnum(1)
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)
end do
end if
end if
!The last thing we do is setup the sub_box_boundaries
call alloc_sub_box(1)
sub_box_num = 1
sub_box_ori(:,:,1) = orient
sub_box_bd(:,1) = box_bd
sub_box_array_bd(1,:,1) = 1
sub_box_array_bd(2,1,1) = atom_num
sub_box_array_bd(2,2,1) = ele_num
end subroutine create
!This subroutine parses the command and pulls out information needed for mode_create
subroutine parse_command()
integer :: arg_pos, ori_pos, i, j, arglen, stat
character(len=100) :: textholder
character(len=8) :: orient_string
!Pull out all required positional arguments
call get_command_argument(2, name, arglen)
if(arglen==0) STOP "Name is missing in mode create"
call get_command_argument(3,element_type, arglen)
if(arglen==0) STOP "Element_type is missing in mode create"
call get_command_argument(4,textholder, arglen)
if(arglen==0) STOP "Lattice Parameter is missing in mode create"
read(textholder, *, iostat=stat) lattice_parameter
if(stat > 0) STOP "Error reading lattice parameter"
call get_command_argument(5, textholder, arglen)
if(arglen==0) STOP "Esize missing in mode create"
read(textholder, *, iostat=stat) esize
if(stat > 0) STOP "Error reading esize"
arg_pos = 6
!Check for optional keywords
do while(.true.)
if(arg_pos > command_argument_count()) exit
!Pull out the next argument which should either be a keyword or an option
call get_command_argument(arg_pos, textholder)
textholder=adjustl(textholder)
arg_pos=arg_pos+1
!Choose what to based on what the option string is
select case(trim(textholder))
!If orient command is passed extract the orientation to numeric array format
case('orient')
do i = 1, 3
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
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
!If the duplicate command is passed then we extract the information on the new bounds.
case('duplicate')
dup_flag = .true.
do i = 1, 3
call get_command_argument(arg_pos, textholder)
read(textholder, *) duplicate(i)
arg_pos = arg_pos + 1
end do
case('origin')
do i = 1, 3
call get_command_argument(arg_pos, textholder)
read(textholder, *) origin(i)
arg_pos = arg_pos + 1
end do
!If a filetype is passed then we add name.ext to the outfiles list
case('xyz')
textholder = trim(adjustl(name)) //'.xyz'
call get_out_file(textholder)
case default
!Check to see if it is an option command, if so then mode_create must be finished
if(textholder(1:1) == '-') then
exit
!Check to see if a filename was passed
elseif(scan(textholder,'.',.true.) > 0) then
call get_out_file(textholder)
end if
end select
end do
!Calculate the lattice periodicity length in lattice units
do i = 1, 3
lattice_space(i) = norm2(orient(i,:))
end do
!Check special periodicity relations
select case(trim(adjustl(element_type)))
case('fcc')
do i = 1,3
!Check if one of the directions is 110
if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(orient(i,3),0.0_dp))).or.&
(is_equal(abs(orient(i,2)), abs(orient(i,3))).and.(is_equal(orient(i,1),0.0_dp))).or.&
(is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(orient(i,2),0.0_dp)))) then
lattice_space(i) = 0.5_dp * lattice_space(i)
!Check if one direction is 112
else if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(abs(orient(i,3)),2.0_dp*abs(orient(i,1))))).or.&
(is_equal(abs(orient(i,2)), abs(orient(i,3))).and.(is_equal(abs(orient(i,1)),2.0_dp*abs(orient(i,2))))).or.&
(is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(abs(orient(i,2)),2.0_dp*abs(orient(i,3))))))&
then
lattice_space(i) = 0.5_dp * lattice_space(i)
end if
end do
end select
!Now normalize the orientation matrix
orient = matrix_normal(orient,3)
!Set lattice_num to 1
lattice_types = 1
!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
max_basisnum = 1
basisnum(1) = 1
call add_atom_type(name, basis_type(1,1)) !If basis command not defined then we use name as the atom_name
basis_pos(:,1) = 0.0_dp
end if
!
end subroutine
subroutine build_with_rhomb(box_in_lat, transform_matrix)
!This subroutine returns all the lattice points in the box in r_lat
!Inputs
integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space
real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space
!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)
real(kind=dp), allocatable :: resize_lat_array(:,:)
logical, allocatable :: lat_points(:,:,:)
logical :: node_in_bd(8)
!Do some value initialization
max_esize = esize
!First find the bounding lattice points (min and max points for the box in each dimension)
numlatpoints = 1
do i = 1, 3
bd_in_lat(2*i-1) = minval(box_in_lat(i,:))
bd_in_lat(2*i) = maxval(box_in_lat(i,:))
numlatpoints = numlatpoints*(bd_in_lat(2*i)-bd_in_lat(2*i-1))
end do
!Allocate the correct lat variable
select case(esize)
!Atomistics
case(2)
allocate(r_atom_lat(3,numlatpoints))
case default
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
!atomistics
do iz = bd_in_lat(5)-5, bd_in_lat(6)+5
do iy = bd_in_lat(3)-5, bd_in_lat(4)+5
do ix = bd_in_lat(1)-5, bd_in_lat(2)+5
v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /)
!Transform point back to real space for easier checking
v = matmul(orient, matmul(transform_matrix,v))
!If within the boundaries
if(in_block_bd(v, box_bd)) then
lat_atom_num = lat_atom_num + 1
r_atom_lat(:,lat_atom_num) = v
end if
end do
end do
end do
else
!If we are working with elements we have to use more complex code
!Initialize finite element
ele(:,:) = (esize-1) * cubic_cell(:,:)
!Make a 3 dimensional array of lattice points. This array is indexed by the integer lattice position.
!A value of true means that the coordinate is a lattice point which is within the box_bd
allocate(lat_points(bd_in_lat(2)-bd_in_lat(1)+10,bd_in_lat(4)-bd_in_lat(3)+10,bd_in_lat(6)-bd_in_lat(5)+10))
lat_points(:,:,:) = .false.
do iz = bd_in_lat(5)-5, bd_in_lat(6)+5
do iy = bd_in_lat(3)-5, bd_in_lat(4)+5
do ix = bd_in_lat(1)-5, bd_in_lat(2)+5
v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /)
!Transform point back to real space for easier checking
v = matmul(orient, matmul(transform_matrix,v))
!If within the boundaries
if(in_block_bd(v, box_bd)) then
lat_points(ix-bd_in_lat(1)+5,iy-bd_in_lat(3)+5,iz-bd_in_lat(5) + 5) = .true.
end if
end do
end do
end do
!Now we redefine bd_in_lat The first 3 indices contains limits for the lat_points array
bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10
bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10
bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10
!Figure out where the starting point is. This is the first piont which fully contains the finite element
outerloop: do iz = 1, bd_in_array(3)
do iy = 1, bd_in_array(2)
do ix = 1, bd_in_array(1)
node_in_bd(8) = .false.
do inod = 1, 8
vlat = ele(:,inod) + (/ ix, iy, iz /)
!Check to see if the node resides at a position containing a lattice point within the box
if(any(vlat > shape(lat_points))) then
continue
else if(lat_points(vlat(1),vlat(2),vlat(3))) then
node_in_bd(inod) = .true.
end if
end do
if(all(node_in_bd)) then
rzero = (/ ix, iy, iz /)
exit outerloop
end if
end do
end do
end do outerloop
!Now build the finite element region
lat_ele_num = 0
lat_atom_num = 0
allocate(r_lat(3,8,numlatpoints/esize))
!Redefined the second 3 indices as the number of elements that fit in the bounds
do i = 1, 3
bd_in_array(3+i) = bd_in_array(i)/esize
end do
!Now start the element at rzero
do inod=1, 8
ele(:,inod) = ele(:,inod) + rzero
end do
do iz = -bd_in_array(6), bd_in_array(6)
do iy = -bd_in_array(5), bd_in_array(5)
do ix = -bd_in_array(4), bd_in_array(4)
node_in_bd(:) = .false.
temp_nodes(:,:,:) = 0.0_dp
temp_lat(:,:) = 0
do inod = 1, 8
vlat= ele(:,inod) + (/ ix*(esize), iy*(esize), iz*(esize) /)
!Transform point back to real space for easier checking
! v = matmul(orient, matmul(transform_matrix,v))
do i = 1,3
v(i) = real(vlat(i) + bd_in_lat(2*i-1) - 5)
end do
temp_nodes(:,1, inod) = matmul(orient, matmul(transform_matrix, v))
temp_lat(:,inod) = vlat
!Check to see if the lattice point values are greater then the array limits
if(any(vlat > shape(lat_points)).or.any(vlat < 1)) then
continue
!If within array boundaries check to see if it is a lattice point
else if(lat_points(vlat(1),vlat(2),vlat(3))) then
node_in_bd(inod) = .true.
end if
end do
if(all(node_in_bd)) then
lat_ele_num = lat_ele_num+1
r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)
!Now set all the lattice points contained within an element to false
do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))
do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:))
do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:))
lat_points(m,n,o) = .false.
end do
end do
end do
end if
end do
end do
end do
!Now figure out how many lattice points could not be contained in elements
allocate(r_atom_lat(3,count(lat_points)))
lat_atom_num = 0
do iz = 1, bd_in_array(3)
do iy = 1, bd_in_array(2)
do ix = 1, bd_in_array(1)
!If this point is a lattice point then save the lattice point as an atom
if (lat_points(ix,iy,iz)) then
v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /)
do i = 1,3
v(i) = v(i) + real(bd_in_lat(2*i-1) - 5)
end do
!Transform point back to real space
v = matmul(orient, matmul(transform_matrix,v))
lat_atom_num = lat_atom_num + 1
r_atom_lat(:,lat_atom_num) = v
end if
end do
end do
end do
end if
end subroutine build_with_rhomb
subroutine error_message(errorid)
integer, intent(in) :: errorid
select case(errorid)
case(1)
STOP "Name is missing."
case(2)
print *, "Element_type is missing"
case(3)
print *, "Lattice_parameter is missing"
case(4)
print *, "Lattice parameter is not in float form"
case(5)
print *, "Esize is missing"
case(6)
print *, "Esize is not in integer form"
case(7)
print *, "Cx(1) should equal Cx(2) in plane centroid finding algorithm. Please double check implementation"
end select
STOP 3
end subroutine error_message
end module mode_create

@ -0,0 +1,102 @@
module mode_merge
!This module contains the code needed for merging several .mb files together
use parameters
use atoms
use io
use subroutines
use elements
character(len=4) :: dim
integer :: in_num
public
contains
subroutine merge
integer :: i
real(kind=dp) :: displace(3), temp_box_bd(6)
!First we parse the merge command
call parse_command
!Now loop over all files and stack them
do i = 1, in_num
displace(:) = 0.0_dp
if ((i==1).or.(trim(adjustl(dim)) == 'none')) then
call read_in(i, displace, temp_box_bd)
call grow_box(temp_box_bd)
else
select case(trim(adjustl(dim)))
case('x')
displace(1) = box_bd(2)
case('y')
displace(2) = box_bd(4)
case('z')
displace(3) = box_bd(6)
end select
call read_in(i, displace, temp_box_bd)
call grow_box(temp_box_bd)
end if
end do
return
end subroutine merge
subroutine parse_command
character(len=100) :: textholder
integer :: i, stat, arglen, arg_pos
!Get dimension to concatenate along
call get_command_argument(2, dim, arglen)
if (arglen == 0) STOP "Missing dim in mode_merge command"
select case(trim(adjustl(dim)))
case('x','y','z','none')
continue
case default
print *, dim, " not accepted as a dimension in mode_merge"
stop 3
end select
!Get number of files to read in
call get_command_argument(3, textholder, arglen)
if (arglen == 0) STOP "Number of files to read missing in mode_merge command"
read(textholder, *, iostat = stat) in_num
if (stat > 0) STOP "Error reading number of files in, must be integer"
!Now loop and pull out all files
do i = 1, in_num
call get_command_argument(3+i, textholder, arglen)
if (arglen == 0) STOP "Fewer files to read in then specified"
!Make sure this file is readable
call get_in_file(textholder)
end do
!Set argpos accordingly
arg_pos = 3+in_num+1
!Now options loop
!Check for optional keywords
do while(.true.)
if(arg_pos > command_argument_count()) exit
!Pull out the next argument which should either be a keyword or an option
call get_command_argument(arg_pos, textholder)
textholder=adjustl(textholder)
arg_pos=arg_pos+1
!Choose what to based on what the option string is
select case(trim(textholder))
case default
!Check to see if it is an option command, if so then mode_create must be finished
if(textholder(1:1) == '-') then
exit
!Check to see if a filename was passed
elseif(scan(textholder,'.',.true.) > 0) then
call get_out_file(textholder)
end if
end select
end do
end subroutine parse_command
end module mode_merge

@ -0,0 +1,8 @@
module parameters
implicit none
integer, parameter :: dp= selected_real_kind(15,307)
real(kind=dp), parameter :: lim_zero = epsilon(1.0_dp), &
lim_large = huge(1.0_dp)
end module parameters

@ -0,0 +1,148 @@
module subroutines
use parameters
use functions
implicit none
integer :: allostat, deallostat
public
contains
!This subroutine is just used to break the code and exit on an error
subroutine read_error_check(para, loc)
integer, intent(in) :: para
character(len=100), intent(in) :: loc
if (para > 0) then
print *, "Read error in ", trim(loc), " because of ", para
stop "Exit with error"
end if
end subroutine
subroutine matrix_inverse(a, n, a_inv)
integer :: i, j, k, piv_loc
integer, intent(in) :: n
real(kind = dp) :: coeff, sum_l, sum_u
real(kind = dp), dimension(n) :: b, x, y, b_piv
real(kind = dp), dimension(n, n) :: l, u, p
real(kind = dp), dimension(n, n), intent(in) :: a
real(kind = dp), dimension(n, n), intent(out) :: a_inv
real(kind = dp), allocatable :: v(:), u_temp(:), l_temp(:), p_temp(:)
l(:, :) = identity_mat(n)
u(:, :) = a(:, :)
p(:, :) = identity_mat(n)
!LU decomposition with partial pivoting
do j = 1, n-1
allocate(v(n-j+1), stat = allostat)
if(allostat /=0 ) then
print *, 'Fail to allocate v in matrix_inverse'
stop
end if
v(:) = u(j:n, j)
if(maxval(abs(v)) < lim_zero) then
print *, 'Fail to inverse matrix', a
stop
end if
piv_loc = maxloc(abs(v), 1)
deallocate(v, stat = deallostat)
if(deallostat /=0 ) then
print *, 'Fail to deallocate v in matrix_inverse'
stop
end if
!partial pivoting
if(piv_loc /= 1) then
allocate( u_temp(n-j+1), p_temp(n), stat = allostat)
if(allostat /=0 ) then
print *, 'Fail to allocate p_temp and/or u_temp in matrix_inverse'
stop
end if
u_temp(:) = u(j, j:n)
u(j, j:n) = u(piv_loc+j-1, j:n)
u(piv_loc+j-1, j:n) = u_temp(:)
p_temp(:) = p(j, :)
p(j, :) = p(piv_loc+j-1, :)
p(piv_loc+j-1, :) = p_temp(:)
deallocate( u_temp, p_temp, stat = deallostat)
if(deallostat /=0 ) then
print *, 'Fail to deallocate p_temp and/or u_temp in matrix_inverse'
stop
end if
if(j > 1) then
allocate( l_temp(j-1), stat = allostat)
if(allostat /= 0) then
print *, 'Fail to allocate l_temp in matrix_inverse'
stop
end if
l_temp(:) = l(j, 1:j-1)
l(j, 1:j-1) = l(piv_loc+j-1, 1:j-1)
l(piv_loc+j-1, 1:j-1) = l_temp(:)
deallocate( l_temp, stat = deallostat)
if(deallostat /=0 ) then
print *, 'Fail to deallocate l_temp in matrix_inverse'
stop
end if
end if
end if
!LU decomposition
do i = j+1, n
coeff = u(i, j)/u(j, j)
l(i, j) = coeff
u(i, j:n) = u(i, j:n)-coeff*u(j, j:n)
end do
end do
a_inv(:, :) = 0.0_dp
do j = 1, n
b(:) = 0.0_dp
b(j) = 1.0_dp
b_piv(:) = matmul(p, b)
!Now we have LUx = b_piv
!the first step is to solve y from Ly = b_piv
!forward substitution
do i = 1, n
if(i == 1) then
y(i) = b_piv(i)/l(i, i)
else
sum_l = 0
do k = 1, i-1
sum_l = sum_l+l(i, k)*y(k)
end do
y(i) = (b_piv(i)-sum_l)/l(i, i)
end if
end do
!then we solve x from ux = y
!backward subsitution
do i = n, 1, -1
if(i == n) then
x(i) = y(i)/u(i, i)
else
sum_u = 0
do k = i+1, n
sum_u = sum_u+u(i, k)*x(k)
end do
x(i) = (y(i)-sum_u)/u(i, i)
end if
end do
!put x into j column of a_inv
a_inv(:, j) = x(:)
end do
return
end subroutine matrix_inverse
end module subroutines
Loading…
Cancel
Save