diff --git a/README.md b/README.md index 5bdd240..036981d 100644 --- a/README.md +++ b/README.md @@ -45,16 +45,6 @@ Default orientation is `[100] [010] [001]`. If this keyword is present then the *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** ``` @@ -68,7 +58,7 @@ Default duplicate is `1 1 1`. This is used to replicate the element along each d **Dimensions** ``` -dimensions dimx dimy dimz +dim 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. @@ -92,6 +82,13 @@ basis basisnum bname bx by bz ``` This function allows you to define a custom basis. `bname bx by bz` must be repeated `basisnum` times. + +**efill** +``` +efill xyz +``` +This command will rerun the creation algorithm with multiple times starting with an esize of `esize` and decreasing it by half on every iteration in an effort to maximize the reduction of degrees of freedom in the system. You must specify which dimensions will be filled. The code accepts `x`, `y`, `z`, `xy`, `yz`, `xz`, and `xyz` specifying which boundaries to fill in. + ### Mode Convert ``` @@ -132,6 +129,11 @@ wrap This will wrap atomic positions back inside the box. Effectively as if periodic boundary conditions are applied so that atoms which exit from one side of the simulation cell enter back in through the other. +###Mode Metric +``` +--metric cmetric nfiles {filepaths} +``` + ## Options Options are additional portions of code which have additional functionality. Options are performed in the order that they appear in the argument list and can be added to any mode. If wanting to use strictly options use `--convert` to specify input and output files. @@ -188,15 +190,15 @@ This option creates a circular planar vacancy cluster of radius `radius` normal `-group select_type group_shape shape_arguments additional keywords` -This option selects a group of either elements, nodes, atoms and applies some transformation to them. +This option selects a group of either elements or atoms and applies some transformation to them. -`select_type` - Either `nodes`, `atoms`, `elements`, `nodes/atoms`, `all`. When using the option `nodes` only nodes which are within the group are selected, `elements` selects elements based on whether the element center is within the group, `nodes/atoms` selects both nodes and atoms for the group. `all` selects elements based on the element center and atoms based on their position. + `select_type` - Either `atoms`, `elements`, or 'both'. `elements` selects elements based on whether the element center is within the group. `both` selects elements based on the element center and atoms based on their position. `group_shape` - Specifies what shape the group takes and dictates which options must be passed. Each shape requires different arguments and these arguments are represented by the placeholder `shape_arguments`. The accepted group shapes and arguments are below: *Block:* -`-group nodes block xlo xhi ylo yhi zlo zhi` +`-group atoms block xlo xhi ylo yhi zlo zhi` This selects a group residing in a block with edges perpendicular to the simulation cell. The block boundaries are given by `xlo xhi ylo yhi zlo zhi`. @@ -204,7 +206,7 @@ This selects a group residing in a block with edges perpendicular to the simulat *Wedge:* -`-group nodes wedge dim1 dim2 bx by bz bw` +`-group atoms wedge dim1 dim2 bx by bz bw` This selects a group which are within a wedge shape. The options are given as follows: `dim1` - The dimension containing the plane normal of the wedge base. `dim2` - The thickness dimension. Wedge groups are currently required to span the entire cell thickness in one dimensions which is normal to the triangular face. This through thickness dimension is dim2. @@ -213,7 +215,7 @@ This selects a group which are within a wedge shape. The options are given as fo *Notch:* -`-group nodes notch dim1 dim2 bx by bz bw tr` +`-group atoms notch dim1 dim2 bx by bz bw tr` This shape is similar to a wedge shape except instead of becoming atomically sharp, it finishes in a rounded tip with tip radius `tr`. Options are as follows. `dim1` - The dimension containing the plane normal of the wedge base. `dim2` - The thickness dimension. Wedge groups are currently required to span the entire cell thickness in one dimensions which is normal to the triangular face. This through thickness dimension is dim2. @@ -221,6 +223,11 @@ This shape is similar to a wedge shape except instead of becoming atomically sha `bw` - Base width `tr` - Tip radius +*Sphere* + +`-group atoms sphere x y z r` +This shape selects all atoms within a sphere centered at `(x,y,z)` with radius `r`. + **Displace:** ``` @@ -240,10 +247,10 @@ This command wraps atoms back into the simulation cell as though periodic bounda **Remesh** ``` -remesh esize lattice_parameter lattice_type +remesh esize ``` -This command remeshes the atoms/elements within the group to the new element size `esize`. Currently only accepts an `esize` of 2 which refines it to full atomistics. When remeshing to atomistics the group can contain any orientations of elements but when remeshing to different finite elements, the group must contain all atoms/elements with the same orientation. `lattice_parameter` is the lattice parameter for the elements within the group and `lattice_type` is the lattice type (integer) that these new elements will be assigned to. +This command remeshes the atoms/elements within the group to the new element size `esize`. **Max** @@ -261,6 +268,29 @@ delete This command deletes all selected atoms and elements within the group. + +**Random** +``` +random n +``` + +This command selects `n` random atoms and `n` random elements within your group bounds. If using group type `atoms` or `elements` then only `n` random atoms or elements are selected. This random atoms/elements then form the new group. + +**Nodes** + +``` +nodes +``` + +This keyword changes the selection criteria when using `elements` or `both` to element nodes instead of element centroids. + +**Flip** + +``` +flip +``` + +This keyword changes the group selection criteria from the atoms/elements inside a region to the atoms/elements outside the group. ### Option overwrite ``` @@ -303,8 +333,14 @@ This command will delete all overlapping atoms within a specific cutoff radius ` This option is primarily used when reading data from non .mb formats. This code simply sets the orientation variable for the specified sub box `sbox`. +### Option redef_box +``` +-redef_box redef_dim {xlo xhi} {ylo yhi} {zlo zhi} +``` +This option allows for the user to redefine box boundaries deleting atoms/elements outside of the new box boundaries. Elements are refined to atmoistics if they intersect the newly defined box boundaries. +The arguments are described below: +`redef_dim` - The dimensions to be redefined. Can be any permutation of `x`, `y`, `z`, `xy`, `yz`, `xz`, `xyz`. The order of the dimensions dictates the order that the new dimensions must be defined **** - ## Position Specification Specifying positions in cacmb can be done through a variety of ways. Examples of each format is shown below. diff --git a/src/Makefile b/src/Makefile index 642fce5..158d9db 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,44 +1,49 @@ -FC=ifort -#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays -FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays -MODES=mode_create.o mode_merge.o mode_convert.o -OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o -OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES) $(OPTIONS) call_option.o sorts.o - -.SUFFIXES: .SUFFIXES: .c .f .f90 .F90 .o +.DEFAULT_GOAL := all -cacmb: $(OBJECTS) - $(FC) $(FFLAGS) $(OBJECTS) parameters.o -o $@ +FC=mpif90 +FFLAGS=-Wall -mcmodel=large -O0 -g -fbacktrace -fcheck=all -ffpe-trap=invalid,zero,overflow,underflow,denormal + +OBJDIR=obj +SRCS := $(wildcard *.f90) +OBJECTS := $(addprefix $(OBJDIR)/,$(SRCS:%.f90=%.o)) + + + +#----------------- DEPENDENCIES -----------------# +# GENERATED USING https://github.com/ZedThree/fort_depend.py **requires python3** +# > pip install fortdepend +# > fortdepend -o Makefile.dep -i mpi -b obj/ +include Makefile.dep + +#----------------- DEFAULTS -----------------# +all: cacmb + +.PHONY: deps +cacmb: $(OBJECTS) $(OBJDIR)/main.o + $(FC) $(FFLAGS) $(OBJECTS) -o $@ + +$(OBJDIR)/%.o: %.f90 + @mkdir -p $(@D) + $(FC) $(FFLAGS) -c -o $@ $< -J$(OBJDIR) .f90.o: $(FC) $(FFLAGS) -c $< + +deps: + @fortdepend -o Makefile.dep -i mpi -b obj -w + +#----------------- CLEAN UP -----------------# + .PHONY: clean -clean: - $(RM) cacmb *.o - -testfuncs: testfuncs.o functions.o subroutines.o - $(FC) testfuncs.o functions.o subroutines.o box.o elements.o -o $@ - -.PHONY: cleantest -cleantest: - $(RM) testfuncs testfuncs.o - -.PHONY: test -test: testfuncs - ./testfuncs - -.PHONY: install -install: cacmb - cp ./cacmb /usr/local/bin - -$(OBJECTS) : parameters.o -atoms.o subroutines.o testfuncs.o box.o : functions.o -main.o io.o $(MODES) $(OPTIONS) : elements.o -call_mode.o : $(MODES) -call_option.o : $(OPTIONS) -elements.o : sorts.o -$(MODES) $(OPTIONS) subroutines.o io.o : atoms.o box.o -$(MODES) main.o : io.o -testfuncs.o elements.o mode_create.o $(OPTIONS) $(MODES): subroutines.o +clean: + $(RM) *.mod *.o + $(RM) $(OBJDIR)/*.mod $(OBJDIR)/*.o CAC + @$(RM) -rf obj/ + +.PHONY: clean-all +clean-all: clean + +# DEBUGGING VARIABLE PRINT +print-% : ; @echo $* = $($*) diff --git a/src/Makefile.dep b/src/Makefile.dep new file mode 100644 index 0000000..54ee834 --- /dev/null +++ b/src/Makefile.dep @@ -0,0 +1,167 @@ +# This file is generated automatically. DO NOT EDIT! + +obj/main : \ + obj/atoms.o \ + obj/box.o \ + obj/caller.o \ + obj/elements.o \ + obj/functions.o \ + obj/io.o \ + obj/main.o \ + obj/mode_calc.o \ + obj/mode_convert.o \ + obj/mode_create.o \ + obj/mode_merge.o \ + obj/mode_metric.o \ + obj/neighbors.o \ + obj/opt_deform.o \ + obj/opt_delete.o \ + obj/opt_disl.o \ + obj/opt_group.o \ + obj/opt_orient.o \ + obj/opt_redef_box.o \ + obj/opt_slip_plane.o \ + obj/parameters.o \ + obj/sorts.o \ + obj/str.o \ + obj/subroutines.o + +obj/atoms.o : \ + obj/functions.o \ + obj/parameters.o + +obj/box.o : \ + obj/functions.o \ + obj/parameters.o + +obj/caller.o : \ + obj/box.o \ + obj/mode_calc.o \ + obj/mode_convert.o \ + obj/mode_create.o \ + obj/mode_merge.o \ + obj/mode_metric.o \ + obj/opt_deform.o \ + obj/opt_delete.o \ + obj/opt_disl.o \ + obj/opt_group.o \ + obj/opt_orient.o \ + obj/opt_redef_box.o \ + obj/opt_slip_plane.o \ + obj/parameters.o + +obj/elements.o : \ + obj/box.o \ + obj/functions.o \ + obj/parameters.o \ + obj/sorts.o \ + obj/subroutines.o + +obj/functions.o : \ + obj/parameters.o + +obj/io.o : \ + obj/atoms.o \ + obj/box.o \ + obj/elements.o \ + obj/parameters.o \ + obj/str.o + +obj/main.o : \ + obj/caller.o \ + obj/elements.o \ + obj/io.o \ + obj/parameters.o + +obj/mode_calc.o : \ + obj/box.o \ + obj/elements.o \ + obj/io.o \ + obj/parameters.o \ + obj/subroutines.o + +obj/mode_convert.o : \ + obj/box.o \ + obj/elements.o \ + obj/io.o \ + obj/parameters.o + +obj/mode_create.o : \ + obj/atoms.o \ + obj/box.o \ + obj/elements.o \ + obj/io.o \ + obj/parameters.o \ + obj/subroutines.o + +obj/mode_merge.o : \ + obj/atoms.o \ + obj/elements.o \ + obj/io.o \ + obj/parameters.o \ + obj/subroutines.o + +obj/mode_metric.o : \ + obj/elements.o \ + obj/io.o \ + obj/neighbors.o \ + obj/parameters.o + +obj/neighbors.o : \ + obj/elements.o \ + obj/functions.o \ + obj/parameters.o \ + obj/subroutines.o + +obj/opt_deform.o : \ + obj/box.o \ + obj/elements.o \ + obj/parameters.o \ + obj/subroutines.o + +obj/opt_delete.o : \ + obj/elements.o \ + obj/neighbors.o \ + obj/parameters.o \ + obj/subroutines.o + +obj/opt_disl.o : \ + obj/box.o \ + obj/elements.o \ + obj/parameters.o \ + obj/subroutines.o + +obj/opt_group.o : \ + obj/box.o \ + obj/elements.o \ + obj/parameters.o \ + obj/subroutines.o + +obj/opt_orient.o : \ + obj/box.o \ + obj/elements.o \ + obj/parameters.o \ + obj/subroutines.o + +obj/opt_redef_box.o : \ + obj/box.o \ + obj/elements.o \ + obj/subroutines.o + +obj/opt_slip_plane.o : \ + obj/elements.o \ + obj/functions.o \ + obj/parameters.o \ + obj/subroutines.o + +obj/parameters.o : + +obj/sorts.o : \ + obj/parameters.o + +obj/str.o : + +obj/subroutines.o : \ + obj/box.o \ + obj/functions.o \ + obj/parameters.o diff --git a/src/box.f90 b/src/box.f90 index c6e7948..9f2c7b2 100644 --- a/src/box.f90 +++ b/src/box.f90 @@ -14,7 +14,6 @@ module 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 - !Below are some simulation parameters which may be adjusted if reading in restart files integer :: timestep=0 real(kind=dp) :: total_time=0.0_dp @@ -74,10 +73,15 @@ module box 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 + if(all(abs(box_bd) < lim_zero)) then + box_bd = temp_box_bd + else + 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 + end if + return end subroutine grow_box @@ -98,4 +102,17 @@ module box end do return end subroutine in_sub_box -end module box \ No newline at end of file + + subroutine reset_box + !This subroutine just resets the box boundary and position + box_bc = "ppp" + box_bd(:) = 0 + end subroutine reset_box + + pure function box_volume() + real(kind = dp) :: box_volume + + box_volume = (box_bd(2) - box_bd(1)) * (box_bd(4) - box_bd(3))*(box_bd(6) - box_bd(5)) + return + end function +end module box diff --git a/src/call_mode.f90 b/src/call_mode.f90 deleted file mode 100644 index fb187c1..0000000 --- a/src/call_mode.f90 +++ /dev/null @@ -1,29 +0,0 @@ -subroutine call_mode(arg_pos,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(out) :: arg_pos - character(len=100), intent(in) :: mode - - select case(mode) - case('--create') - call create(arg_pos) - case('--convert') - call convert(arg_pos) - case('--merge') - call merge(arg_pos) - 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 diff --git a/src/call_option.f90 b/src/call_option.f90 deleted file mode 100644 index ead7af0..0000000 --- a/src/call_option.f90 +++ /dev/null @@ -1,40 +0,0 @@ -subroutine call_option(option, arg_pos) - use parameters - use opt_disl - use opt_group - use opt_orient - use opt_delete - use box - implicit none - - integer, intent(inout) :: arg_pos - character(len=100), intent(in) :: option - - select case(trim(adjustl(option))) - case('-dislgen', '-disloop','-vacancydisloop') - call dislocation(option, arg_pos) - case('-group') - call group(arg_pos) - case('-ow') - arg_pos = arg_pos + 1 - case('-wrap') - arg_pos = arg_pos + 1 - case('-orient') - call orient(arg_pos) - case('-unorient') - call unorient - arg_pos = arg_pos + 1 - case('-boundary') - arg_pos=arg_pos+1 - call get_command_argument(arg_pos, box_bc) - arg_pos=arg_pos+1 - bound_called = .true. - case('-sbox_ori') - call sbox_ori(arg_pos) - case('-delete') - call run_delete(arg_pos) - case default - print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.' - stop 3 - end select -end subroutine call_option \ No newline at end of file diff --git a/src/caller.f90 b/src/caller.f90 new file mode 100644 index 0000000..df29501 --- /dev/null +++ b/src/caller.f90 @@ -0,0 +1,91 @@ +module caller + !this module just calls modes and options + + use mode_create + use mode_convert + use mode_merge + use mode_metric + use mode_calc + use parameters + use opt_disl + use opt_group + use opt_orient + use opt_deform + use opt_delete + use opt_redef_box + use opt_slip_plane + use box + + + implicit none + public + contains + subroutine call_mode(arg_pos) + !This code is used to parse the command line argument for the mode information and calls the required + !mode module. + + integer, intent(out) :: arg_pos + + select case(mode) + case('--create') + call create(arg_pos) + case('--convert') + call convert(arg_pos) + case('--merge') + call merge(arg_pos) + case('--metric') + call metric(arg_pos) + case('--calc') + call calc(arg_pos) + 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 + + subroutine call_option(option, arg_pos) + integer, intent(inout) :: arg_pos + character(len=100), intent(in) :: option + + select case(trim(adjustl(option))) + case('-disl','-dislgen', '-disloop','-vacancydisloop') + call dislocation(option, arg_pos) + case('-group') + call group(arg_pos) + case('-ow') + arg_pos = arg_pos + 1 + case('-wrap') + arg_pos = arg_pos + 1 + case('-orient') + call orient_opt(arg_pos) + case('-unorient') + call unorient + arg_pos = arg_pos + 1 + case('-boundary') + arg_pos=arg_pos+1 + call get_command_argument(arg_pos, box_bc) + arg_pos=arg_pos+1 + bound_called = .true. + case('-sbox_ori') + call sbox_ori(arg_pos) + case('-deform') + call deform(arg_pos) + case('-delete') + call run_delete(arg_pos) + case('-set_cac') + arg_pos=arg_pos +3 + case('-set_types') + arg_pos = arg_pos + 3 + atom_types + case('-redef_box') + call redef_box(arg_pos) + case('-slip_plane') + call run_slip_plane(arg_pos) + case default + print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.' + stop 3 + end select + end subroutine call_option + +end module caller diff --git a/src/elements.f90 b/src/elements.f90 index 0dffecc..a70af56 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -11,24 +11,29 @@ 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(:), sbox_ele(:) !Element size + integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:), tag_ele(:) !Element size real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array + !Element result data structures + real(kind=dp), allocatable :: force_node(:,:,:,:), virial_node(:,:,:,:), energy_node(:,:,:) integer, save :: ele_num !Number of elements integer, save :: node_num !Total number of nodes + integer, save :: node_atoms !Count of all basis atoms at nodes summed over all nodes !Data structure used to represent atoms integer, allocatable :: type_atom(:)!atom type - integer, allocatable :: sbox_atom(:) + integer, allocatable :: sbox_atom(:), tag_atom(:) real(kind =dp),allocatable :: r_atom(:,:) !atom position integer :: atom_num=0 !Number of atoms + !Atom result data structures information + real(kind=8), allocatable :: force_atom(:,:), virial_atom(:,:), energy_atom(:) !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) + 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) integer :: cubic_faces(4,6) !Below are lattice type arrays which provide information on the general form of the elements. @@ -36,16 +41,19 @@ module elements 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) !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) - real(kind=dp) :: lapa(10) !Additional module level variables we need logical :: wrap_flag + + !flags for data variables + logical :: vflag public contains @@ -87,12 +95,33 @@ module elements 0.0_dp, 0.5_dp, 0.5_dp, & 0.5_dp, 0.0_dp, 0.5_dp /), & shape(fcc_mat)) + + !Initialize the bcc primitive cell + bcc_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, & + 0.5_dp, -0.5_dp, 0.5_dp, & + 1.0_dp, 0.0_dp, 1.0_dp, & + 0.5_dp, 0.5_dp, 0.5_dp, & + -0.5_dp, 0.5_dp, 0.5_dp, & + 0.0_dp, 0.0_dp, 1.0_dp, & + 0.5_dp, 0.5_dp, 1.5_dp, & + 0.0_dp, 1.0_dp, 1.0_dp /), & + shape(bcc_cell)) + + bcc_mat = reshape((/ 0.5_dp, -0.5_dp, 0.5_dp, & + 0.5_dp, 0.5_dp, 0.5_dp, & + -0.5_dp, 0.5_dp, 0.5_dp /), & + shape(bcc_mat)) + + + call matrix_inverse(fcc_mat,3,fcc_inv) + call matrix_inverse(bcc_mat,3,bcc_inv) max_basisnum = 0 basisnum(:) = 0 ng_node(:) = 0 node_num = 0 + node_atoms = 0 ele_num = 0 atom_num = 0 end subroutine lattice_init @@ -112,6 +141,7 @@ module elements real(kind=dp), dimension(3,max_ng_node) :: adjustVar adjustVar(:,:) = 0.0_dp + vflag = .false. select case(trim(ele_type)) @@ -130,6 +160,9 @@ module elements end if cell_mat(:, 1:8) = fcc_cell + adjustVar(:,1:8) cell_mat(:,1:8) = lapa * ((esize-1)*matmul(orient_mat, cell_mat(:,1:8))) + case('bcc') + cell_mat(:,1:8) = bcc_cell + cell_mat(:,1:8) = lapa* ((esize-1)*matmul(orient_mat, cell_mat(:,1:8))) case default print *, "Element type ", trim(ele_type), " currently not accepted" stop @@ -146,7 +179,7 @@ module elements !Allocate element arrays if(n > 0) then - allocate(type_ele(n), size_ele(n), lat_ele(n), sbox_ele(n), r_node(3,max_basisnum, max_ng_node,n), & + allocate(type_ele(n), tag_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 @@ -156,7 +189,7 @@ module elements if(m > 0) then !Allocate atom arrays - allocate(type_atom(m), sbox_atom(m), r_atom(3,m), stat=allostat) + allocate(type_atom(m), sbox_atom(m), tag_atom(m), r_atom(3,m), stat=allostat) if(allostat > 0) then print *, "Error allocating atom arrays in elements.f90 because of: ", allostat stop @@ -175,84 +208,124 @@ module elements !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, 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, 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 - 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) + + !First check to make sure if it is allocated + if (allocated(size_ele)) then + + !Figure out the size of the atom and element arrays + ele_size = size(size_ele) + + !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_size+buffer_size)) + temp_int(1:ele_size) = lat_ele(1:ele_size) + temp_int(ele_size+1:) = 0 + call move_alloc(temp_int, lat_ele) + + allocate(temp_int(n+ele_size+buffer_size)) + temp_int(1:ele_size) = tag_ele(1:ele_size) + temp_int(ele_size+1:) = 0 + call move_alloc(temp_int, tag_ele) + + allocate(temp_int(n+ele_size+buffer_size)) + temp_int(1:ele_size) = size_ele(1:ele_size) + temp_int(ele_size+1:) = 0 + call move_alloc(temp_int, size_ele) + + allocate(temp_int(n+ele_size+buffer_size)) + temp_int(1:ele_size) = sbox_ele(1:ele_size) + temp_int(ele_size+1:) = 0 + call move_alloc(temp_int, sbox_ele) + + allocate(char_temp(n+ele_size+buffer_size)) + char_temp(1:ele_size) = type_ele(1:ele_size) + call move_alloc(char_temp, type_ele) + + allocate(temp_ele_real(3, max_basisnum, max_ng_node, n+ele_size+buffer_size)) + temp_ele_real(:,:,:,1:ele_size) = r_node(:,:,:,1:ele_size) + temp_ele_real(:,:,:,ele_size+1:) = 0.0_dp + call move_alloc(temp_ele_real, r_node) + end if + else + call alloc_ele_arrays(n,0) 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_int(m+atom_num+buffer_size)) - temp_int(1:atom_size) = sbox_atom - temp_int(atom_size+1:) = 0 - call move_alloc(temp_int, sbox_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) + if (allocated(type_atom)) then + atom_size = size(type_atom) + if (m+atom_num > atom_size) then + allocate(temp_int(m+atom_size+buffer_size)) + temp_int(1:atom_size) = type_atom + temp_int(atom_size+1:) = 0 + call move_alloc(temp_int, type_atom) + + allocate(temp_int(m+atom_size+buffer_size)) + temp_int(1:atom_size) = tag_atom + temp_int(atom_size+1:) = 0 + call move_alloc(temp_int, tag_atom) + + allocate(temp_int(m+atom_size+buffer_size)) + temp_int(1:atom_size) = sbox_atom + temp_int(atom_size+1:) = 0 + call move_alloc(temp_int, sbox_atom) + + allocate(temp_real(3,m+atom_size+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 + else + call alloc_ele_arrays(0,m) end if end subroutine - subroutine add_element(type, size, lat, sbox, r) + subroutine add_element(tag, type, size, lat, sbox, r) !Subroutine which adds an element to the element arrays - integer, intent(in) :: size, lat, sbox + integer, intent(in) :: size, lat, sbox, tag character(len=100), intent(in) :: type real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node) + integer :: newtag + ele_num = ele_num + 1 + node_num = node_num + ng_node(lat) + node_atoms = node_atoms + ng_node(lat)*basisnum(lat) + + if (tag==0) then + newtag = ele_num !If we don't assign a tag then pass the tag as the ele_num + else + newtag = tag + end if + !Check to see if we need to grow the arrays call grow_ele_arrays(1,0) + tag_ele(ele_num) = newtag 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) + r_node(:,:,:,ele_num) = r(:,:,:) end subroutine add_element - subroutine add_atom(type, sbox, r) + subroutine add_atom(tag, type, sbox, r) !Subroutine which adds an atom to the atom arrays - integer, intent(in) :: type, sbox + integer, intent(in) :: type, sbox, tag real(kind=dp), intent(in), dimension(3) :: r + integer :: newtag + atom_num = atom_num+1 + if(tag==0) then + newtag = atom_num !If we don't assign a tag then pass the tag as the atom_num + else + newtag = tag + end if !Check to see if we need to grow the arrays call grow_ele_arrays(0,1) + tag_atom(atom_num) = tag type_atom(atom_num) = type r_atom(:,atom_num) = r(:) sbox_atom(atom_num) = sbox @@ -295,12 +368,12 @@ module elements integer :: i - max_ng_node = 0 - do i=1, n select case(trim(adjustl(element_types(i)))) case('fcc') ng_node(i) = 8 + case('bcc') + ng_node(i) = 8 end select if(ng_node(i) > max_ng_node) max_ng_node = ng_node(i) @@ -309,10 +382,14 @@ module elements subroutine set_max_esize !This subroutine sets the maximum esize - max_esize=maxval(size_ele) + if(allocated(size_ele)) then + max_esize=maxval(size_ele) + else + max_esize = 2 + end if end subroutine - subroutine interpolate_atoms(type, esize, lat_type, r_in, type_interp, r_interp) + subroutine interpolate_atoms(type, esize, lat_type, r_in, type_interp, r_interp, eng, f, v, data_interp) !This subroutine returns the interpolated atoms from the elements. !Arguments @@ -322,15 +399,19 @@ module elements 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 + real(kind = dp), optional, intent(in) :: eng(max_basisnum, max_ng_node), f(3, max_basisnum, max_ng_node), & + v(6, max_basisnum, max_ng_node) + 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 - real(kind=dp), allocatable :: a_shape(:) + real(kind=dp) :: a_shape(max_ng_node) real(kind=dp) :: t, s, r !Initialize some variables r_interp(:,:) = 0.0_dp type_interp(:) = 0 + if(present(data_interp)) data_interp = 0.0_dp ia = 0 !Define bnum based on the input lattice type. If lat_type=0 then we are interpolating lattice points which means @@ -346,8 +427,7 @@ module elements end select select case(trim(adjustl(type))) - case('fcc') - allocate(a_shape(8)) + 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) @@ -363,6 +443,12 @@ module elements type_interp(ia) = basis_type(ibasis,lat_type_temp) r_interp(:,ia) = r_interp(:,ia) + a_shape(inod) * r_in(:,ibasis,inod) + if(present(data_interp)) then + !If data is present then interpolate data arrays as well + data_interp(1,ia) = data_interp(1,ia) + eng(ibasis, inod)*a_shape(inod) + data_interp(2:4,ia) = data_interp(2:4,ia) + f(:, ibasis, inod)*a_shape(inod) + data_interp(5:10,ia) = data_interp(5:10,ia) + v(:, ibasis, inod)*a_shape(inod) + end if end do end do @@ -381,7 +467,7 @@ module elements 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) + real(kind=8), intent(out) :: shape_fun(:) 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 @@ -446,19 +532,21 @@ module elements !We go from largest index to smallest index just to make sure that we don't miss !accidentally overwrite values which need to be deleted do i = num, 1, -1 + node_num = node_num - ng_node(lat_ele(sorted_index(i))) if(sorted_index(i) == ele_num) then r_node(:,:,:,sorted_index(i)) = 0.0_dp type_ele(sorted_index(i)) ='' size_ele(sorted_index(i)) = 0 lat_ele(sorted_index(i)) = 0 sbox_ele(sorted_index(i)) = 0 + tag_ele(sorted_index(i)) = 0 else - node_num = node_num - ng_node(lat_ele(sorted_index(i))) r_node(:,:,:,sorted_index(i)) = r_node(:,:,:,ele_num) type_ele(sorted_index(i)) = type_ele(ele_num) size_ele(sorted_index(i)) = size_ele(ele_num) lat_ele(sorted_index(i)) = lat_ele(ele_num) sbox_ele(sorted_index(i)) = sbox_ele(ele_num) + tag_ele(sorted_index(i)) = tag_ele(ele_num) end if ele_num = ele_num - 1 end do @@ -481,8 +569,7 @@ module elements max_bd(:) = -huge(1.0_dp) min_bd(:) = huge(1.0_dp) - - do i = 1, atom_num +do i = 1, atom_num do j = 1, 3 if (r_atom(j,i) > max_bd(j)) max_bd(j) = r_atom(j,i) + tol if (r_atom(j,i) < min_bd(j)) min_bd(j) = r_atom(j,i) - tol @@ -616,17 +703,17 @@ module elements esize = size_ele(ie) select case(iface) case(1) - pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp /) + pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, -10.0_dp**(-2.0_dp) /) case(2) - pos = (/ real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /) + pos = (/ real(esize-1,dp)/2.0_dp, -10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp /) case(3) - pos = (/ (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) + pos = (/ (esize-1)+10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) case(4) - pos = (/ real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /) + pos = (/ real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp /) case(5) - pos = (/ -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) + pos = (/ -10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) case(6) - pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp /) + pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**(-2.0_dp) /) end select !Now transform it to real space and adjust it to the position of the element in the first node. @@ -648,4 +735,131 @@ module elements end select end subroutine + subroutine lattice_map(in_bnum, in_btypes, in_ngnodes, in_lapa, lat_type) + !This subroutine maps an input lattice type to either a new lattice type or an existing one depending on basis_type and + !number of nodes at the atoms + + integer, intent(in) :: in_ngnodes, in_bnum, in_btypes(10) !Input variables + real(kind=dp), intent(in) :: in_lapa + integer, intent(out) :: lat_type + + integer j, ibasis + + lat_type = 0 + lat_loop:do j = 1, lattice_types + !Check all the lattice level variables + if ((basisnum(j) == in_bnum).and.(ng_node(j) == in_ngnodes).and.(is_equal(lapa(j),in_lapa))) then + !Now check lattice level variables + do ibasis = 1, basisnum(j) + if(basis_type(ibasis,j) /= in_btypes(ibasis)) cycle lat_loop + end do + lat_type = j + exit lat_loop + end if + end do lat_loop + + !If it doesn't match an existing lattice type we add it + if( lat_type == 0) then + lattice_types = lattice_types + 1 + basisnum(lattice_types) = in_bnum + basis_type(:,lattice_types) = in_btypes + ng_node(lattice_types) = in_ngnodes + lapa(lattice_types) = in_lapa + lat_type = lattice_types + end if + + end subroutine lattice_map + + subroutine get_interp_pos(i,j,k, ie, rout) + !This returns the position of an interpolated basis from an element ie. + !i, j, k should be in natural coordinates + + integer, intent(in) :: i, j, k + real(kind=dp), dimension(3,max_basisnum), intent(out) :: rout + + integer :: ie, ibasis, inod + real(kind=dp) :: a_shape(8), r, s, t + + r = (1.0_dp*(i-1)-(size_ele(ie)-1)/2)/(1.0_dp*(size_ele(ie)-1)/2) + s = (1.0_dp*(j-1)-(size_ele(ie)-1)/2)/(1.0_dp*(size_ele(ie)-1)/2) + t = (1.0_dp*(k-1)-(size_ele(ie)-1)/2)/(1.0_dp*(size_ele(ie)-1)/2) + rout(:,:) = 0 + do ibasis = 1, basisnum(lat_ele(ie)) + do inod = 1, ng_node(lat_ele(ie)) + call rhombshape(r,s,t,a_shape) + rout(:,ibasis) = rout(:,ibasis) + a_shape(inod) * r_node(:,ibasis,inod,ie) + end do + end do + + end subroutine get_interp_pos + + subroutine alloc_dat_arrays(n,m) + !This subroutine used to provide initial allocation for the atom and element data 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 + if (allocated(force_node)) then + deallocate(force_node, virial_node, energy_node) + end if + allocate(force_node(3,max_basisnum, max_ng_node, n), & + virial_node(6,max_basisnum, max_ng_node, n), & + energy_node(max_basisnum,max_ng_node,n), & + stat=allostat) + if(allostat > 0) then + print *, "Error allocating element data arrays in mode_metric because of:", allostat + stop + end if + + end if + + if (m > 0) then + if (allocated(force_atom)) then + deallocate(force_atom, virial_atom, energy_atom) + end if + allocate(force_atom(3, m), & + virial_atom(6, m), & + energy_atom(m), & + stat=allostat) + if(allostat > 0) then + print *, "Error allocating atom data arrays in mode_metric because of:", allostat + stop + end if + end if + + end subroutine + + subroutine add_atom_data(ia, eng, force, virial) + !Function which sets the atom data arrays + integer, intent(in) :: ia + real(kind=dp), intent(in) :: eng, force(3), virial(6) + + energy_atom(ia) = eng + force_atom(:,ia) = force(:) + virial_atom(:,ia) = virial(:) + vflag = .true. + return + end subroutine add_atom_data + + subroutine add_element_data(ie, eng, force, virial) + !Function which sets the element data arrays + integer, intent(in) :: ie + real(kind=dp), intent(in) :: eng(max_basisnum, max_ng_node), & + force(3,max_basisnum, max_ng_node), & + virial(6,max_basisnum,max_ng_node) + vflag = .true. + energy_node(:,:,ie) = eng + force_node(:,:,:,ie) = force + virial_node(:,:,:,ie) = virial + return + end subroutine add_element_data + + subroutine reset_data + !This function resets all of the data arrays for the elements and atoms + atom_num = 0 + ele_num = 0 + node_num = 0 + end subroutine reset_data + end module elements diff --git a/src/functions.f90 b/src/functions.f90 index 6fed2b7..ebd9a50 100644 --- a/src/functions.f90 +++ b/src/functions.f90 @@ -271,4 +271,132 @@ END FUNCTION StrDnCase norm_dis(1:3) = (rk - rl) norm_dis(4) = norm2(rk-rl) end function + + pure function matinv3(A) result(B) + !! Performs a direct calculation of the inverse of a 3×3 matrix. + real(kind=dp), intent(in) :: A(3,3) !! Matrix + real(kind=dp) :: B(3,3) !! Inverse matrix + real(kind=dp) :: detinv + + if(abs(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)& + - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)& + + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)) < lim_zero) then + B(:,:) = 0 + return + else + ! Calculate the inverse determinant of the matrix + + detinv = 1/(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)& + - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)& + + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)) + + ! Calculate the inverse of the matrix + B(1,1) = +detinv * (A(2,2)*A(3,3) - A(2,3)*A(3,2)) + B(2,1) = -detinv * (A(2,1)*A(3,3) - A(2,3)*A(3,1)) + B(3,1) = +detinv * (A(2,1)*A(3,2) - A(2,2)*A(3,1)) + B(1,2) = -detinv * (A(1,2)*A(3,3) - A(1,3)*A(3,2)) + B(2,2) = +detinv * (A(1,1)*A(3,3) - A(1,3)*A(3,1)) + B(3,2) = -detinv * (A(1,1)*A(3,2) - A(1,2)*A(3,1)) + B(1,3) = +detinv * (A(1,2)*A(2,3) - A(1,3)*A(2,2)) + B(2,3) = -detinv * (A(1,1)*A(2,3) - A(1,3)*A(2,1)) + B(3,3) = +detinv * (A(1,1)*A(2,2) - A(1,2)*A(2,1)) + end if + end function + + pure function transpose3(A) result(B) + !!Transposes matrix A + real(kind=dp), intent(in) :: A(3,3) + real(kind=dp) :: B(3,3) + + integer :: i, j + forall(i =1:3,j=1:3) B(i,j) = A(j,i) + + end function transpose3 + + pure function sqrt3(A) result(B) + !This calculates the square of matrix A fulfilling the equation B*B = A according to + !the algorithm by Franca1989 + + real(kind=dp), intent(in) :: A(3,3) + real(kind=dp) :: B(3,3) + + real(kind=dp) :: Ione, Itwo, Ithree, l, k, phi, Asq(3,3), lambda, Bone, Btwo, Bthree, p + + !Step 1 is calculating the invariants of C + Ione = A(1,1) + A(2,2) + A(3,3) + Asq = matmul(A,A) + Itwo = 0.5_dp *(Ione*Ione - (Asq(1,1) + Asq(2,2) + Asq(3,3))) + Ithree = (A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)& + - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)& + + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)) + + if (Ithree < 0) then + B(:,:)=0.0_dp + return + end if + !Check for an isotropic matrix + k = Ione*Ione - 3*Itwo + if (k < lim_zero) then + lambda = sqrt(Ione/3.0_dp) + B = lambda*identity_mat(3) + else + l = Ione*(Ione*Ione - 9.0_dp/2.0_dp * Itwo) + 27.0_dp/2.0_dp * Ithree + p = l/(k**(1.5_dp)) + + if (p > 1.0 ) then + B(:,:) = 0.0_dp + return + end if + + if ((p< -1).or.(p>1)) then + B(:,:) = 0.0_dp + return + end if + phi = acos(p) + lambda = sqrt(1.0_dp/3.0_dp * (Ione + 2*sqrt(k)*cos(phi/3))) + + !Now calculate invariantes of B + Bthree = sqrt(Ithree) + if((-lambda*lambda + Ione + 2*Ithree/lambda) < 0) then + B(:,:) = 0.0_dp + return + end if + Bone = lambda + sqrt(- lambda*lambda + Ione + 2*Ithree/lambda) + Btwo = (Bone*Bone - Ione)/2.0_dp + + !Now calculate B + if(abs(Bone*Btwo -Bthree) < lim_zero) then + B(:,:) = 0.0_dp + return + end if + B = 1/(Bone*Btwo - Bthree) *(Bone*Bthree*identity_mat(3) + (Bone*Bone - Btwo)*A - matmul(A,A)) + end if + end function sqrt3 + + pure function permutation(i,j,k) result(e) + !Calculates the permutation tensor + integer, intent(in) :: i,j,k + integer :: e + + if ( ((i==1).and.(j==2).and.(k==3)).or. & + ((i==2).and.(j==3).and.(k==1)).or. & + ((i==3).and.(j==1).and.(k==2))) then + e=1 + else if( ((i==2).and.(j==1).and.(k==3)).or. & + ((i==1).and.(j==3).and.(k==2)).or. & + ((i==3).and.(j==2).and.(k==1))) then + e=-1 + else + e=0 + end if + end function permutation + + pure function evtogp(virial) + real(kind=dp), dimension(6), intent(in) :: virial + real(kind=dp), dimension(6) :: evtogp + + evtogp = virial * 1e21_dp * 1.602176565e-19_dp + + end function + end module functions diff --git a/src/io.f90 b/src/io.f90 index d178baf..f8dd0b4 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -4,13 +4,14 @@ module io use parameters use atoms use box + use str implicit none integer :: outfilenum = 0, infilenum = 0 - character(len=100) :: outfiles(100), infiles(100) + character(len=100) :: outfiles(100), infiles(100), in_lattice_type='' logical :: force_overwrite - + real(kind=dp) :: in_lapa=0.0 public contains @@ -59,7 +60,7 @@ module io cycle end if select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:)) - case('xyz', 'lmp', 'vtk', 'mb', 'restart') + case('xyz', 'lmp', 'vtk', 'mb', 'restart', 'dump') outfilenum=outfilenum+1 outfiles(outfilenum) = temp_outfile exit @@ -104,9 +105,11 @@ module io call write_pycac(outfiles(i)) case('cac') call write_lmpcac(outfiles(i)) + case('dump') + call write_ldump(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" + " is not accepted for writing. Please select from: xyz, lmp, vtk, mb, restart, cac and try again" stop end select @@ -119,30 +122,41 @@ module io !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 :: i, inod, ibasis + integer :: i, inod, ibasis, outn open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') !Write total number of atoms + elements - write(11, '(i16)') node_num+atom_num + write(11, '(i16)') node_atoms+atom_num !Write comment line write(11, '(a)') "#Node + atom file created using cacmb" + outn=0 !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, '(i16, 3f23.15)') basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i) + write(11, '(2i16, 3f23.15)') basis_type(ibasis,lat_ele(i)), 1, r_node(:,ibasis,inod,i) + outn = outn + 1 end do end do end do + if(outn /= node_atoms) then + print *, "outn", outn, " doesn't equal node_atoms ", node_atoms + end if + !Write atom positions do i = 1, atom_num - write(11, '(i16, 3f23.15)') type_atom(i), r_atom(:,i) + write(11, '(2i16, 3f23.15)') type_atom(i), 0, r_atom(:,i) + outn = outn + 1 end do + if((outn-node_atoms) /= atom_num) then + print *, "outn", (outn-node_atoms), " doesn't equal atom_num ", atom_num + end if + !Finish writing close(11) end subroutine write_xyz @@ -161,8 +175,10 @@ module io !Calculate total atom number write_num = atom_num do i = 1,ele_num - if(type_ele(i) == 'fcc') write_num = write_num + basisnum(lat_ele(i))*size_ele(i)**3 + if(type_ele(i) == 'fcc') write_num = write_num + size_ele(i)**3 + if(type_ele(i) == 'bcc') 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 @@ -196,7 +212,7 @@ module io 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') + case('fcc','bcc') do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 interp_num = interp_num+1 call apply_periodic(r_interp(:,iatom)) @@ -206,6 +222,88 @@ module io end do end subroutine write_lmp + subroutine write_ldump(file) + !This subroutine will only work if element data is defined + character(len = *), intent(in) :: file + integer :: write_num, i, iatom + logical :: write_dat + integer :: type_interp(max_basisnum*max_esize**3), interp_num + real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), data_interp(10, max_basisnum*max_esize**3) + + + +1 format('ITEM: TIMESTEP'/i16) +2 format('ITEM: NUMBER OF ATOMS' /i16) +3 format('ITEM: BOX BOUNDS ', 2a1, ' ', 2a1, ' ', 2a1 / & + 2f23.15 / 2f23.15 / 2f23.15) +4 format('ITEM: ATOMS id type x y z energy fx fy fz s11 s22 s33 s23 s13 s12') +5 format('ITEM: ATOMS id type x y z') + + open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') + + print *, max_basisnum, max_esize + !Write header information + write(11,1) timestep + + !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 + if(type_ele(i) == 'bcc') write_num = write_num + size_ele(i)**3 + end do + !Write total number of atoms + write(11,2) write_num + !Write box information + write(11,3) box_bc(1:1), box_bc(1:1), box_bc(2:2), box_bc(2:2), box_bc(3:3), box_bc(3:3), box_bd(:) + + !Now pick if we are interpolating data or not + if(allocated(force_node).or.allocated(force_atom)) then + write(11, 4) + write_dat = .true. + else + write(11, 5) + write_dat = .false. + end if + + if (write_dat) then + do i = 1, atom_num + write(11, '(2i16, 13f23.15)') i, type_atom(i), r_atom(:,i), energy_atom(i), force_atom(:,i), virial_atom(:,i) + end do + else + do i = 1, atom_num + write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i) + end do + end if + + !Write refined element atomic positions + interp_num = 0 + do i = 1, ele_num + if(write_dat) then + call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp, & + energy_node(:,:,i), force_node(:,:,:,i), virial_node(:,:,:,i), data_interp) + else + call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp) + end if + select case(trim(adjustl(type_ele(i)))) + case('fcc','bcc') + if(write_dat) then + do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 + interp_num = interp_num+1 + call apply_periodic(r_interp(:,iatom)) + write(11, '(2i16, 13f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom), & + data_interp(:,iatom) + end do + else + do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 + interp_num = interp_num+1 + call apply_periodic(r_interp(:,iatom)) + write(11, '(2i16, 3f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom) + end do + end if + end select + end do + end subroutine write_ldump + subroutine write_lmpcac(file) !This subroutine writes out a .lmp style dump file character(len=100), intent(in) :: file @@ -339,6 +437,7 @@ module io end do close(11) + !Now we write the vtk file for the elements open(unit=11, file='cg_'//trim(adjustl(file)), action='write', status='replace',position='rewind') write(11,1) write(11,2) @@ -357,6 +456,7 @@ module io write(11,5) ele_num do i = 1, ele_num if(trim(adjustl(type_ele(i))) == 'fcc') write(11, '(i16)') 12 + if(trim(adjustl(type_ele(i))) == 'bcc') write(11, '(i16)') 12 end do write(11,12) ele_num write(11,20) @@ -375,55 +475,32 @@ module io !NOTE: This code doesn't work for arbitrary number of basis atoms per node. It assumes that the !each element has only 1 atom type at the node. character(len=100), intent(in) :: file - integer :: interp_max, i, j, inod, ibasis, ip, unique_index(50), unique_size(50), unique_num, & + integer :: interp_max, i, j, inod, ibasis, ip, unique_index(50), unique_size(50), unique_type(50), unique_num, & etype - real(kind=dp) :: box_vec(3) + real(kind=dp) :: box_vec(3), masses(10) 1 format('time' / i16, f23.15) 2 format('number of elements' / i16) 3 format('number of nodes' / i16) -4 format('element types' / i16) 5 format('number of atoms' / i16) -6 format('number of grains' / i16) 7 format('boundary ' / 3a1) 8 format('box bound' / 6f23.15) 9 format('box length' / 3f23.15) 10 format('box matrix') 11 format(3f23.15) 12 format('coarse-grained domain') -13 format('ie ele_type grain_ele lat_type_ele'/ 'ip ibasis x y z') -14 format('atomistic domain' / 'ia grain_atom type_atom x y z') -15 format('maximum lattice periodicity length' / 3f23.15) -16 format('Number of lattice types and atom types '/ 2i16) -17 format('lattice type IDs') -18 format('lattice types for grains') -19 format('max nodes per element' / i16) +13 format('ie basis_num ng_node esize'/ 'ip ibasis type x y z') +14 format('atomistic domain' / 'ia type_atom x y z') +19 format('max nodes per element and basis per nodes' / 2i16) 20 format('max interpo per element' / i16) 21 format('atom types to elements') open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') + !Below writes the header information for the restart file write(11,1) timestep, total_time write(11,2) ele_num - !Below writes the header information for the restart file - - - !First figure out all of the unique element types - unique_num = 0 - unique_index(:) = 0 - eleloop:do i = 1, ele_num - do j =1 , unique_num - if ( ( size_ele(i) == size_ele( unique_index(j) ) ).and. & - ( lat_ele(i) == lat_ele(unique_index(j)) ) ) then - cycle eleloop - end if - end do - unique_num = unique_num + 1 - unique_index(unique_num) = i - unique_size(unique_num) = size_ele(i) - end do eleloop - !Calculate the max number of atoms per element select case(max_ng_node) case(8) @@ -431,31 +508,19 @@ module io case default interp_max = 0 end select + write(11,20) interp_max + !Write write(11,3) node_num - write(11,19) max_ng_node - write(11,4) unique_num + write(11,19) max_ng_node, max_basisnum write(11,5) atom_num - write(11,6) 1 !Grain_num is ignored - write(11,16) lattice_types, atom_types - write(11,21) do i = 1, atom_types - write(11,*) i, type_to_name(i) + call atommass(type_to_name(i),masses(i)) end do + write(11,*) "masses: " + write(11, *) (masses(i), i = 1, atom_types) write(11,7) box_bc(1:1), box_bc(2:2), box_bc(3:3) - write(11,18) - write(11,'(2i16)') 1,1 !This is another throwaway line that is meaningless - write(11,17) - !This may have to be updated in the future but currently the only 8 node element is fcc - do i = 1, lattice_types - select case(ng_node(i)) - case(8) - write(11, *) i, 'fcc' - end select - end do - write(11,15) 1.0_dp, 1.0_dp, 1.0_dp !Another throwaway line that isn't needed write(11,8) box_bd - write(11,9) box_bd(2)-box_bd(1), box_bd(4) - box_bd(3), box_bd(6)-box_bd(5) write(11,10) !Current boxes are limited to being rectangular do i = 1,3 @@ -463,35 +528,18 @@ module io box_vec(i) = box_bd(2*i) - box_bd(2*i-1) write(11,11) box_vec end do - !We write this as box_mat ori and box_mat current - do i = 1,3 - box_vec(:) = 0.0_dp - box_vec(i) = box_bd(2*i) - box_bd(2*i-1) - write(11,11) box_vec - end do !write the element information if(ele_num > 0) then write(11,12) - - do i = 1, unique_num - write(11,'(3i16)') i, size_ele(unique_index(i))-1, basis_type(1,lat_ele(unique_index(i))) - end do ip = 0 write(11,13) do i = 1, ele_num - !Figure out the ele type - do j = 1, unique_num - if ( unique_size(j) == size_ele(i)) then - etype = j - exit - endif - end do - write(11, '(4i16)') i, etype, 1, basis_type(1,lat_ele(i)) + write(11, '(4i16)') i, basisnum(lat_ele(i)), 2, (size_ele(i)-1) do inod = 1, ng_node(lat_ele(i)) do ibasis = 1, basisnum(lat_ele(i)) ip = ip + 1 - write(11, '(2i16, 3f23.15)') ip, ibasis, r_node(:, ibasis, inod, i) + write(11, '(3i16, 3f23.15)') ip, ibasis, basis_type(ibasis, lat_ele(i)), r_node(:, ibasis, inod, i) end do end do end do @@ -501,7 +549,7 @@ module io if(atom_num /= 0) then write(11,14) do i = 1, atom_num - write(11, '(3i16, 3f23.15)') i, 1, type_atom(i), r_atom(:,i) + write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i) end do end if @@ -543,13 +591,13 @@ module io !Write out atoms first do i = 1, atom_num - write(11,*) i, type_atom(i), sbox_atom(i), r_atom(:,i) + write(11,*) tag_atom(i), type_atom(i), sbox_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), sbox_ele(i), type_ele(i) + write(11, *) tag_ele(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) @@ -579,34 +627,41 @@ module io 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 an existing file to read." + stop 3 + end if - !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 + select case(temp_infile(scan(temp_infile,'.',.true.)+1:)) + case('restart', 'mb', 'cac') + infilenum=infilenum+1 + infiles(infilenum) = temp_infile + case('out') + if(atom_types == 0) then + print *, "Please run -set_types command prior to running code requiring reading in pycac_*.out files" + stop 3 end if - - select case(temp_infile(scan(temp_infile,'.',.true.)+1:)) - case('restart', '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, restart." - read(*,*) temp_infile - + select case(trim(adjustl(mode))) + case('--calc', '--convert','--metric') + infilenum = infilenum+1 + infiles(infilenum) = temp_infile + case default + print *, "Files of type .out cannot be used with mode ", trim(adjustl(mode)) + stop 3 end select - end do + + 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, restart, cac, or out." + stop 3 + end select 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 + !This subroutine reads in file i integer, intent(in) :: i real(kind=dp), dimension(3), intent(in) :: displace @@ -618,9 +673,13 @@ module io call read_mb(infiles(i), displace, temp_box_bd) case('restart') call read_pycac(infiles(i), displace, temp_box_bd) + case('cac') + call read_lmpcac(infiles(i), displace, temp_box_bd) + case('out') + call read_pycac_out(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" + " is not accepted for reading. Please select from: mb,restart,cac,out and try again" stop end select @@ -655,7 +714,7 @@ module io 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 - + call grow_box(temp_box_bd) !Read in the number of sub_boxes and allocate the variables read(11, *) n @@ -738,19 +797,19 @@ module io do i = 1, in_atoms read(11,*) j, type, sbox, r(:) r = r+newdisplace - call add_atom(new_type_to_type(type), sbox+sub_box_num, r) + call add_atom(j, new_type_to_type(type), sbox+sub_box_num, r) end do !Read the elements do i = 1, in_eles - read(11, *) l, type, size, sbox, etype + read(11, *) j, type, size, sbox, etype do inod = 1, ng_node(type) do ibasis =1, basisnum(type) - read(11,*) j, k, r_innode(:, ibasis, inod) + read(11,*) k, l, r_innode(:, ibasis, inod) r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace end do end do - call add_element(etype, size, new_lattice_map(type), sbox+sub_box_num, r_innode) + call add_element(j, etype, size, new_lattice_map(type), sbox+sub_box_num, r_innode) end do !Close the file being read @@ -773,48 +832,42 @@ module io real(kind=dp), dimension(3), intent(in) :: displace real(kind = dp), dimension(6), intent(out) :: temp_box_bd - integer :: i, inod, ibasis, j, k, in_eles, in_atoms, ele_types, in_lat_num, in_atom_types, & - atom_type_map(10), etype_map(10), etype, lat_type, new_lattice_map(10), & - atom_type - real(kind=dp) :: newdisplace(3), r_in(3,1,8), r_in_atom(3), new_displace(3) + integer :: i, inod, ibasis, j, k, l, in_eles, in_atoms, ele_types, in_lat_num, in_atom_types, & + atom_type_map(100), etype_map(100), etype, lat_type, new_lattice_map(100), & + atom_type, stat + real(kind=dp) :: newdisplace(3), r_in(3,1,8), r_in_atom(3), atomic_masses(10) character(len=100) :: textholder, in_lattype_map(10) character(len=2) :: atomic_element !First open the file open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') - !Disregard unneeded information - do i = 1, 3 - read(11,*) textholder - end do + !Read the timestep information + read(11,*) textholder + read(11,*) timestep, total_time !Read element number + read(11,*) textholder read(11,*) in_eles - !Discard info and read ng_max_node - do i = 1, 5 + !Discard info and read ng_max_node and max_basisnum + do i = 1, 5 read(11,*) textholder end do - read(11,*) max_ng_node - - !Read element types (only needed inside this subroutine) - read(11,*) textholder - read(11,*) ele_types + read(11,*) max_ng_node, max_basisnum !Read in atom num read(11,*) textholder read(11,*) in_atoms - !read in lattice_types and atom types - do i = 1,3 - read(11,*) textholder - end do - read(11,*) in_lat_num, in_atom_types - + !read in atom masses + read(11, *) textholder + read(11, '(a)') textholder + j = tok_count(textholder) + read(textholder, *) (atomic_masses(i), i=1, j) - !Read define atom_types by name - read(11,*) textholder + !Read define atom_types by mass do i = 1, in_atom_types - read(11,*) j, atomic_element + call atommassspecies(atomic_masses(i), atomic_element) call add_atom_type(atomic_element, atom_type_map(i)) end do @@ -822,22 +875,6 @@ module io read(11,*) textholder read(11,*) box_bc - !Disregard useless info - do i = 1, 3 - read(11,*) textholder - end do - - !Read in lat_type map - do i = 1, in_lat_num - read(11,*) j, in_lattype_map(i) - ng_node(lattice_types+i) = 8 !Only cubic elements are currently supported in pyCAC - end do - - !Disregard useless info - do i =1 , 3 - read(11,*) textholder - end do - !Read box boundaries and displace them if necessary read(11,*) temp_box_bd(:) do i = 1, 3 @@ -865,64 +902,21 @@ module io sub_box_bd(:, sub_box_num+1) = temp_box_bd !Read in more useless info - do i = 1, 10 + do i = 1, 9 read(11,*) textholder end do - !Now read the ele_type to size and lat map - do i = 1, ele_types - read(11,*) j, etype_map(i) - etype_map(i) = etype_map(i) + 1 - end do - - - !Now set up the lattice types. In this code it assumes that lattice_type 1 maps to - !atom type 1 because it only allows 1 atom per basis in pyCAC at the moment. - do i = 1, in_lat_num - basisnum(lattice_types+i) = 1 - basis_type(1,lattice_types+i) = atom_type_map(i) - end do - - !Figure out the lattice type maps in case we have repeated lattice_types - k = lattice_types + 1 - new_lattice_map(:) = 0 - new_loop:do i = 1, in_lat_num - old_loop:do j = 1, lattice_types - !First check all the lattice level variables - if ((basisnum(lattice_types+i) == basisnum(j)).and. & - (ng_node(lattice_types+i) == ng_node(j))) then - !Now check the basis level variables - do ibasis =1, basisnum(i) - if(basis_type(ibasis,lattice_types+i) /= basis_type(ibasis,j)) then - cycle old_loop - end if - end do - new_lattice_map(i) = j - cycle new_loop - end if - end do old_loop - new_lattice_map(i) = k - k = k+1 - end do new_loop - - !Read more useless data - read(11,*) textholder - - !set max values and allocate variables - max_basisnum = maxval(basisnum) - max_ng_node = maxval(ng_node) - call grow_ele_arrays(in_eles, in_atoms) - !Now start reading the elements if(in_eles > 0) then + read(11,*) textholder read(11,*) textholder do i = 1, in_eles read(11,*) j, etype, k, lat_type do inod = 1, 8 - read(11, *) j, k, r_in(:,1,inod) + read(11, *) k, l, r_in(:,1,inod) r_in(:,1,inod) = r_in(:,1,inod) + newdisplace end do - call add_element(in_lattype_map(lat_type), etype_map(etype), new_lattice_map(lat_type), sub_box_num + 1, r_in) + call add_element(j, in_lattype_map(lat_type), etype_map(etype), new_lattice_map(lat_type), sub_box_num + 1, r_in) end do end if @@ -935,9 +929,13 @@ module io end if do i = 1, in_atoms - read(11,*) j, k, atom_type, r_in_atom(:) + read(11,*, iostat=stat) j, k, atom_type, r_in_atom(:) + if(stat > 0) then + print *, j + stop + end if r_in_atom = r_in_atom + newdisplace - call add_atom(atom_type_map(atom_type), sub_box_num + 1, r_in_atom) + call add_atom(j,atom_type_map(atom_type), sub_box_num + 1, r_in_atom) end do !Close file close(11) @@ -949,4 +947,274 @@ module io call set_max_esize end if end subroutine read_pycac + + subroutine read_pycac_out(file, displace, temp_box_bd) + !This subroutine reads in the pyCAC dump file + + + !Arguments + character(len=100), intent(in) :: file + real(kind=dp), dimension(3), intent(in) :: displace + real(kind=dp), dimension(6), intent(out) :: temp_box_bd + + !Internal Variables + integer :: i, in_eles, in_atoms, inbtypes(10), lat_type, ia, ie, inod, & + id, type_node, ilat, esize, tag, type, bnum, n, ibasis, ip + real(kind=dp) :: newdisplace(3), ra(3), in_lapa, ea, fa(3), va(6), & + ee(10,8), fe(3,10,8), ve(6,10,8), re(3,10,8) + character(len=100) :: textholder, fcc + character(len=1000) :: line + + + open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') + + !Now initialize some important variables if they aren't defined + if (max_basisnum==0) max_basisnum = 10 + if (max_ng_node==0) max_ng_node=8 + fcc="fcc" + + !Skip header comment lines + read(11, *) textholder + read(11, *) textholder + + !Read the timestep + read(11, *) textholder, timestep + + !Read atom number and element number and grow element arrays by needed amount + read(11,*) textholder, in_atoms, textholder, in_eles + call grow_ele_arrays(in_eles, in_atoms) + call alloc_dat_arrays(in_eles, in_atoms) + + !Read boundary information + read(11,*) textholder, box_bc(1:1), box_bc(2:2), box_bc(3:3), temp_box_bd(:) + do i = 1, 3 + if (abs(displace(i)) > lim_zero) then + newdisplace(i) = displace(i) - temp_box_bd(2*i-1) + else + newdisplace(i)=displace(i) + end if + 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 + + call grow_box(temp_box_bd) + + !Allocate sub_box boundaries + if (sub_box_num == 0) then + call alloc_sub_box(1) + else + call grow_sub_box(1) + end if + + !Because orientations and other needed sub_box information isn't really + !present within the .cac file we just default a lot of it. + sub_box_ori(:,:,sub_box_num+1) = identity_mat(3) + sub_box_bd(:, sub_box_num+1) = temp_box_bd + sub_box_num = sub_box_num + 1 + + if(in_atoms > 0 ) then + !Read atom header + read(11,*) textholder + do ia = 1, in_atoms + read(11,'(a)') line(:) + read(line,*) tag, type, ra(:), ea, fa(:), va(:) + call add_atom(tag, type, sub_box_num, ra) + call add_atom_data(atom_num, ea, fa, va) + end do + + end if + + if(in_eles > 0) then + !Read element and node headers + read(11,*) textholder + read(11,*) textholder + + !read element information, currently only 8 node elements with 1 basis + do ie =1, in_eles + read(11,*) tag, n, bnum, esize + inbtypes(:) = 0 + do inod =1, n*bnum + read(11,*) ip, ibasis, inbtypes(ibasis), re(:,ibasis,ip), ee(ibasis,ip), fe(:,ibasis,ip), ve(:,ibasis,ip) + end do + call lattice_map(bnum, inbtypes, n, 1.0_dp, lat_type) + call add_element(tag, fcc, esize+1, lat_type, sub_box_num, re) + call add_element_data(ele_num, ee, fe, ve) + end do + end if + call set_max_esize + return + end subroutine + + + + subroutine read_lmpcac(file, displace, temp_box_bd) + !This subroutine is used to read .cac files which are used with the lammpsCAC format + character(len=100), intent(in) :: file + real(kind=dp), dimension(3), intent(in) :: displace + real(kind = dp), dimension(6), intent(out) :: temp_box_bd + + character(len=100) :: textholder, element_type + character(len=2) :: atom_species + integer :: i, j, k, ele_in, type_in, type_map(10), in_basis, node_types(10,8), inod, ibasis, in_basis_types(10), esize, & + lat_type + real(kind=dp) :: mass, r_in(3,10,8), lat_vec(3,3), in_ori(3,3), newdisplace(3) + + !First check to make sure that we have set the needed variables + if(is_equal(in_lapa,0.0_dp).or.(in_lattice_type=='')) then + print *, "Please use set_cac to set needed parameters to read in .cac file" + stop 3 + end if + !Open the file + open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') + + !Now initialize some important variables if they aren't defined + if (max_basisnum==0) max_basisnum = 10 + if (max_ng_node==0) max_ng_node=8 + + !Read header information + read(11, *) textholder + + !Read number of elements + read(11, *) ele_in, textholder + read(11, *) type_in, textholder + + !Read box_boundaries + read(11,*) temp_box_bd(1:2), textholder + read(11,*) temp_box_bd(3:4), textholder + read(11,*) temp_box_bd(5:6), textholder + + !Shift the box boundaries if needed + do i = 1, 3 + if (abs(displace(i)) > lim_zero) then + newdisplace(i) = displace(i) - temp_box_bd(2*i-1) + else + newdisplace(i)=displace(i) + end if + 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 + + !Grow box boundaries + call grow_box(temp_box_bd) + + !Allocate sub_box + if (sub_box_num == 0) then + call alloc_sub_box(1) + else + call grow_sub_box(1) + end if + + !Because orientations and other needed sub_box information isn't really + !present within the .cac file we just default a lot of it. + sub_box_ori(:,:,sub_box_num+1) = identity_mat(3) + sub_box_bd(:, sub_box_num+1) = temp_box_bd + sub_box_num = sub_box_num + 1 + + + !Read useless information + read(11,*) textholder + + !Read atomic masses + do i = 1, type_in + read(11,*) j, mass, textholder + call ATOMMASSSPECIES(mass, atom_species) + call add_atom_type(atom_species, type_map(i)) + end do + + !Read useless info + read(11,*) textholder + + !Start the reading loop + do i = 1, ele_in + read(11,*) j, element_type, in_basis, esize + !Check to see if we need to grow the max_basis_num + select case(trim(adjustl(element_type))) + case('Eight_Node') + !Read in all the data + do j = 1, 8*in_basis + read(11, *) inod, ibasis, in_basis_types(ibasis), r_in(:,ibasis,inod) + end do + + !Now calculate the lattice vectors and shift the nodal points from the corners to the center of the unit cell + !Please check the nodal numbering figure in the readme in order to understand which nodes are used for the + !calculation + lat_vec(:,1) = (r_in(:,1,2) - r_in(:,1,1))/(2*esize) + lat_vec(:,2) = (r_in(:,1,4) - r_in(:,1,1))/(2*esize) + lat_vec(:,3) = (r_in(:,1,5) - r_in(:,1,1))/(2*esize) + + !Now shift all the nodal positions + select case(trim(adjustl(in_lattice_type))) + case('fcc','FCC') + do ibasis = 1, in_basis + r_in(:,ibasis,1) = r_in(:,ibasis,1) + lat_vec(:,1) + lat_vec(:,2) + lat_vec(:,3) + newdisplace + r_in(:,ibasis,2) = r_in(:,ibasis,2) - lat_vec(:,1) + lat_vec(:,2) + lat_vec(:,3) + newdisplace + r_in(:,ibasis,3) = r_in(:,ibasis,3) - lat_vec(:,1) - lat_vec(:,2) + lat_vec(:,3) + newdisplace + r_in(:,ibasis,4) = r_in(:,ibasis,4) + lat_vec(:,1) - lat_vec(:,2) + lat_vec(:,3) + newdisplace + r_in(:,ibasis,5) = r_in(:,ibasis,5) + lat_vec(:,1) + lat_vec(:,2) - lat_vec(:,3) + newdisplace + r_in(:,ibasis,6) = r_in(:,ibasis,6) - lat_vec(:,1) + lat_vec(:,2) - lat_vec(:,3) + newdisplace + r_in(:,ibasis,7) = r_in(:,ibasis,7) - lat_vec(:,1) - lat_vec(:,2) - lat_vec(:,3) + newdisplace + r_in(:,ibasis,8) = r_in(:,ibasis,8) + lat_vec(:,1) - lat_vec(:,2) - lat_vec(:,3) + newdisplace + end do + case default + print *, in_lattice_type, " is not an accepted lattice type. Please select from: fcc" + end select + !Now map it to either an existing or new lattice type + call lattice_map(in_basis, in_basis_types, 8, in_lapa, lat_type) + + !Now add the element + call add_element(0,in_lattice_type, esize, lat_type, sub_box_num, r_in(:,1:max_basisnum,1:max_ng_node)) + + case('Atom') + read(11, *) inod, ibasis, in_basis_types(ibasis), r_in(:,1,1) + call add_atom(0,in_basis_types(ibasis), sub_box_num, r_in(:,1,1)) + end select + end do + + end subroutine read_lmpcac + + subroutine set_cac(apos) + !This code parses input values + integer, intent(in) :: apos + integer :: arglen, arg_pos + + character(len=100) :: textholder + + arg_pos = apos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) then + print *, "Missing lattice parameter for set_input_lat" + end if + read(textholder,*) in_lapa + print *, in_lapa + + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) then + print *, "Missing lattice type for set_input_lat" + end if + read(textholder,*) in_lattice_type + print *, in_lattice_type + + end subroutine set_cac + + subroutine set_types(apos) + !This code + integer, intent(in) :: apos + integer :: i, j,arglen, arg_pos, ntypes + + character(len=100) :: textholder + + arg_pos = apos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) stop "Missing numtypes in io" + read(textholder,*) ntypes + + do i=1,ntypes + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + call add_atom_type(textholder, j) + end do + + return + end subroutine set_types end module io diff --git a/src/main.f90 b/src/main.f90 index 958a3c0..e508ece 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -15,6 +15,7 @@ program main use parameters use elements use io + use caller integer :: i, end_mode_arg, arg_num, arg_pos @@ -60,6 +61,11 @@ program main !This lets us know if we need to wrap atomic positions back into the cell case('-wrap') wrap_flag=.true. + + case('-set_cac') + call set_cac(i) + case('-set_types') + call set_types(i) end select end do !Determine if a mode is being used and what it is. The first argument has to be the mode @@ -68,7 +74,8 @@ program main argument = trim(adjustl(argument)) if (argument(1:2) == '--') then - call call_mode(end_mode_arg, argument) + mode = argument + call call_mode(end_mode_arg) end if !Check to make sure a mode was called @@ -105,11 +112,13 @@ program main if(bound_called) call def_new_box !Check to make sure a file was passed to be written out and then write out - ! Before building do a check on the file - if (outfilenum == 0) then - argument = 'none' - call get_out_file(argument) + ! Before building do a check on the file + if ((trim(adjustl(mode)) /= "--metric").and.(trim(adjustl(mode)) /= "--calc"))then + if ((outfilenum == 0)) then + argument = 'none' + call get_out_file(argument) + end if + call write_out end if - call write_out end program main diff --git a/src/mode_calc.f90 b/src/mode_calc.f90 new file mode 100644 index 0000000..497c90b --- /dev/null +++ b/src/mode_calc.f90 @@ -0,0 +1,95 @@ +module mode_calc + !This mode is used to calculate various quantities based on input information + use parameters + use io + use subroutines + use elements + use box + + character(len=100) :: calc_opt + real(kind=dp), allocatable :: calculated(:) + public + contains + subroutine calc(arg_pos) + !Main calling subroutine for mode_create + integer, intent(out) :: arg_pos + + print *, '------------------------Mode Calc----------------------------' + + !First parse command + call parse(arg_pos) + + print *, "Calculating ", trim(adjustl(calc_opt)), " for ", ele_num, " elements and ", atom_num, " atoms." + !Now call the correct calc function based on calc_opt + select case(trim(adjustl(calc_opt))) + case('tot_virial') + allocate(calculated(6)) + call calc_tot_virial + case default + print *, trim(adjustl(calc_opt)), " is not accepted as a calc option in mode_calc" + stop 3 + end select + end subroutine calc + + subroutine parse(arg_pos) + !This parses the mode calc options + integer, intent(out) :: arg_pos + + character(len = 100) :: infile + integer:: arglen + real(kind=dp) :: temp_box_bd(6) + + call get_command_argument(2, infile, arglen) + if (arglen == 0 ) stop "Missing calc option in mode calc" + 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) + + call get_command_argument(3, calc_opt, arglen) + if (arglen == 0 ) stop "Missing calc option in mode calc" + + arg_pos = 4 + end subroutine parse + + subroutine calc_tot_virial + !Calculate the the total box pressure in GPa + + integer :: i, j, ibasis, inod + real(kind=dp) :: avg_virial(6) + + !First check to make sure that the virial was set for the atoms/elements + if(.not.vflag) then + print *, "Virial data has not been sent/may not be available with your current input file " + stop 3 + end if + + !Sum the atom virials + calculated = 0 + do i = 1, atom_num + do j = 1, 6 + calculated(j) = calculated(j) + virial_atom(j, i) + end do + end do + + !Sum the nodal virials + do i = 1, ele_num + avg_virial(:) = 0 + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + do j = 1,6 + avg_virial(j) = avg_virial(j) + virial_node(j,ibasis,inod,i)/(basisnum(lat_ele(i))*ng_node(lat_ele(i))) + end do + end do + end do + + !Now add the total virial from the element + calculated = calculated + avg_virial*(esize**3.0_dp) + end do + + !Now calculate the total box virial and convert to GPa + calculated = evtogp(calculated)/box_volume() + + print *, "Total virial is calculated as : (v11, v22, v33, v32, v31, v21)" + print *, calculated + end subroutine +end module mode_calc diff --git a/src/mode_convert.f90 b/src/mode_convert.f90 index 3bb8a3b..eb89994 100644 --- a/src/mode_convert.f90 +++ b/src/mode_convert.f90 @@ -14,6 +14,7 @@ module mode_convert character(len=100) :: infile real(kind = dp) :: temp_box_bd(6) !First read in the file + temp_box_bd(:) = 0.0_dp 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) @@ -21,4 +22,4 @@ module mode_convert arg_pos = 3 end subroutine convert -end module mode_convert \ No newline at end of file +end module mode_convert diff --git a/src/mode_create.f90 b/src/mode_create.f90 index ed60daf..7994f55 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -10,14 +10,15 @@ module mode_create implicit none - character(len=100) :: name, element_type + 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) 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 + basis_pos(3,10), esize_nums, esize_index(10) + logical :: dup_flag, dim_flag, efill real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:) + integer, allocatable :: elat(:) public contains @@ -26,7 +27,7 @@ module mode_create integer, intent(out) :: arg_pos - integer :: i, ibasis, inod + integer :: i, ibasis, inod, ei, curr_esize real(kind=dp), allocatable :: r_node_temp(:,:,:) print *, '-----------------------Mode Create---------------------------' @@ -45,6 +46,7 @@ module mode_create basisnum = 0 lat_ele_num = 0 lat_atom_num = 0 + efill = .false. !First we parse the command call parse_command(arg_pos) @@ -63,10 +65,17 @@ module mode_create do i = 1, 8 box_vert(:,i) = duplicate(:)*esize*lattice_space(:)*cubic_cell(:,i) + (origin(:)/lattice_parameter) end do - !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 + !Now get the rotated box vertex positions in lattice space. Should be integer units and get the new maxlen + select case(trim(adjustl(element_type))) + case('fcc') + box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1 + maxbd = maxval(matmul(orient,matmul(fcc_mat,box_lat_vert)),2) + case('bcc') + box_lat_vert = int(matmul(bcc_inv, matmul(orient_inv, box_vert)))+1 + maxbd = maxval(matmul(orient,matmul(bcc_mat,box_lat_vert)),2) + end select + !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) @@ -83,7 +92,12 @@ module mode_create box_vert(:,i) = (cubic_cell(:,i)*box_len(:) + origin(:))/lattice_parameter end do !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 + select case(trim(adjustl(element_type))) + case('fcc') + box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1 + case('bcc') + box_lat_vert = int(matmul(bcc_inv, matmul(orient_inv, box_vert)))+1 + end select !Now get the box_bd in lattice units do i = 1, 3 @@ -105,10 +119,10 @@ module mode_create end do end do do i = 1,3 - box_bd(2*i) = maxval(r_node_temp(i,:,:))+10.0_dp**-6.0_dp - box_bd(2*i-1) = minval(r_node_temp(i,:,:)) - 10.0_dp**-6.0_dp + box_bd(2*i) = maxval(r_node_temp(i,:,:))+10.0_dp**(-6.0_dp) + box_bd(2*i-1) = minval(r_node_temp(i,:,:)) - 10.0_dp**(-6.0_dp) end do - call add_element(element_type, esize, 1, 1, r_node_temp) + call add_element(0,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 @@ -117,10 +131,11 @@ module mode_create !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('bcc') + call build_with_rhomb(box_lat_vert, bcc_mat) case default - print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", & + print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ",& "element type" stop 3 end select @@ -135,7 +150,7 @@ module mode_create if(lat_atom_num > 0) then do i = 1, lat_atom_num do ibasis = 1, basisnum(1) - call add_atom(basis_type(ibasis, 1), 1, (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis)) + call add_atom(0,basis_type(ibasis, 1), 1, (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis)) end do end do deallocate(r_atom_lat) @@ -148,7 +163,8 @@ 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, 1, r_node_temp) + + call add_element(0,element_type, elat(i), 1, 1, r_node_temp) end do end if end if @@ -248,6 +264,8 @@ module mode_create end do end do + case('efill') + efill=.true. case default !If it isn't an option then you have to exit arg_pos = arg_pos -1 @@ -272,13 +290,20 @@ module mode_create 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 + 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 + case('bcc') + do i = 1, 3 + !Check if the direction is 111 + if((is_equal(abs(orient(i,1)),abs(orient(i,2)))).and.(is_equal(abs(orient(i,2)),abs(orient(i,3))))) then + lattice_space(i) = 0.5_dp * lattice_space(i) end if end do end select @@ -313,14 +338,19 @@ module mode_create 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) + integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,8), rzero(3), efill_size, & + vlat(3), temp_lat(3,8), m, n, o, j, k, nump_ele, efill_temp_lat(3,8), filzero(3), bd_ele_lat(6),& + efill_ele(3,8), ebd(6) + real(kind=dp) :: v(3), temp_nodes(3,1,8), r(3), centroid_bd(6) logical, allocatable :: lat_points(:,:,:) - logical :: node_in_bd(8) + logical :: node_in_bd(8), add, lat_points_ele(esize,esize,esize), intersect_bd(3) !Do some value initialization max_esize = esize + do i = 1,3 + centroid_bd(2*i) = -huge(1.0_dp) + centroid_bd(2*i-1) = huge(1.0_dp) + end do !First find the bounding lattice points (min and max points for the box in each dimension) numlatpoints = 1 @@ -415,8 +445,9 @@ module mode_create !Now build the finite element region lat_ele_num = 0 lat_atom_num = 0 - allocate(r_lat(3,8,numlatpoints/esize)) - + allocate(r_lat(3,8,numlatpoints/esize), elat(numlatpoints/esize)) + + ele(:,:) = (esize-1) * cubic_cell(:,:) !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 @@ -444,17 +475,19 @@ module mode_create !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 + 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 we are on the first round of element building then we can just add the element if all(node_in_bd) is + !true if(all(node_in_bd)) then lat_ele_num = lat_ele_num+1 r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) - + elat(lat_ele_num) = esize !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,:)) @@ -464,11 +497,83 @@ module mode_create end do end do + !If any nodes are in the boundary and we want to efill then run the efill code + else if(any(node_in_bd).and.efill) then + + !Pull out the section of the lat points array + lat_points_ele(:,:,:)=.false. + do i = 1,3 + if (minval(temp_lat(i,:)) size(lat_points,i)) then + bd_ele_lat(2*i) = size(temp_lat(i,:)) + else + bd_ele_lat(2*i) = maxval(temp_lat(i,:)) + end if + end do + + lat_points_ele(1:(bd_ele_lat(2)-bd_ele_lat(1)),1:(bd_ele_lat(4)-bd_ele_lat(3)),& + 1:(bd_ele_lat(6)-bd_ele_lat(5)))= lat_points(bd_ele_lat(1):bd_ele_lat(2), & + bd_ele_lat(3):bd_ele_lat(4), & + bd_ele_lat(5):bd_ele_lat(6)) + !Now start looping through elements and try to fit as many as you can + efill_size = esize-2 + i=1 + j=1 + k=1 + nump_ele = count(lat_points_ele) + do i = 1, 3 + filzero(i) = bd_ele_lat(2*i-1) -1 + end do + do while(efill_size>min_efillsize) + !First check whether there are enough lattice points to house the current element size + efill_ele=cubic_cell*(efill_size-1) + if (nump_ele < efill_size**3) then + efill_size = efill_size - 2 + else + ze: do k = 1, (esize-efill_size) + ye: do j = 1, (esize-efill_size) + xe: do i = 1, (esize-efill_size) + do inod = 1,8 + vlat = efill_ele(:,inod) + (/ i, j, k /) + if (.not.lat_points_ele(vlat(1),vlat(2),vlat(3))) cycle xe + do o = 1,3 + v(o) = real(vlat(o) + filzero(o) + bd_in_lat(2*o-1) -5) + end do + temp_nodes(:,1, inod) = matmul(orient, matmul(transform_matrix, v)) + efill_temp_lat(:,inod) = vlat + end do + + do o = 1,3 + ebd(2*o-1) = minval(efill_temp_lat(o,:)) + ebd(2*o) = maxval(efill_temp_lat(o,:)) + end do + lat_ele_num = lat_ele_num+1 + r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) + elat(lat_ele_num) = efill_size + nump_ele = nump_ele - efill_size**3 + !Now set all the lattice points contained within an element to false + do o = ebd(5), ebd(6) + do n = ebd(3), ebd(4) + do m = ebd(1), ebd(2) + lat_points(m+filzero(1),n+filzero(2),o+filzero(3)) = .false. + lat_points_ele(m,n,o) = .false. + end do + end do + end do + end do xe + end do ye + end do ze + efill_size = efill_size-2 + end if + 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 @@ -518,6 +623,5 @@ module mode_create STOP 3 end subroutine error_message - - + end module mode_create diff --git a/src/mode_merge.f90 b/src/mode_merge.f90 index afc9bba..242ecd2 100644 --- a/src/mode_merge.f90 +++ b/src/mode_merge.f90 @@ -25,6 +25,7 @@ module mode_merge shift_flag = .false. shift_vec(:) = 0.0_dp + temp_box_bd(:) = 0.0_dp !First we parse the merge command call parse_command(arg_pos) @@ -41,7 +42,6 @@ module mode_merge 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') @@ -53,7 +53,6 @@ module mode_merge end select call read_in(i, displace, temp_box_bd) - call grow_box(temp_box_bd) end if if(shift_flag) call shift(new_starts, i) @@ -168,4 +167,4 @@ module mode_merge end if end subroutine shift -end module mode_merge \ No newline at end of file +end module mode_merge diff --git a/src/mode_metric.f90 b/src/mode_metric.f90 new file mode 100644 index 0000000..dbe910e --- /dev/null +++ b/src/mode_metric.f90 @@ -0,0 +1,249 @@ +module mode_metric + !This mode is used to calculate continuum metrics for the j + + use parameters + use io + use elements + use neighbors + + implicit none + + integer :: nfiles + character(len=100) :: metric_type + real(kind=dp) :: rc_off + + !Save reference positions + integer :: np, npreal, nmet + real(kind=dp), allocatable :: r_zero(:,:), r_curr(:,:), met(:,:) + + public + contains + subroutine metric(arg_pos) + !This is the main calling subroutine for the metric code + integer, intent(out) :: arg_pos + character(len=100) :: infile, outfile + + integer :: i, ibasis, inod, np_temp, ppos + real(kind=dp), dimension(6) :: temp_box_bd + + !These are the variables containing the cell list information + integer, dimension(3) :: cell_num + integer, allocatable :: num_in_cell(:,:,:), which_cell(:,:) + integer, allocatable :: cell_list(:,:,:,:) + + !Parse the command arguments + call parse_command(arg_pos) + + !Now read the first file + call read_in(1, (/ 0.0_dp, 0.0_dp, 0.0_dp /), temp_box_bd) + np = atom_num + max_basisnum*max_ng_node*ele_num + allocate(r_zero(3,atom_num+max_basisnum*max_ng_node*ele_num), & + r_curr(3,atom_num+max_basisnum*max_ng_node*ele_num)) + r_zero(:,:) = -huge(1.0_dp) + + !Set up the met variable for the user desired metric + select case(trim(adjustl(metric_type))) + case('def_grad') + allocate(met(9, np)) + case('microrotation') + allocate(met(4,np)) + end select + + !Now set the reference positions + call convert_positions(r_zero, npreal) + + !Now calculate the neighbor list for the reference configuration + call calc_neighbor(5.0_dp, r_zero, np) + + !Reset element and box + call reset_data + call reset_box + + !Now loop over new files + do i = 2, nfiles + call read_in(i, (/ 0.0_dp, 0.0_dp, 0.0_dp /), temp_box_bd) + call convert_positions(r_curr, np_temp) + if (npreal /= np_temp) then + print *, "Error in mode_metric where number of points in ", i, "th file is ", np_temp, " and number of points in" & + // "reference file is", npreal + end if + call calc_metric + !Now create the output file num and write out to xyz format + ppos = scan(trim(infiles(i)),".", BACK= .true.) + if ( ppos > 0 ) then + outfile = infiles(i)(1:ppos)//'xyz' + else + outfile = infiles(i)//'.xyz' + end if + call write_metric_xyz(outfile) + call reset_data + call reset_box + end do + end subroutine metric + + subroutine parse_command(arg_pos) + !This subroutine parses the arguments for mode metric + integer, intent(out) :: arg_pos + + integer :: i, arglen + character(len=100) :: textholder + logical :: file_exists + + !First read the metric to be used + call get_command_argument(2,metric_type,arglen) + if (arglen == 0) stop "Incomplete mode metric command, check documentation" + select case(trim(adjustl(metric_type))) + case("microrotation", "def_grad") + continue + case default + print *, "Mode metric does not accept metric ", metric_type, ". Please select from: microrotation, def_grad" + stop 3 + end select + + !Now read the cutoff radius + call get_command_argument(3,textholder,arglen) + if (arglen == 0) stop "Incomplete mode metric command, check documentation" + read(textholder, *) rc_off + + !Now read the number of files to read and allocate the variables + call get_command_argument(4, textholder, arglen) + if (arglen == 0) stop "Incomplete mode metric command, check documentation" + read(textholder, *) nfiles + + !Now read the files to be read + do i = 1, nfiles + call get_command_argument(4+i, textholder, arglen) + call get_in_file(textholder) + end do + + arg_pos = 5+nfiles + return + end subroutine parse_command + + subroutine calc_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) + + !Loop over all points + do ip = 1, np + eta(:,:) = 0.0_dp + omega(:,:) = 0.0_dp + def_grad(:,:) = 0.0_dp + 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) + oldrij = r_zero(:,nei) - r_zero(:,ip) + + !Calculate eta and omega + do i = 1,3 + do j = 1,3 + omega(i,j) = omega(i,j) + rij(i) * oldrij(j) + eta(i,j) = eta(i,j) + oldrij(i) * oldrij(j) + end do + end do + end do + + eta_inv=matinv3(eta) + def_grad=matmul(omega,eta_inv) + + select case(trim(adjustl(metric_type))) + case('def_grad') + k = 1 + do i = 1,3 + do j = 1, 3 + met(k, ip) = def_grad(i,j) + k = k + 1 + end do + end do + case('microrotation') + met(:,ip) = 0.0_dp + if(.not.all(def_grad == 0)) then + !Now calculate microrotation + ftf = matmul(transpose3(def_grad), def_grad) + U = sqrt3(ftf) + if(.not.all(abs(U) < lim_zero)) then + R = matmul(def_grad, matinv3(U)) + Rskew = 0.5_dp * ( R - transpose3(R)) + do k =1,3 + do j = 1,3 + do i = 1,3 + met(k,ip) = met(k,ip) -0.5*permutation(i,j,k)*Rskew(i,j) + end do + end do + end do + met(4,ip) = norm2(met(1:3,ip)) + end if + end if + end select + end do + return + end subroutine + + subroutine convert_positions(rout, npoints) + !This subroutine just converts current atom and element arrays to a single point based form + real(kind=dp), dimension(3,atom_num+max_ng_node*max_basisnum*ele_num), intent(inout) :: rout + integer, intent(out) :: npoints + + integer :: i, inod, ibasis + + npoints=0 + + if(atom_num > 0) then + do i = 1, atom_num + rout(:,tag_atom(i)) = r_atom(:,i) + npoints = npoints + 1 + end do + end if + + if (ele_num > 0) then + do i = 1, ele_num + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + rout(:, atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis) & + = r_node(:,ibasis,inod,i) + npoints = npoints + 1 + end do + end do + end do + end if + + end subroutine convert_positions + + subroutine write_metric_xyz(outfile) + character(len=100), intent(in) :: outfile + + integer :: i, inod, ibasis + real(kind = dp) :: r(3), eng + open (unit=11, file=trim(adjustl(outfile)), action='write', position='rewind', status = 'replace') + !Write the header + write(11,*) npreal + + select case(metric_type) + case('def_grad') + write(11,*) "type element x y z F11 F12 F13 F21 F22 F23 F31 F32 F33" + case('microrotation') + write(11,*) "type element x y z micro1 micro2 micro3 norm2(micro)" + end select + + if(atom_num > 0) then + do i = 1, atom_num + write(11,*) type_atom(i), 0, r_atom(:,i), met(:,tag_atom(i)) + end do + end if + + if (ele_num > 0) then + do i = 1, ele_num + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + write(11,*) basis_type(ibasis,lat_ele(i)), 1, r_node(:,ibasis,inod,i), & + met(:, atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis) + end do + end do + end do + end if + end subroutine write_metric_xyz + +end module mode_metric diff --git a/src/neighbors.f90 b/src/neighbors.f90 new file mode 100644 index 0000000..0faebca --- /dev/null +++ b/src/neighbors.f90 @@ -0,0 +1,142 @@ +module neighbors + + use parameters + use elements + use subroutines + use functions + + integer, allocatable :: nei_list(:,:), nei_num(:) + 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) + !This subroutine builds a cell list based on rc_off + + !----------------------------------------Input/output variables------------------------------------------- + + integer, intent(in) :: numinlist !The number of points within r_list + + real(kind=dp), dimension(3,numinlist), intent(in) :: r_list !List of points to be used for the construction of + !the cell list. + + real(kind=dp), intent(in) :: rc_off ! Cutoff radius which dictates the size of the cells + + integer, dimension(3), intent(inout) :: cell_num !Number of cells in each dimension. + + integer, allocatable, intent(inout) :: num_in_cell(:,:,:) !Number of points within each cell + + integer, allocatable, intent(inout) :: cell_list(:,:,:,:) !Index of points from r_list within each cell. + + integer, dimension(3,numinlist), intent(out) :: which_cell !The cell index for each point in r_list + + !----------------------------------------Begin Subroutine ------------------------------------------- + + integer :: i, j, cell_lim, c(3) + real(kind=dp) :: box_len(3) + integer, allocatable :: resize_cell_list(:,:,:,:) + + !First calculate the number of cells that we need in each dimension + do i = 1,3 + box_len(i) = box_bd(2*i) - box_bd(2*i-1) + cell_num(i) = int(box_len(i)/(rc_off/2))+1 + end do + + !Initialize/allocate variables + cell_lim = 10 + allocate(num_in_cell(cell_num(1),cell_num(2),cell_num(3)), cell_list(cell_lim, cell_num(1), cell_num(2), cell_num(3))) + + !Now place points within cell + do i = 1, numinlist + !Check to see if the current point is a filler point and if so just skip it + if(r_list(1,i) < -huge(1.0_dp)+1) cycle + + !c is the position of the cell that the point belongs to + do j = 1, 3 + c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off/2)) + 1 + end do + + !Place the index in the correct position, growing if necessary + num_in_cell(c(1),c(2),c(3)) = num_in_cell(c(1),c(2),c(3)) + 1 + if (num_in_cell(c(1),c(2),c(3)) > cell_lim) then + allocate(resize_cell_list(cell_lim+10,cell_num(1),cell_num(2),cell_num(3))) + resize_cell_list(1:cell_lim, :, :, :) = cell_list + resize_cell_list(cell_lim+1:, :, :, :) = 0 + call move_alloc(resize_cell_list, cell_list) + end if + + cell_list(num_in_cell(c(1),c(2),c(3)),c(1),c(2),c(3)) = i + which_cell(:,i) = c + end do + + return + end subroutine build_cell_list + + subroutine calc_neighbor(rc_off, r_list, n) + !This function populates the neighbor list in this module + + real(kind=dp), intent(in) :: rc_off + integer, intent(in) :: n + real(kind=dp), dimension(3,n) :: r_list + + integer :: i, c(3), ci, cj, ck, num_nei, nei + !Variables for cell list code + integer, dimension(3) ::cell_num + integer, allocatable :: num_in_cell(:,:,:), cell_list(:,:,:,:) + integer :: which_cell(3,n) + + !First reallocate the neighbor list codes + if (allocated(nei_list)) then + deallocate(nei_list,nei_num) + end if + + allocate(nei_list(100,n),nei_num(n)) + + !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) + + !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) < -Huge(-1.0_dp)+1) cycle + + !c is the position of the cell that the point + 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 + 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) + nei = cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci) + + !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 + if((nei /= i)) then + + !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 + + nei_num(i) = nei_num(i) + 1 + nei_list(nei_num(i), i) = nei + + end if + end if + end do + end do + end do + end do + + end do pointloop + + return + end subroutine calc_neighbor + +end module neighbors diff --git a/src/opt_deform.f90 b/src/opt_deform.f90 new file mode 100644 index 0000000..c9be654 --- /dev/null +++ b/src/opt_deform.f90 @@ -0,0 +1,98 @@ +module opt_deform + !This module constains the deform option which applies a uniaxial strain to the system + use parameters + use subroutines + use elements + use box + + implicit none + + + real(kind=dp), save :: applied_strain + integer, save :: sdim + + public + contains + + subroutine deform(arg_pos) + !This subroutine applies the simulation box deformation to the system + + integer, intent(inout) :: arg_pos + + character(len=1) :: dims(3) + integer :: i, j, k + real(kind=dp) :: frac_atom(atom_num), frac_node(max_basisnum, max_ng_node, ele_num) + + !initialize some variables + dims(1) = 'x' + dims(2) = 'y' + dims(3) = 'z' + + !First parse the input command + call parse_deform(arg_pos) + print *, '-----------------------Option Deform------------------------' + !Now convert all positions in the specified dimension to fractional coordinates + do i = 1, atom_num + frac_atom(i) = (r_atom(sdim, i) - box_bd(2*sdim-1))/(box_bd(2*sdim)-box_bd(2*sdim-1)) + end do + do i = 1, ele_num + do j = 1, ng_node(lat_ele(i)) + do k = 1, basisnum(lat_ele(i)) + frac_node(k,j,i) = (r_node(sdim,k,j,i) - box_bd(2*sdim-1))/(box_bd(2*sdim)-box_bd(2*sdim-1)) + end do + end do + end do + + print *, "Original box bounds in ", dims(sdim), " are ", box_bd(2*sdim-1:2*sdim) + box_bd(2*sdim) = box_bd(2*sdim) + applied_strain + print *, "New box bounds are ", box_bd(2*sdim-1:2*sdim) + + !Now reassign the positions + do i = 1, atom_num + r_atom(sdim,i) = frac_atom(i)*(box_bd(2*sdim)-box_bd(2*sdim-1)) + box_bd(2*sdim-1) + end do + do i = 1, ele_num + do j = 1, ng_node(lat_ele(i)) + do k = 1, basisnum(lat_ele(i)) + r_node(sdim,k,j,i) = frac_node(k,j,i)*(box_bd(2*sdim)-box_bd(2*sdim-1)) + box_bd(2*sdim-1) + end do + end do + end do + + end subroutine deform + + subroutine parse_deform(arg_pos) + + integer, intent(inout) :: arg_pos + integer :: arg_len + character(len=100) :: textholder + + !Pull out the dimension to be strained + arg_pos=arg_pos+1 + call get_command_argument(arg_pos, textholder, arg_len) + if (arg_len == 0) stop "Missing dim in deform command" + + select case(trim(adjustl(textholder))) + case('x','X') + sdim = 1 + case('y','Y') + sdim = 2 + case('z','Z') + sdim = 3 + case default + print *, "Dimension ", trim(adjustl(textholder)), " is not accepted. Please select either x, y, or z" + end select + + !Now pull out the strain vector, currently only accepts a real number by which to reduce the simulation cell size by in + !that dimension + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arg_len) + if (arg_len == 0) stop "Missing strain in deform command" + read(textholder, *) applied_strain + + arg_pos = arg_pos + 1 + + end subroutine parse_deform + + +end module opt_deform diff --git a/src/opt_delete.f90 b/src/opt_delete.f90 index 09cefb1..c3bcc9e 100644 --- a/src/opt_delete.f90 +++ b/src/opt_delete.f90 @@ -3,6 +3,7 @@ module opt_delete use parameters use subroutines use elements + use neighbors implicit none diff --git a/src/opt_disl.f90 b/src/opt_disl.f90 index 19f8ed8..1833daa 100644 --- a/src/opt_disl.f90 +++ b/src/opt_disl.f90 @@ -9,6 +9,7 @@ module opt_disl implicit none integer :: vloop_size ! Number of atoms to remove in planar vacancy cluster + integer :: iline, islip_plane real(kind=dp), dimension(3) :: line, slip_plane, centroid!dislocation line, slip plane vectors, centroid, real(kind=dp) :: burgers(3) !burgers vector of loop real(kind=dp) :: poisson, char_angle, lattice_parameter!Poisson ratio and character angle, lattice_parameter for burgers vector @@ -31,6 +32,9 @@ module opt_disl print *, '--------------------Option Dislocation-----------------------' select case(trim(adjustl(option))) + case('-disl') + call parse_disl(arg_pos) + call disl case('-dislgen') call parse_dislgen(arg_pos) call dislgen @@ -43,6 +47,146 @@ module opt_disl end select end subroutine dislocation + subroutine parse_disl(arg_pos) + + !Parse disl command + + integer, intent(inout) :: arg_pos + + integer :: i,arglen + character(len=100) :: textholder + character(len=1) :: parse_dim + + !Parse all of the commands + arg_pos = arg_pos + 1 + line(:) = 0.0_dp + + call get_command_argument(arg_pos, parse_dim, arglen) + if (arglen==0) STOP "Missing line dim in disl command" + select case(parse_dim) + case('x','X') + iline=1 + case('y','Y') + iline=2 + case('z','Z') + iline=3 + end select + + arg_pos=arg_pos + 1 + call get_command_argument(arg_pos, parse_dim, arglen) + if (arglen==0) STOP "Missing plane dim in disl command" + select case(parse_dim) + case('x','X') + islip_plane=1 + case('y','Y') + islip_plane=2 + case('z','Z') + islip_plane=3 + end select + do i = 1, 3 + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing centroid in disl command" + call parse_pos(i, textholder, centroid(i)) + end do + + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing parameter in disl command" + read(textholder, *) b + + + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing character angle in disl command" + read(textholder, *) char_angle + + + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing poisson in disl command" + read(textholder, *) poisson + + arg_pos = arg_pos + 1 + end subroutine parse_disl + + subroutine disl + !This subroutine creates the actual dislocation + + integer :: i, sub_box, inod, ibasis, ix, iy, iz + real(kind=dp) :: ss_ori(3,3), ss_inv(3,3), be, bs, slipx(3), disp_transform(3,3), inv_transform(3,3), & + actan, r(3), disp(3) + + print *, "Dislocation with centroid ", centroid, " is inserted" + + !Calculate screw and edge burgers vectors + be = sin(char_angle*pi/180.0_dp)*b + bs = cos(char_angle*pi/180.0_dp)*b + + !Map box dimensions to dislocation dimensions, iz is along the dislocation line and iy is along the slip plane + iz = iline + iy = islip_plane + do i = 1, 3 + if ((i /= iy).and.(i /=iz)) then + ix = i + exit + end if + end do + + if(atom_num > 0) then + do i = 1, atom_num + r=r_atom(:,i) - centroid + if (r(ix) == 0) then + actan=pi/2 + else + actan = datan2(r(iy),r(ix)) + end if + + if ((r(ix)**2 + r(iy)**2) == 0) cycle + + !This is the elastic displacement field for dislocation according to Hirth and Lowe + disp(ix) = be/(2.0_dp*pi) * (actan + (r(ix)*r(iy))/(2.0_dp*(1.0_dp-poisson)*(r(ix)**2.0_dp + r(iy)**2.0_dp))) + disp(iy) = -be/(2.0_dp*pi)*((1.0_dp-2.0_dp*poisson)/(4.0_dp-4.0_dp*poisson) * & + log(r(ix)**2.0_dp + r(iy)**2.0_dp) & + + (r(ix)**2.0_dp - r(iy)**2.0_dp)/(4.0_dp*(1.0_dp-poisson)& + *(r(ix)**2.0_dp+r(iy)**2.0_dp))) + disp(iz) = bs/(2.0_dp*pi) * actan + + r_atom(:,i) = r_atom(:,i) + disp + end do + end if + + if(ele_num > 0) then + do i = 1, ele_num + do inod=1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + r = r_node(:,ibasis,inod,i)-centroid + if (r(ix) == 0) then + actan = pi/2 + else + actan = datan2(r(iy),r(ix)) + end if + + if ((r(ix)**2 + r(iy)**2) == 0) cycle + + !This is the elastic displacement field for dislocation according to Hirth and Lowe + disp(ix) = be/(2.0_dp*pi)*(actan + (r(ix)*r(iy))/(2.0_dp*(1.0_dp-poisson)*(r(ix)**2.0_dp + r(iy)**2.0_dp))) + disp(iy) = -be/(2.0_dp*pi)*((1.0_dp-2.0_dp*poisson)/(4.0_dp-4.0_dp*poisson) * & + log(r(ix)**2.0_dp + r(iy)**2.0_dp) & + + (r(ix)**2.0_dp - r(iy)**2.0_dp)/(4.0_dp*(1.0_dp-poisson)& + *(r(ix)**2.0_dp+r(iy)**2.0_dp))) + disp(iz) = bs/(2.0_dp*pi) * actan + + + r_node(:,ibasis,inod,i) = r_node(:,ibasis,inod,i) + disp + end do + end do + end do + end if + + + end subroutine disl + subroutine parse_dislgen(arg_pos) !Parse dislgen command @@ -123,6 +267,7 @@ module opt_disl call matrix_inverse(ss_ori, 3, ss_inv) !Apply the rotation + print *, ss_inv disp_transform = matmul(sub_box_ori(:,:,i), ss_inv) call matrix_inverse(disp_transform,3,inv_transform) @@ -183,8 +328,6 @@ module opt_disl end do end if - !Now make sure all atoms are wrapped back into the simulation cell - call wrap_atoms end subroutine dislgen @@ -398,8 +541,6 @@ module opt_disl end do return - !Now make sure all atoms are wrapped back into the simulation cell - call wrap_atoms end subroutine !******************************************************** @@ -561,7 +702,6 @@ module opt_disl !Now reset the list for the scanning algorithm delete_num = 0 - !Now scan over all atoms again and find the closest vloop_size number of atoms to the initial atom !that reside on the same plane. If loop_radius is > 0 then we build a circular vacancy cluster if(loop_radius > 0) then do i = 1, atom_num diff --git a/src/opt_group.f90 b/src/opt_group.f90 index dd5a1fb..603bf6e 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -8,10 +8,10 @@ module opt_group use box implicit none - integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2 + integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape - real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), tip_radius, bwidth - logical :: displace, delete, max_remesh, refine + real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth, shell_thickness + logical :: displace, delete, max_remesh, refine, group_nodes, flip, efill, refinefill integer, allocatable :: element_index(:), atom_index(:) @@ -22,31 +22,58 @@ module opt_group !Main calling function for the group option integer, intent(inout) :: arg_pos - print *, '-----------------------Option Group-------------------------' + print *, '------------------------------------------------------------' + print *, 'Option Group' + print *, '------------------------------------------------------------' group_ele_num = 0 group_atom_num = 0 remesh_size=0 + random_num=0 + group_type=0 + notsize=0 displace=.false. delete=.false. max_remesh=.false. refine = .false. + flip = .false. if(allocated(element_index)) deallocate(element_index) if(allocated(atom_index)) deallocate(atom_index) call parse_group(arg_pos) - call get_group !Now call the transformation functions for the group - if(displace) call displace_group + if(refine) then + call get_group + call refine_group + end if + + if(refinefill) then + call get_group + call refinefill_group + end if + + if(displace)then + call get_group + call displace_group + end if - if(remesh_size > 0) call remesh_group + if(delete)then + call get_group + call delete_group + end if - if(delete) call delete_group + if(remesh_size > 0)then + call get_group + call remesh_group + end if - if(refine) call refine_group + if(group_type > 0) then + call get_group + call change_group_type + end if end subroutine group @@ -67,7 +94,7 @@ module opt_group continue case default print *, "Select_type ", trim(adjustl(type)), " is not an accept group selection criteria. ", & - "Please select from atoms, nodes, or both." + "Please select from atoms, elements, or both." end select arg_pos = arg_pos + 1 @@ -158,7 +185,7 @@ module opt_group arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing tip radius in group notch command" - read(textholder,*) tip_radius + read(textholder,*) radius !Calculate the vertex positions vertices(:,1) = centroid @@ -188,7 +215,7 @@ module opt_group !Now update the centroid so that the desire tip diameter matches the wedge with if (bwidth > 0) then - centroid(dim1) = centroid(dim1) + 2*tip_radius*(H)/bwidth + centroid(dim1) = centroid(dim1) + 2*radius*(H)/bwidth end if !Read the ID type shape for group case('id') @@ -280,7 +307,7 @@ module opt_group case('elements','element') if (group_ele_num > 0) then - print *, "Elements specifier used more than once in group id command with type both, either use type ", & + print *, "Elements specifier used more than once in group id command with type both, either use type ",& "elements or include atoms identifier" stop 3 @@ -300,6 +327,47 @@ module opt_group if(i ==1) arg_pos = arg_pos + 1 end do end select + + case('sphere') + !First extract the sphere centroid + do i = 1, 3 + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing sphere centroid in group command" + call parse_pos(i, textholder, centroid(i)) + end do + !Now get the tip radius + arg_pos=arg_pos+1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing sphere radius in group command" + read(textholder, *) radius + + case('shell') + + !First extract the shell centroid + do i = 1, 3 + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing shell centroid in group command" + call parse_pos(i, textholder, centroid(i)) + end do + + !Now get the radius + arg_pos=arg_pos+1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing shell radius in group command" + read(textholder, *) radius + + !Now get the shell thickness + arg_pos=arg_pos+1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing shell thickness in group command" + read(textholder, *) shell_thickness + + case('all') + !Do nothing if the shape is all + continue + case default print *, "Group shape ", trim(adjustl(shape)), " is not currently accepted. Please check documentation ", & "for accepted group shapes." @@ -325,6 +393,8 @@ module opt_group end do case('refine') refine=.true. + case('refinefill') + refinefill=.true. case('remesh') arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) @@ -334,6 +404,28 @@ module opt_group max_remesh =.true. case('delete') delete=.true. + case('nodes') + group_nodes=.true. + case('random') + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) stop "Missing number of random atoms in group command" + read(textholder, *) random_num + case('flip') + flip=.true. + case('efill') + efill=.true. + case('type') + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) stop "Missing atom type for group" + call add_atom_type(textholder, group_type) + case('notsize') + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if(arglen ==0) stop "Missing notsize size" + read(textholder, *) notsize + print *, "Ignoring elements with size ", notsize case default !If it isn't an available option to opt_disl then we just exit exit @@ -344,34 +436,38 @@ module opt_group subroutine get_group !This subroutine finds all elements and/or atoms within the group boundaries !specified by the user. - integer :: i, j, inod, ibasis + integer :: i, j, inod, ibasis, temp, node_in_out(max_ng_node) integer, allocatable :: resize_array(:) - real(kind=dp) :: r_center(3) + real(kind=dp) :: r_center(3), rand select case(trim(adjustl(shape))) case('block') print *, "Group has block shape with boundaries: ", block_bd case ('wedge') print *, "Group has wedge shape with dim1", dim1, "and dim2", dim2, "and vertices ", vertices + case ('notch') + print *, "Group has notch shape with dim1", dim1, "and dim2", dim2, " tip radius ", radius, "and vertices ", & + vertices case('id') print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms." return + case('sphere') + print *, "Group has sphere shape with centroid ", centroid(:), " and radius ", radius + case('all') + print *, "Group has all of type ", type end select - !Allocate variables to arbitrary size + !Reset group if needed + if(allocated(element_index)) deallocate(element_index,atom_index) + + group_ele_num = 0 + group_atom_num = 0 allocate(element_index(1024), atom_index(1024)) !Check the type to see whether we need to find the elements within the group select case(trim(adjustl(type))) case('elements', 'both') do i = 1, ele_num - r_center(:) = 0.0_dp - do inod = 1, ng_node(lat_ele(i)) - do ibasis = 1, basisnum(lat_ele(i)) - r_center = r_center + r_node(:,ibasis,inod,i)/(basisnum(lat_ele(i))*ng_node(lat_ele(i))) - end do - end do - - if (in_group(r_center)) then + if(in_group_ele(size_ele(i), lat_ele(i), r_node(:,:,:,i))) then group_ele_num = group_ele_num + 1 if(group_ele_num > size(element_index)) then allocate(resize_array(size(element_index) + 1024)) @@ -381,15 +477,29 @@ module opt_group end if element_index(group_ele_num) = i - end if - end do - end select + end if + end do - !Check the type to see if we need to find the atoms within the group + if(random_num > 0) then + !If we have the random option enabled then we select random_num number of elements from the group and overwrite + !the group with those elements + do i = 1, random_num + call random_number(rand) + j = i + floor((group_ele_num+1-i)*rand) + temp = element_index(j) + element_index(j) = element_index(i) + element_index(i) = temp + end do + + group_ele_num = random_num + end if + + end select + !Check the type to see if we need to find the atoms within the group select case(trim(adjustl(type))) case('atoms','both') do i = 1, atom_num - if(in_group(r_atom(:,i))) then + if(in_group(r_atom(:,i)).neqv.flip) then group_atom_num = group_atom_num + 1 if (group_atom_num > size(atom_index)) then allocate(resize_array(size(atom_index) + 1024)) @@ -401,6 +511,19 @@ module opt_group atom_index(group_atom_num) = i end if end do + if(random_num > 0) then + !If we have the random option enabled then we select random_num number of atom from the group and overwrite + !the group with those atom + do i = 1, random_num + call random_number(rand) + j = i + floor((group_atom_num+1-i)*rand) + temp = atom_index(j) + atom_index(j) = atom_index(i) + atom_index(i) = temp + end do + + group_atom_num = random_num + end if end select j = 0 @@ -450,7 +573,7 @@ module opt_group end subroutine displace_group subroutine refine_group - !This command is used to remesh the group to a desired element size + !This command is used to refine the group to full atomistics integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3) @@ -469,7 +592,7 @@ module opt_group !here as well to make sure they are in the box do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3 call apply_periodic(r_interp(:,j)) - call add_atom(type_interp(j), sbox_ele(ie), r_interp(:,j)) + call add_atom(0,type_interp(j), sbox_ele(ie), r_interp(:,j)) end do end do !Once all atoms are added we delete all of the elements @@ -477,7 +600,110 @@ module opt_group print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms." end if - end subroutine + end subroutine refine_group + + subroutine refinefill_group + !This command is used to refine the group to full atomistics + + integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num, m, n, o, esize, & + ele(3,8), new_ele_num, ibasis, inod, vlat(3), nump_ele, added_points + real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum,max_ng_node), ravg(3), ratom(3,max_basisnum) + logical :: lat_points(max_esize, max_esize, max_esize) + + !Refining to atoms + if(group_ele_num > 0) then + orig_atom_num = atom_num + new_ele_num = 0 + !Estimate number of atoms we are adding, this doesn't have to be exact + add_atom_num = group_ele_num*basisnum(lat_ele(element_index(1)))*size_ele(element_index(1))**3 + call grow_ele_arrays(0,add_atom_num) + do i = 1, group_ele_num + ie = element_index(i) + !Find all possible elements that we can make while making sure they aren't in the group + lat_points(1:size_ele(ie),1:size_ele(ie),1:size_ele(ie)) = .true. + + !Now calculate the number of elemets which are available for remeshing + nump_ele = size_ele(ie)**3 + do o =1, size_ele(ie) + do n = 1, size_ele(ie) + do m =1, size_ele(ie) + call get_interp_pos(m,n,o,i,rfill(:,:,1)) + + ravg(:) = 0 + do ibasis = 1, basisnum(lat_ele(ie)) + ravg = ravg + rfill(:,ibasis, 1)/basisnum(lat_ele(ie)) + end do + + if( in_group(ravg)) then + nump_ele = nump_ele - 1 + end if + end do + end do + end do + + !Now start the remeshing loop for the element + esize = size_ele(ie) - 2 + added_points=0 + do while(esize > min_efillsize) + if(nump_ele < min_efillsize**3) then + exit + else if (nump_ele < esize**3) then + esize = esize - 2 + else + ele = cubic_cell*(esize-1) + do o = 1, size_ele(ie) - esize + do n = 1, size_ele(ie) - esize + latloop:do m = 1, size_ele(ie) - esize + do inod = 1, ng_node(lat_ele(ie)) + vlat = ele(:,inod) + (/ m, n, o /) + if (.not.lat_points(vlat(1), vlat(2),vlat(3))) cycle latloop + call get_interp_pos(vlat(1), vlat(2), vlat(3), ie, rfill(:,:,inod)) + end do + + !Check to make sure all lattice points exist for the current element + if(any(.not.lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1))) cycle latloop + if (.not.in_group_ele(esize, lat_ele(ie), rfill)) then + nump_ele=nump_ele - esize**3 + lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1) = .false. + call add_element(0,type_ele(ie), esize, lat_ele(ie), sbox_ele(ie), rfill) + new_ele_num = new_ele_num + 1 + added_points = added_points + esize**3 + end if + + end do latloop + end do + end do + esize=esize-2 + end if + end do + !Now add the leftover lattice points as atoms + do o = 1, size_ele(ie) + do n = 1, size_ele(ie) + do m = 1, size_ele(ie) + if(lat_points(m,n,o)) then + call get_interp_pos(m,n,o, ie, ratom(:,:)) + do ibasis = 1, basisnum(lat_ele(ie)) + call apply_periodic(ratom(:,ibasis)) + call add_atom(0, basis_type(ibasis,lat_ele(ie)), sbox_ele(ie), ratom(:,ibasis)) + added_points=added_points + 1 + end do + end if + end do + end do + end do + + if (added_points /= (size_ele(ie)**3)) then + + print *, "Element ", ie, " is refined incorrectly in refinefill" + end if + end do + !Once all atoms are added we delete all of the elements + call delete_elements(group_ele_num, element_index) + print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms and ", new_ele_num, & + " elements." + end if + + end subroutine refinefill_group subroutine remesh_group !This command is used to remesh the group to a desired element size @@ -705,7 +931,8 @@ module opt_group !Add the element, for the sbox we just set it to the same sbox that we get the orientation !from. In this case it is from the sbox of the first atom in the group. new_ele = new_ele+1 - call add_element(remesh_ele_type, working_esize, ilat, sbox_atom(atom_index(1)),r_new_node) + call add_element(0,remesh_ele_type, working_esize, ilat, & + sbox_atom(atom_index(1)),r_new_node) end if end if @@ -726,7 +953,7 @@ module opt_group lat_points(ix,iy,iz) = .false. r = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat) new_atom=new_atom+1 - call add_atom(basis_type(1,ilat), is, r) + call add_atom(0,basis_type(1,ilat), is, r) end if end do end do @@ -755,11 +982,40 @@ module opt_group call delete_elements(group_ele_num, element_index) end subroutine delete_group + + subroutine change_group_type + !This subroutine changes all atoms and nodes at atoms within a group to a specific type + integer :: i, j, ltype,ibasis, inod, basis_type(10) + + print *, "Changing ", group_atom_num, " atoms and ", group_ele_num, " elements to atom type ", group_type + + !Change all atom group types + do i = 1, group_atom_num + j = atom_index(i) + type_atom(j) = group_type + end do + + !Map to a new lattice type for all element + do i =1, group_ele_num + j = element_index(i) + ltype = lat_ele(j) + do ibasis = 1, basisnum(ltype) + basis_type(ibasis) = group_type + end do + call lattice_map(basisnum(ltype), basis_type, ng_node(ltype), lapa(ltype), lat_ele(j)) + end do + + end subroutine change_group_type + + subroutine split_group_elements + ! + end subroutine split_group_elements + function in_group(r) !This subroutine determines if a point is within the group boundaries real(kind=dp), intent(in) :: r(3) - real(kind=dp) :: r_notch + real(kind=dp) :: rnorm integer :: dim3, i logical :: in_group @@ -779,9 +1035,95 @@ module opt_group if (r(dim1) < centroid(dim1)) in_group = .false. end if - r_notch = sqrt((r(dim1) - centroid(dim1))**2 + (r(dim3)-centroid(dim3))**2) - in_group = in_group.or.(r_notch < tip_radius) + rnorm = sqrt((r(dim1) - centroid(dim1))**2 + (r(dim3)-centroid(dim3))**2) + in_group = in_group.or.(rnorm < radius) - end select + case('sphere') + rnorm = norm2(r(:) - centroid(:)) + if (rnorm <= radius) then + in_group = .true. + else + in_group = .false. + end if + case('shell') + rnorm = norm2(r(:) - centroid(:)) + if ((rnorm >= radius).and.(rnorm<=(radius + shell_thickness))) then + in_group = .true. + else + in_group = .false. + end if + case('all') + in_group = .true. + end select end function in_group + + function in_group_ele(esize, elat, rn) + !This figures out if the elements are in the group boundaries + real(kind=dp), intent(in) :: rn(3,max_basisnum, max_ng_node) + integer, intent(in) :: esize, elat + logical :: in_group_ele + + integer :: i, inod, ibasis, node_in_out(max_ng_node) + real(kind=dp) :: r_center(3) + + in_group_ele=.false. + + if(trim(adjustl(shape)) == 'shell') then + node_in_out(:) = -1 + !First calculate whether each element node is within the shell region, inside the shell sphere, or outside the + !shell region + nodeloop:do inod = 1, ng_node(elat) + r_center(:)=0.0_dp + do ibasis = 1, basisnum(elat) + r_center(:)= r_center(:) + rn(:,ibasis,inod)/basisnum(elat) + end do + + if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(size_ele(i)/=notsize)) then + node_in_out(inod) = 2 + exit nodeloop + end if + + shape ='sphere' + if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(esize/=notsize)) then + node_in_out(inod) = 1 + else + node_in_out(inod) = 0 + end if + shape='shell' + end do nodeloop + + !If any nodes are within the shell region, or if the shell region interescts an element then add it to the group + if(any(node_in_out == 2).or.(any(node_in_out==1).and.(any(node_in_out==0)))) then + in_group_ele=.true. + return + end if + + else if(.not.(group_nodes)) then + r_center(:) = 0.0_dp + do inod = 1, ng_node(elat) + do ibasis = 1, basisnum(elat) + r_center = r_center + rn(:,ibasis,inod)/(basisnum(elat)*ng_node(elat)) + end do + end do + + if ((in_group(r_center).neqv.flip).and.(esize/= notsize)) then + in_group_ele=.true. + return + end if + + else if(group_nodes) then + r_center(:) = 0.0_dp + do inod = 1, ng_node(elat) + do ibasis = 1, basisnum(elat) + if ((in_group(rn(:,ibasis,inod)).neqv.flip).and.(esize/=notsize)) then + in_group_ele=.true. + return + end if + end do + end do + end if + + + end function in_group_ele + end module opt_group diff --git a/src/opt_orient.f90 b/src/opt_orient.f90 index 4ec4b4c..fd42cb1 100644 --- a/src/opt_orient.f90 +++ b/src/opt_orient.f90 @@ -11,18 +11,19 @@ module opt_orient real(kind=dp), save :: new_orient(3,3) real(kind=dp), dimension(6) :: orig_box_bd real(kind=dp), allocatable :: orig_sub_box_ori(:,:,:) + character(len=3) :: orig_box_bc public contains - subroutine orient(arg_pos) + subroutine orient_opt(arg_pos) integer, intent(inout) :: arg_pos - integer :: i, ibasis, inod - logical :: isortho, isrighthanded + integer :: i, j, k, ibasis, inod + logical :: isortho, isrighthanded, matching real(kind=dp) :: inv_sub_box_ori(3,3,sub_box_num) - character(len=3) :: old_box_bc + !First parse the orient command call parse_orient(arg_pos) @@ -63,12 +64,25 @@ module opt_orient !Save original box boundaries orig_box_bd = box_bd - !Now find new box boundaries, have to temporarily define the box as shrink wrapped for def new box to work - old_box_bc = box_Bc - box_bc = 'sss' - call def_new_box - box_bc = old_box_bc - end subroutine orient + !Now find new box boundaries, if any orientations are the same we leave them as they are. If they are different then we have + !to shrink wrap them + + orig_box_bc = box_bc + do i = 1,3 + matching=.true. + sbox_loop:do j = 1, sub_box_num + do k = 1, 3 + if(.not.is_equal(orig_sub_box_ori(i,k,j), new_orient(i,k))) then + matching = .false. + exit sbox_loop + end if + end do + end do sbox_loop + if(.not.matching) box_bc(i:i)='s' + end do + + call def_new_box + end subroutine orient_opt subroutine parse_orient(arg_pos) !This parses the orient option @@ -129,8 +143,9 @@ module opt_orient end do end do - !Restore original box boundaries + !Restore original box boundaries and box BC box_bd = orig_box_bd + box_bc = orig_box_bc end subroutine unorient subroutine sbox_ori(arg_pos) diff --git a/src/opt_redef_box.f90 b/src/opt_redef_box.f90 new file mode 100644 index 0000000..f6330e0 --- /dev/null +++ b/src/opt_redef_box.f90 @@ -0,0 +1,131 @@ +module opt_redef_box + + use box + use elements + use subroutines + implicit none + + character(len=3) :: redef_dim, new_bc + real(kind=dp) :: new_bd(6) + public + contains + + subroutine redef_box(arg_pos) + !This is the main calling function for opt_redef_box + integer, intent(inout) :: arg_pos + integer :: i, inod, ibasis, iatom, delete_list(atom_num), delete_num, type_interp(max_basisnum*max_esize**3) + real(kind=dp):: r_interp(3, max_basisnum*max_esize**3) + logical :: node_out(8) + + !First parse the argument + call parse_redef_box(arg_pos) + + print *, '------------------------------------------------------------' + print *, 'Option redef_box' + print *, '------------------------------------------------------------' + + !Now first filter atoms that don't fit in the new box bounds and delete them + delete_num = 0 + do i = 1, atom_num + if(.not.in_block_bd(r_atom(:,i),new_bd)) then + delete_num = delete_num + 1 + delete_list(delete_num) = i + end if + end do + call delete_atoms(delete_num, delete_list(1:delete_num)) + + !Now loop over all elements + delete_num = 0 + delete_list(:) = 0 + do i = 1, ele_num + !Determine if all nodes are within the new boundaries + node_out(:) = .false. + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + if(.not.in_block_bd(r_node(:,ibasis,inod,i), new_bd)) then + node_out(inod) = .true. + end if + end do + end do + + !If all nodes are out just add the element to the delete list + if(all(node_out)) then + delete_num = delete_num +1 + delete_list(delete_num) = i + + !If any nodes are out we add the element to the delete list, but then loop over the interpoalted atoms to figure out + !which ones fit inside the boundary to keep the box rectangular + else if (any(node_out)) then + delete_num = delete_num +1 + delete_list(delete_num) = i + + call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp) + + !loop over all interpolated atoms and add them to the system + do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 + if(in_block_bd(r_interp(:,iatom), new_bd)) then + call add_atom(0,type_interp(iatom), sbox_ele(i), r_interp(:,iatom)) + end if + end do + end if + end do + + call delete_elements(delete_num, delete_list(1:delete_num)) + + print *, "Old box_bd: ", box_bd, " is redefined to new box boundaries: ", new_bd + box_bd=new_bd + box_bc = new_bc + + end subroutine redef_box + + subroutine parse_redef_box(arg_pos) + !Parse the command + integer, intent(inout) :: arg_pos + + integer :: i, j, arglen + character(len=100) textholder + + + !First read in the dimensions that we are redefining + redef_dim = '' + arg_pos=arg_pos+1 + call get_command_argument(arg_pos, redef_dim, arglen) + select case(trim(adjustl(redef_dim))) + case('x','y','z','xy','yx','xz','zx','yz','zy','xyz','yxz','xzy','zyx','zxy','yzx') + continue + case default + print *, "Incorrect redef_dim ", redef_dim, "please select any permuation of x, y, z, xy, yz, xz, xyz" + stop 3 + end select + + !Now read in the new dimensions + new_bd = box_bd + new_bc = box_bc + do i = 1, 3 + select case(trim(adjustl(redef_dim(i:i)))) + case('x') + j = 1 + case('y') + j = 2 + case('z') + j = 3 + case default + exit + end select + + arg_pos=arg_pos +1 + call get_command_argument(arg_pos, textholder, arglen) + if(arglen == 0) stop "Missing a box dimension in opt_redef_box" + call parse_pos(j, textholder,new_bd(2*j-1)) + + arg_pos=arg_pos +1 + call get_command_argument(arg_pos, textholder, arglen) + if(arglen == 0) stop "Missing a box dimension in opt_redef_box" + call parse_pos(j, textholder,new_bd(2*j)) + new_bc(j:j) = 's' + end do + + arg_pos = arg_pos + 1 + end subroutine parse_redef_box + +end module opt_redef_box diff --git a/src/opt_slip_plane.f90 b/src/opt_slip_plane.f90 new file mode 100644 index 0000000..80da930 --- /dev/null +++ b/src/opt_slip_plane.f90 @@ -0,0 +1,176 @@ +module opt_slip_plane + use parameters + use elements + use functions + use subroutines + + implicit none + + integer :: sdim + real(kind=dp) :: spos + logical :: efill + + public + contains + + subroutine run_slip_plane(arg_pos) + !Main calling function for the slip_plane option + integer, intent(inout) :: arg_pos + + integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3), ibasis + + integer, allocatable :: slip_eles(:), temp_int(:) + real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node), ratom(3,max_basisnum), & + maxp, minp + + integer :: type_interp(max_basisnum*max_esize**3) + logical :: lat_points(max_esize,max_esize, max_esize) + + + print *, '---------------------Option Slip_Plane----------------------' + + !Initialize variables + efill = .false. + slip_enum = 0 + old_atom_num = atom_num + + !!Parse the argument + call parse(arg_pos) + + + !If we are running the efill code then we have to initiate some variables + if(efill) then + new_ele_num = 0 + end if + allocate(slip_eles(1024)) + !Now loop over all elements, find which ones intersect + do ie = 1, ele_num + if( (spos < maxval(r_node(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)),ie))).and. & + (spos > minval(r_node(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)),ie)))) then + slip_enum = slip_enum + 1 + if (slip_enum > size(slip_eles)) then + allocate(temp_int(size(slip_eles)+1024)) + temp_int(1:size(slip_eles)) = slip_eles + temp_int(size(slip_eles)+1:) = 0 + call move_alloc(temp_int, slip_eles) + end if + slip_eles(slip_enum) = ie + + !If we aren't efilling then just refine the element + if(.not.efill) then + call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) + do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3 + call apply_periodic(r_interp(:,ia)) + call add_atom(0, type_interp(ia), sbox_ele(ie), r_interp(:,ia)) + end do + !If we are efilling then the code is slightly more complex + else + !First populate the lat points array + lat_points(1:size_ele(ie),1:size_ele(ie), 1:size_ele(ie)) = .true. + + !Now start trying to remesh the region, leaving the slip plane as a discontinuity + esize = size_ele(ie) - 2 + nump_ele = size_ele(ie)**3 + do while(esize > min_efillsize) + if(nump_ele < esize**3) then + esize = esize - 2 + else + ele = cubic_cell*(esize -1) + do o = 1, size_ele(ie) - esize + do n = 1, size_ele(ie) - esize + latloop:do m = 1, size_ele(ie) - esize + do inod = 1, ng_node(lat_ele(ie)) + vlat = ele(:,inod) + (/ m, n, o /) + if (.not.lat_points(vlat(1), vlat(2),vlat(3))) cycle latloop + call get_interp_pos(vlat(1), vlat(2), vlat(3), ie, rfill(:,:,inod)) + end do + + !Check to make sure all lattice points exist for the current element + if(any(.not.lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1))) cycle latloop + !Check to see if the plane intersects this element if not then add it + maxp = maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)))) + minp = minval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)))) + if(.not.(spos < maxp).and.(spos > minp))then + nump_ele = nump_ele - esize**3 + lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1) = .false. + call add_element(0, type_ele(ie), esize, lat_ele(ie), sbox_ele(ie), rfill) + new_ele_num = new_ele_num + 1 + end if + end do latloop + end do + end do + end if + esize= esize-2 + end do + ! Now add the leftover lattice points as atoms + do o = 1, size_ele(ie) + do n = 1, size_ele(ie) + do m = 1, size_ele(ie) + if(lat_points(m,n,o)) then + call get_interp_pos(m,n,o, ie, ratom(:,:)) + do ibasis = 1, basisnum(lat_ele(ie)) + call apply_periodic(r_atom(:,ibasis)) + call add_atom(0, basis_type(ibasis,lat_ele(ie)), sbox_ele(ie), ratom(:,ibasis)) + end do + end if + end do + end do + end do + end if + end if + end do + + !Once we finish adding elements delete the old ones + call delete_elements(slip_enum, slip_eles(1:slip_enum)) + + !Output data + if(.not.efill) then + print *, "We refine ", slip_enum, " elements into ", atom_num - old_atom_num , " atoms" + else + print *, "We refine ", slip_enum, " elements into ", atom_num - old_atom_num , " atoms and ", new_ele_num, " elements" + end if + + end subroutine run_slip_plane + + subroutine parse(arg_pos) + !This subroutine parses the input arguments to the mode + integer, intent(inout) :: arg_pos + + integer :: arglen + character(len = 100) :: textholder + + !First read the dimension + arg_pos = arg_pos +1 + call get_command_argument(arg_pos,textholder, arglen) + if(arglen == 0) stop "Incorrect slip_plane command. Please check documentation for correct format" + + !Check to make sure that the dimension is correct + select case(trim(adjustl(textholder))) + case('x','X') + sdim = 1 + case('y','Y') + sdim = 2 + case('z','Z') + sdim = 3 + case default + print *, "Error: dimension ", trim(adjustl(textholder)), " is not accepted. Please select from x, y, or z" + end select + + !Now parse the position of the slip plane + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if(arglen == 0) stop "Incorrect slip_plane command. Please check documentation for correct format" + call parse_pos(sdim, textholder, spos) + + !Now check to see if efill was passed + arg_pos = arg_pos + 1 + if(.not.(arg_pos > command_argument_count())) then + call get_command_argument(arg_pos, textholder, arglen) + if(arglen == 0) stop "Incorrect slip_plane command. Please check documentation for correct format" + if(trim(adjustl(textholder)) == "efill") then + arg_pos = arg_pos +1 + efill = .true. + end if + end if + end subroutine parse +end module opt_slip_plane diff --git a/src/parameters.f90 b/src/parameters.f90 index f261552..c589b4b 100644 --- a/src/parameters.f90 +++ b/src/parameters.f90 @@ -3,7 +3,8 @@ module parameters implicit none !Default precision - integer, parameter :: dp= selected_real_kind(15,307) + integer, parameter :: dp= selected_real_kind(15,307), & + min_efillsize = 11 !Parameters for floating point tolerance real(kind=dp), parameter :: lim_zero = epsilon(1.0_dp), & lim_large = huge(1.0_dp), & @@ -12,5 +13,7 @@ module parameters !Numeric constants real(kind=dp), parameter :: pi = 3.14159265358979323846_dp - + + !Mode type that is being called + character(len=100) :: mode end module parameters diff --git a/src/str.f90 b/src/str.f90 new file mode 100644 index 0000000..2d0ed73 --- /dev/null +++ b/src/str.f90 @@ -0,0 +1,33 @@ +module str + + !this module has some string manipulation commands + public + contains + + pure function tok_count(text) + !counts number of tokens in a string + character(len = *), intent(in) :: text + integer :: tok_count + integer :: i, j + logical :: in_tok + + j = len(trim(adjustl(text))) + in_tok = .false. + tok_count = 0 + do i = 1, j + !This checks if it is a white space character which is the delimiter + if(trim(adjustl(text(i:i))) == ' ') then + !If previously we were in token and the current character is the delimiter + !Then we are no longer in the token + if(in_tok) in_tok = .false. + + !If the character isn't a white space character and we previously weren't in the token then set in_tok + !to true and increment token count + else if(.not.in_tok) then + in_tok = .true. + tok_count = tok_count + 1 + end if + end do + return + end function tok_count +end module str diff --git a/src/subroutines.f90 b/src/subroutines.f90 index 1015968..b272c16 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -198,66 +198,6 @@ module subroutines end do end subroutine - subroutine build_cell_list(numinlist, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell) - !This subroutine builds a cell list based on rc_off - - !----------------------------------------Input/output variables------------------------------------------- - - integer, intent(in) :: numinlist !The number of points within r_list - - real(kind=dp), dimension(3,numinlist), intent(in) :: r_list !List of points to be used for the construction of - !the cell list. - - real(kind=dp), intent(in) :: rc_off ! Cutoff radius which dictates the size of the cells - - integer, dimension(3), intent(inout) :: cell_num !Number of cells in each dimension. - - integer, allocatable, intent(inout) :: num_in_cell(:,:,:) !Number of points within each cell - - integer, allocatable, intent(inout) :: cell_list(:,:,:,:) !Index of points from r_list within each cell. - - integer, dimension(3,numinlist), intent(out) :: which_cell !The cell index for each point in r_list - - !----------------------------------------Begin Subroutine ------------------------------------------- - - integer :: i, j, cell_lim, c(3) - real(kind=dp) :: box_len(3) - integer, allocatable :: resize_cell_list(:,:,:,:) - - !First calculate the number of cells that we need in each dimension - do i = 1,3 - box_len(i) = box_bd(2*i) - box_bd(2*i-1) - cell_num(i) = int(box_len(i)/(rc_off/2))+1 - end do - - !Initialize/allocate variables - cell_lim = 10 - allocate(num_in_cell(cell_num(1),cell_num(2),cell_num(3)), cell_list(cell_lim, cell_num(1), cell_num(2), cell_num(3))) - - !Now place points within cell - do i = 1, numinlist - !c is the position of the cell that the point belongs to - do j = 1, 3 - c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off/2)) + 1 - end do - - !Place the index in the correct position, growing if necessary - num_in_cell(c(1),c(2),c(3)) = num_in_cell(c(1),c(2),c(3)) + 1 - if (num_in_cell(c(1),c(2),c(3)) > cell_lim) then - allocate(resize_cell_list(cell_lim+10,cell_num(1),cell_num(2),cell_num(3))) - resize_cell_list(1:cell_lim, :, :, :) = cell_list - resize_cell_list(cell_lim+1:, :, :, :) = 0 - call move_alloc(resize_cell_list, cell_list) - end if - - cell_list(num_in_cell(c(1),c(2),c(3)),c(1),c(2),c(3)) = i - which_cell(:,i) = c - end do - - return - end subroutine build_cell_list - - subroutine check_right_ortho(ori, isortho, isrighthanded) !This subroutine checks whether provided orientations in the form: ! | x1 x2 x3 | @@ -301,4 +241,5 @@ module subroutines return end subroutine check_right_ortho + end module subroutines