Latest updates

.gitlab-ci.yml

@@ -0,0 +1,27 @@
+image: python:3.8-buster
+
+before_script:
+  - pip install -r requirements.txt
+
+stages: 
+  - tests
+  - deploy
+
+pagestests:   
+  stage: tests   
+  script:   
+    - mkdocs build --strict --verbose --site-dir test   
+  artifacts:     
+    paths:     
+      - test   
+  except:   
+    - development
+pages:
+  stage: deploy
+  script:
+  - mkdocs build --strict --verbose
+  artifacts:
+    paths:
+    - public
+  only:
+  - development

README.md

@@ -3,378 +3,7 @@
 
 This is a tool for building models in CAC. Commands and usage options are below. This code is intended to follow the atomsk code fairly closely.
 
-## Modes
+Documentation is available at: [aselimov.gitlab.io/cacmb](https://aselimov.gitlab.io/cacmb)
 
-The modes follow similarly to the modes you find when using atomsk. The modes will be listed below alongside their syntax and other usage instructions. As a note, if a mode is being used then it has to come first.
+A paper describing the use of this tool is currently in preparation.
 
-### Mode Create
-
-``` 
-cacmb --create name element_type lattice_parameter esize
-```
-
-Mode create has the following parameters:
-
-`name` - User defined name that either defines the atom type or the lattice type if using the basis option
-
-`element_type` - Specifies which element type to use, this dictates the crystal being build. Current   acceptable options for element_type are:
-
-* FCC - Uses the Rhombohedral primitive fcc unit cell as the finite element.
-
-`lattice_parameter` - The lattice parameter for the crystal structure.
-
-`esize` - Number of atoms per edge of the finite element. A value of 2 signifies full atomistic resolution and is the lowest number acceptable.
-
-**Example**
-
-```
-cacmb --create Cu fcc 3.615 11
-```
-
-Creates a copper element with a lattice parameter of 3.615 with 11 atoms per side
-
-#### Optional keywords
-
-**Orient**
-
-``` orient [hkl] [hkl] [hkl] 
-orient [hkl] [hkl] [hkl]
-```
-
-Default orientation is `[100] [010] [001]`. If this keyword is present then the user must provide the orientation matrix in form `[hkl] [hkl] [hkl]`.
-
-*Example:* `orient [-112] [110] [-11-1]` 
-
-**Duplicate**
-
-```
-duplicate numx numy numz
-```
-
-Default duplicate is `1 1 1`. This is used to replicate the element along each dimensions. This cannot be used if the keyword dimensions is included. By default jagged edges along boundaries are filled if duplicate is greater than `1 1 1`.
-
-*Example:* `duplicate 10 10 10`
-
-**Dimensions**
-
-```
-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.
-
-Example: `dimensions 100 100 100` 
-
-**Origin**
-
-```
-origin x y z
-```
-
-Default origin is `0 0 0`. This command just sets the origin for where the simulation cell starts building.
-
-*Example:* `origin 10 0 1`
-
-**Basis**
-
-```
-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
-
-```
-cacmb --convert infile outfile
-```
-
-This mode converts a file `infile` to a file of  `outfile`. The extensions determine the conversion process.
-
-### Mode Merge
-
-```
-cacmb --merge dim N infiles outfile
-```
-
-This mode merges multiple data files and creates one big simulation cell. The parameters are:
-
-`N` - The number of files which are being read
-
-`dim` - the dimension they are to be stacked along, can be either `x`, `y`, or `z`. If the argument `none` is passed then the cells are just overlaid. Future options will include a delete overlap command.
-
-**Additional options:**
-
-**Shift**
-
-```
-shift x y z
-```
-
-If the shift command is passed to mode merge then each file after the first file in the merge command is displaced by the vector `[x, y, z]`. This is additive so if you are merging three files and this command is passed then the second file is shifted by `[x,y,z]` and the third file is shifted by `2*[x,y,z]`. 
-
-Example: `cacmb --merge z 2 Cu.mb Cu2.mb Cu3.mb Cumerged.mb shift 2 0 0` will shift the atomic and element positions in the `Cu2.mb` file by `[2,0,0]` and the positions in `Cu3.mb` by `[4,0,0]`.
-
-**Wrap**
-
-```
-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.
-
-### Option dislgen
-
-```
--dislgen [ijk] [hkl] x y z burgers char_angle poisson 
-```
-
-
-
-This options adds an arbitrarily oriented dislocation into your model based on user inputs using the volterra displacement fields. The options are below
-
-`[ijk]` - The vector for the line direction
-
-`[hkl]` - The vector for the slip plane 
-
-`burgers` - The magnitude of the burgers vector for the dislocation to be inserted
-
-`x y z` - The position of the dislocation centroid
-
-`char_angle` - Character angle of the dislocation (0 is screw and 90 is edge)
-
-`poisson` - Poisson's ratio used for the displacement field.
-
-### Option disloop
-
-````
--disloop loop_normal radius x y z bx by bz poisson
-````
-
-This option imposes the displacement field for a dislocation in order to create a loop. This loop is unstable and will close if stress isn't applied.
-
-`loop_normal` - The box dimension which defines the normal to the loop plane. As of now this dimension must be a closed back direction, meaning that for fcc a box dimension has to be of the (111) family of planes. Either `x`, `y`, or `z`.
-
-`n` - The number of atoms to delete on the loop plane
-
-`x y z` - The centroid of the loop. 
-
-`bx by bz` - The burgers vector for the dislocation
-
-`poisson` - Poisson ratio for continuum solution
-
-### Option vacancy_disloop
-
-```
--vacancydisloop loop_normal radius x y z 
-```
-
-This option creates a circular planar vacancy cluster of radius `radius` normal to the `loop_normal` centered on position `x y z`. Upon relaxing or energy minimization this cluster should become a prismatic dislocation loop.
-
-### Option Group
-
-`-group select_type group_shape shape_arguments additional keywords`
-
-This option selects a group of either elements or atoms and applies some transformation to them.
-
-    `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 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`. 
-
-`additional keywords`- Represents the various transformations which can be performed on a group. These additional keywords are given below.
-
-*Wedge:*
-
-`-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.
-`bx by bz` - Centroid of the center of the base
-`bw` - Base width
-
-*Notch:*
-
-`-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.
-`bx by bz` - Centroid of the center of the base
-`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:**
-
-```
-displace x y z 
-```
-
-This is similar to the mode merge `-shift` argument. This simply shift atoms and elements within the group by a vector `[x,y,z]`.
-
-**Wrap** 
-
-```
-wrap	
-```
-
-This command wraps atoms back into the simulation cell as though periodic boundary conditions are being used. 
-
-**Remesh**
-
-```
-remesh esize 
-```
-
-This command remeshes the atoms/elements within the group to the new element size `esize`.
-
-**Max**
-
-```
-max
-```
-
-This command attempts to reduce the degrees of freedom in the model by replacing them with graded elements. This code works by starting at elements with size `esize` and then checks all degrees of freedom to see which ones can be replaced by inserting the element. It then iterates over elements of `esize-2` to `esize=2` which is full atomic resolution.
-
-**Delete**
-
-```
-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 
-
-```
--ow
-```
-
-If this option is passed then all files are automatically overwritten without asking the user.
-
-### Option boundary
-
-``` 
--boundary box_bc
-```
-
-This allows the user to specify the boundary conditions for the model being outputted. The format is a 3 character string with `p` indicating periodic and `s` indicating shrink-wrapped.
-
-**Example:** `-boundary psp`
-
-### Option Delete
-
-```
--delete keywords
-```
-
-Delete requires the usage of additional keywords to specify which delete action will be taken. These additional keywords are below:
-
-**overlap**
-
-```
--delete overlap rc_off
-```
-
-This command will delete all overlapping atoms within a specific cutoff radius `rc_off`. This currently does not affect elements. 
-
-### Option sbox_ori
-
-```
--sbox_ori sbox [hkl] [hkl] [hkl]
-```
-
-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.
-
-`val` - Where `val` is a number, then that value in Angstroms is used as the position. As an example, `11.1` would be read in as a position of 11.1 $\AA$.
-
-`-inf` - This specifies the lower box boundary in the specific dimension. The vector `-inf -inf -inf`  specifies the bottom corner of the simulation cell which also acts as the simulation cell origin. The vector `-inf 10 3` instead puts only the x position at the simulation cell origin.
-
-`inf` - Similar to `-inf` but references the upper boundary of the box in that dimension
-
-`inf-val` - Using a minus sign reduces the position from the **upper boundary** by `val`. `inf-10` would be at a distance of $10 \AA$ from the upper boundary in that dimension.
-
-`inf+val` - This increases the position from the **lower boundary**. `inf+10` would be a position $10\AA$ from the lower boundary within the cell.
-
-`inf*val` - This gives you a fractional position in the simulation cell. As an example `inf*0.5` gives you the center point of the simulation cell. 
-
-`rand` - Returns a random position that lies within the simulation cell.
-
-`rand[val1:val2]` - returns a random position that lies within the range
-
-`rande[facenum]` - Returns a random position in an interelement boundary which is offset of the element face `facenum`. Face numbers are based on the which vertices comprise the face. Vertex numbers are shown in the figure below for the primitive fcc unit cell which is what the fcc rhombohedral element is based from. The face numbers are:
-
-Face 1: [1,2,3,4]
-
-Face 2: [1,2,6,5]
-
-Face 3: [2,3,7,6]
-
-Face 4: [3,4,8,7]
-
-Face 5: [1,4,8,5]
-
-Face 6: [5,6,7,8]
-
-Image for vertex numbers is: 
-
-![](/home/alexselimov/Documents/CACmb/Numbered_element.png)

docs/Examples/semi-coherent.md

@@ -0,0 +1,24 @@
+# CuNi Bicrystal Interface
+This example builds a bicrystal of CuNi with a semi-coherent interface. The interface region is rendered at full atomistic resolution and the bulk of the layers are coarse-grained.
+
+```sh
+#!/bin/bash
+
+#Build Copper CG region and atom pad
+cacmb --create Cu fcc 3.615 25 orient [-21-1] [011] [11-1] duplicate 7.44 6 4 Cu_cg.mb -ow 
+cacmb --create Cu fcc 3.615 2 orient [-21-1] [011] [11-1] duplicate 93 75 3 Cu_atom.mb -ow
+#Build Nickel CG region and atom pad
+cacmb --create Ni fcc 3.52 25 orient [-21-1] [011] [11-1] duplicate 7.6409 6.16 4 Ni_cg.mb -ow 
+cacmb --create Ni fcc 3.52 2 orient [-21-1] [011] [11-1] duplicate 95.511 77 3 Ni_atom.mb -ow
+#Merge all into one block keeping full atomistic resolution at the interface
+cacmb --merge z 4  Ni_cg.mb Ni_atom.mb Cu_atom.mb Cu_cg.mb -boundary pps cac_in.restart cac_in.xyz  -ow 
+
+#Delete leftover files
+rm Ni_cg.mb Ni_atom.mb Cu_atom.mb Cu_cg.mb 
+
+```
+
+Below are the .xyz file render which shows the atomistic and coarse-grained regions and the common neighbor analysi of the interface.
+
+![xyz bicrystal render](../img/CuNi_xyz.png) 
+![interface structure for bicrystal model](../img/CuNi_int.png) 

docs/Getting_Started.md

@@ -0,0 +1,48 @@
+# Getting Started
+
+The CACmb tool separates commands into two general groups, modes and options. 
+Modes are the main commands which conduct primary operations, these are listed below:
+
+* [Mode Convert](Modes/convert.md)
+* [Mode Create](Modes/create.md) 
+* [Mode Merge](Modes/merge.md) 
+
+To start a cacmb run, it is required to run one of the available modes and only one mode can be run at a time.
+Additionally one can use options in order to adjust the model during any of the mode runs and unlimited number of options may be used.
+While single line commands may be used satisfactorily to build a model, in order to build more complex models it is required to run multiple commands.
+
+Currently CACmb only supports rhombehedral FCC finite elements with faces aligned to {111} slip planes.
+This finite element is the primitive unit cell for the FCC crystal and is shown below:
+
+![Rhombohedral Element](img/rhomb.png) 
+
+## Simple Example 1
+For example, to build a model with an atomistic region bound by a coarse-grained region the following commands should be run:
+
+```
+cacmb --create Cu fcc 3.615 15 duplicate 4 2 1 orient [112] [1-10] [11-1] Cu_cg.mb -ow
+cacmb --create Cu fcc 3.615 2 duplicate 30 15 2 orient [112] [1-10] [11-1] Cu_at.mb -ow
+cacmb --merge z 3 Cu_cg.mb Cu_at.mb Cu_cg.mb model.xyz -ow
+```
+
+When visualizing the model.xyz file using ovito we see our desired model (where the red atoms represent the element nodes and blue atoms are atomistic regions):
+
+![Simple example 1 render](img/simple_example_1.png) 
+
+**Important Note: In order to fully utilize the cacmb capabilites use the .mb format when running multiple model build steps**
+
+A description of each build step is below:
+
+1. Build a coarse grained model with elements that have 15 atoms per edge. Orientation is set to `[112] [1-10] [11-1]`. We save it to `Cu_cg.mb` and pass `-ow` which is the flag to overwrite existing files with that name.
+2. Build an atomistic model with the same dimensionsin the x and y directions as the coarse grained model. *A general way to match dimensions is by setting the duplicate value of the atoms to `esize*duplicate/2`. So for the x dimension we do `15*4/2` which is 30 and for the y dimension `15*2/2` which is 15.
+3. Merge the models by stacking the coarse-grained and atomistic models. We start the stacking sequence from the bottom so the first region is the `Cu_cg.mb`, followed by the `Cu_at.mb` model in the middle, followed by another `Cu_cg.mb` region. We output it to model.xyz which renders the atoms and nodes of the model for visualization.
+
+This whole process can be easily run as a shell script as long as cacmb is on your path.
+
+```sh
+#!/bin/sh
+cacmb --create Cu fcc 3.615 15 duplicate 4 2 1 orient [112] [1-10] [11-1] Cu_cg.mb -ow
+cacmb --create Cu fcc 3.615 2 duplicate 30 15 2 orient [112] [1-10] [11-1] Cu_at.mb -ow
+cacmb --merge z 3 Cu_cg.mb Cu_at.mb Cu_cg.mb model.xyz -ow
+```
+

docs/Modes/convert.md

@@ -0,0 +1,13 @@
+#Convert
+
+```
+cacmb --convert infile outfile(s)
+```
+
+This mode converts a file `infile` to different files `outfile(s)`. 
+The extensions on both files determine which read and write subroutines are called.
+This mode can also be used to apply options to a file by setting `infile` and `outfile` to the same value such as:
+```
+cacmb --convert cac.mb cac.mb {-options}
+```
+

docs/Modes/create.md

@@ -0,0 +1,113 @@
+# Create
+ 
+``` 
+cacmb --create name element_type lattice_parameter esize
+```
+
+Mode create has the following parameters:
+
+`name` - User defined name that either defines the atom type or the lattice type if using the basis option. If the basis command is not called then name must be an element.
+
+`element_type` - Specifies which element type to use, this dictates the crystal being build. Current acceptable options for element_type are:
+
+* FCC - Uses the Rhombohedral primitive fcc unit cell as the finite element.
+
+`lattice_parameter` - The lattice parameter for the crystal structure.
+
+`esize` - Number of atoms per edge of the finite element. A value of 2 signifies full atomistic resolution and is the lowest number acceptable.
+
+**Example**
+
+```
+cacmb --create Cu fcc 3.615 11
+```
+
+Creates a copper fcc element with a lattice parameter of 3.615 with 11 atoms per side
+
+## Optional keywords
+
+Below are optional keywords that must directly follow after the esize command and before any options or filenames are passed
+
+### Orient
+
+``` 
+orient [hkl] [hkl] [hkl]
+```
+
+The default orientation that is built is `[100] [010] [001]`.
+If this keyword is present then the user must provide the orientation matrix in form `[hkl] [hkl] [hkl]` without spaces.
+Currently only accepts whole number values for `[hkl]`.
+
+*Example:* `orient [-112] [110] [-11-1]` 
+
+### Duplicate
+
+```
+duplicate numx numy numz
+```
+
+The default duplicate is `1 1 1`. 
+This is used to replicate elements along each dimension. 
+The unit of duplication is the lattice periodicity length times the size of elements being duplicated.
+`numx numy numz` dicate the number of duplications units in the x,y, and z dimensions.
+This keyword **cannot** be used with the `dim` command.
+
+*Example:* `duplicate 10 10 10`
+
+### Dim
+
+```
+dim dimx dimy dimz
+```
+ 
+There is no default `dim` as `duplicate` is the default option. 
+This command assigns a box with user-assigned dimensions and fills it with the desired element size.
+**Note:** using dim may not result in periodic cells unless the origin of the cell is shifted from 0,0,0.
+
+*Example:* `dimensions 100 100 100` 
+
+### Origin
+
+```
+origin x y z
+```
+
+This keyword is used to set the origin of the simulation cell.
+When using the `duplicate` command, the origin is set to the minimum position of the rotated element. 
+When using the `dim` command the default origin is `(0,0,0)`.
+When using the `dim` command for building it may become necessary to shift the origin of the simulation cell in order to produce a periodic simulation.
+This origin can be retrieved by first running `--create` with `duplicate 1 1 1`.
+
+*Example:* `origin 10 0 1`
+
+### Basis
+
+```
+basis basisnum bname bx by bz
+```
+
+This function allows you to define a custom basis to be used at every lattice point. 
+The parameters for this keyword are:
+
+-   `basisnum` - The number of atoms in the basis
+-   `bname` - The type of the basis atom (e.g. Cu or Si)
+-   `bx by bz` - The position of the basis atom.
+
+`bname bx by bz` must be repeated `basisnum` times for this command.
+
+*Example:* 
+`basis 2 Si 0 0 0 Si 1.35775 1.35775 1.35775` when used with the `fcc` element produces a diamond structure with Si.
+
+### efill
+```
+efill 
+```
+This command attempts to maximize the degree of coarse-graining by iterating through esizes smaller than the user specified when filling in the jagged boundaries to create a periodic block.
+This command will iterate through element sizes from the user specified `esize-2` to a minimum esize of 11. 
+Below is an example of a model without efill and a model with efill.
+
+*Model without efill* 
+![Model without efill command](../img/not_efilled_vtk.png) 
+
+*Model with efill* 
+![Model with efill command](../img/efilled_vtk.png) 

docs/Modes/merge.md

@@ -0,0 +1,28 @@
+# Merge
+
+```
+cacmb --merge dim N infiles 
+```
+
+This mode merges multiple data files and creates one big simulation cell. The parameters are:
+
+`dim` - the dimension they are to be stacked along, can be either `x`, `y`, or `z`. If the argument `none` is passed then the cells are just overlaid. 
+
+`N` - The number of files which are being read
+
+ `infiles` - The input files which are to be merged. There must be `N` of them.
+
+
+## Additional options:
+
+### Shift
+
+```
+shift x y z
+```
+
+If the shift command is passed to mode merge then each file after the first file in the merge command is displaced by the vector `[x, y, z]`. This is additive so if you are merging three files and this command is passed then the second file is shifted by `[x,y,z]` and the third file is shifted by `2*[x,y,z]`. 
+
+Example: `cacmb --merge z 2 Cu.mb Cu2.mb Cu3.mb Cumerged.mb shift 2 0 0` will shift the atomic and element positions in the `Cu2.mb` file by `[2,0,0]` and the positions in `Cu3.mb` by `[4,0,0]`.
+
+

docs/Options/boundary.md

@@ -0,0 +1,9 @@
+# Boundary
+
+``` 
+-boundary box_bc
+```
+
+This allows the user to specify the boundary conditions for the model being outputted. The format is a 3 character string with `p` indicating periodic and `s` indicating shrink-wrapped.
+
+**Example** `ppp` - is a fully periodic model, `pss` is periodic in the x dimension and shrink wrapped in y and z.

docs/Options/deform.md

@@ -0,0 +1,7 @@
+# Deform
+```
+-deform dim dl
+```
+This command is used to apply a uniaxial strain unto a simulation cell. 
+The argument `dim` takes the values of `x, y, or z` and the argument `dl` is the change in length to be applied (in Angstroms.
+

docs/Options/delete.md

@@ -0,0 +1,9 @@
+# Delete Overlap
+ ```
+-delete overlap rc_off
+ ```
+ 
+ This delete option removes overlapping atoms.
+ For every neighboring pair of atoms, the distance between them is checked.
+ If the distance is less than `rc_off` then the atom with the higher index in the data arrays is deleted.
+

docs/Options/disl.md

@@ -0,0 +1,52 @@
+# Dislocation
+
+There are various options for creating dislocations. 
+These are listed below.
+
+## dislgen
+ 
+```
+-dislgen [ijk] [hkl] x y z burgers char_angle poisson 
+```
+
+This options adds an arbitrarily oriented dislocation into your model based on user inputs using the volterra displacement fields. The options are below
+
+`[ijk]` - The vector for the line direction
+
+`[hkl]` - The vector for the slip plane 
+
+`burgers` - The magnitude of the burgers vector for the dislocation to be inserted
+
+`x y z` - The position of the dislocation centroid
+
+`char_angle` - Character angle of the dislocation (0 is screw and 90 is edge)
+
+`poisson` - Poisson's ratio used for the displacement field.
+
+## disloop
+
+````
+-disloop loop_normal radius x y z bx by bz poisson
+````
+
+This option imposes the displacement field for a dislocation in order to create a loop. This loop is unstable and will close if stress isn't applied.
+
+`loop_normal` - The box dimension which defines the normal to the loop plane. As of now this dimension must be a close packed direction, meaning that for fcc a box dimension has to be of the (111) family of planes. Either `x`, `y`, or `z`.
+
+`radius` - The radius of the loop in Angstroms.
+
+`x y z` - The centroid of the loop. 
+
+`bx by bz` - The burgers vector for the dislocation
+
+`poisson` - Poisson ratio for continuum solution
+
+## vacancy disloop
+
+```
+-vacancydisloop loop_normal radius x y z 
+```
+
+This option creates a circular planar vacancy cluster of radius `radius` normal to the `loop_normal` centered on position `x y z`. Upon relaxing or energy minimization this cluster should become a prismatic dislocation loop.
+
+

docs/Options/group.md

@@ -0,0 +1,129 @@
+# Group
+
+```
+-group select_type group_shape shape_arguments
+```
+
+This option selects a group of either elements or atoms and applies some transformation to them.
+
+    `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:
+
+## Group Shapes
+These are the allowed group shapes
+
+### Block
+```
+-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`. 
+
+`additional keywords`- Represents the various transformations which can be performed on a group. These additional keywords are given below.
+
+### Wedge
+```
+-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.
+`bx by bz` - Centroid of the center of the base
+`bw` - Base width
+
+### Notch
+```
+-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.
+`bx by bz` - Centroid of the center of the base
+`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`.
+
+## Group selection operators
+
+The following are a list of group operations which alter how the atoms/elements in the group are selected.
+
+
+### 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. 
+These random atoms/elements then form the new group.
+
+*Example:*
+```
+group atoms block -inf inf -inf inf inf*0.5 inf random 10
+```
+The above example selects 10 random atoms from the specified group.
+These atoms will all have z values which place them in the top half of the cell.
+
+### Node
+```
+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.
+
+### Notsize
+```
+notsize esize
+```
+Notsize filters out elements of size `esize` when determining the group. For example, if the user runs `notsize 25` and the group bounds contain elements with esizes of 13, 15, 25 (from the efill command for example) then only the elements with esize 13 and 15 are selected as part of the group.
+
+
+## Group modification operators
+
+These modifiers operate on selected atoms and elements to adjust their properties.
+
+### Displace 
+
+```
+displace x y z 
+```
+This operation shifts all grouped atoms/nodes by a vector `(x,y,z)`.
+
+### Delete
+
+```
+delete
+```
+This command deletes all selected atoms and elements within the group.
+
+### Type
+```
+type atom_type
+```
+This command changes all atoms and basis atoms at element nodes for the group elements and atoms to type `atom_type`. 
+`atom-type` should be two characters which describe the atomic element (such as `Cu` or `Ni`)
+
+### Refine
+```
+refine
+```
+This command refines all elements within the group to full atomistic resolution.
+
+

docs/Options/orient.md

@@ -0,0 +1,18 @@
+# Orient 
+
+``` 
+-orient [hkl] [hkl] [hkl]
+```
+
+This command transforms the simulation cell into the specified orientation. 
+This can be useful when attempting to view the model along different slip planes.
+In general the dimensions will not be periodic as the box will not have flat surfaces along the x, y, and z dimensions.
+*Example:* `-orient [-112] [110] [-11-1]` 
+
+This command also comes with a sub command `unorient`.
+## Unorient
+```
+-unorient
+```
+This removes the orientation applied by an orient command during the current cacmb run.
+

docs/Options/ow.md

@@ -0,0 +1,7 @@
+# Overwrite
+```
+-ow
+```
+
+If this option is passed then all files are automatically overwritten without asking the user.
+

docs/Options/redef_box.md

@@ -0,0 +1,8 @@
+# Redefine Box Boundaries 
+
+```
+-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

docs/Options/sbox_ori.md

@@ -0,0 +1,20 @@
+# Define Orientation
+
+```
+-sbox_ori sbox [hkl] [hkl] [hkl]
+```
+
+This option is primarily used when reading data from non .mb formats.
+This code sets the orientation of the elements and atoms that are part of the sub_box `sbox` to the user provided orientation.
+This doesn't change anything about the model, but is required in order to accurately use the `dislgen` or `orient` options when operating on non .mb files.
+
+*Example:* 
+
+```
+cacmb --create Cu fcc 3.615 2 orient [112] [1-10] [11-1] duplicate 10 10 10 cac.lmp
+cacmb --convert cac.lmp cac_out.lmp -sbox_ori 1 [112] [1-10] [11-1] -dislgen [112] [11-1] inf*0.5 inf*0.5 inf*0.5 2.556 90 0.3 
+```
+The above command attempts to insert a dislocation using the `dislgen` option. 
+In order to do the coordinate transformation, the code needs to know what the orientation of the model originally and so the sbox_ori command sets the sbox information. 
+To better understand how the sub_boxes work, please view sub_box in the Getting Started guide.
+

docs/Options/slip_plane.md

@@ -0,0 +1,14 @@
+# Slip Plane
+```
+-slip_plane dim pos 
+```
+
+This command forces a slip plane at position `pos` along dimension `dim` by finding which elements that plane intersects and refining them to full atomistics.
+
+## Optional Arguments 
+
+### efill
+```
+-slip_plane dim pos efill
+```
+The efill option attempts to remesh the elements that are refined while ensuring an element discontinuity on the prescribed slip plane

docs/Options/wrap.md

@@ -0,0 +1,8 @@
+# Wrap
+```
+-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.
+

docs/index.md

@@ -0,0 +1,31 @@
+# CACmb: The Concurent Atomistic-Continuum Model Builder
+
+CACmb is a tool used to build coarse-grained atomistic models for use with the PyCAC and LammpsCAC toolkits.
+This tool also allows for the creation of equivalent atomistic models to be used in lammps throught the .lmp data format.
+
+![cacmb image](img/demo.gif) 
+
+## Requirements
+CACmb requires only a fortran compiler and does not need any external libraries. 
+This code has been tested using both gfortran and ifort fortran compilers solely on linux.
+
+## Build instructions
+CACmb can be built by the following commands:
+```
+git clone https://gitlab.com/aselimov/cacmb
+cd cacmb/src
+make
+```
+You can also install cacmb to /usr/local/bin by running:
+```
+sudo make install
+```
+Otherwise make sure your cacmb executable is available on your system path.
+
+Please view the [Getting Started](Getting_Started.md) guide for an introduction and basic usage instructions.
+
+## License
+This tool is available under the MIT license which is printed below, a full copy of which is available at the cacmb [Gitlab repository.](https://gitlab.com/aselimov/cacmb/-/blob/development/LICENSE)
+
+
+

mkdocs.yml

@@ -0,0 +1,10 @@
+site_name: "CACmb: The Concurrent Atomistic Continuum Model Builder"
+site_url: https://aselimov.gitlab.io/cacmb
+site_dir: public
+theme: readthedocs
+
+use_directory_urls: false 
+
+extra_javascript:     
+    - https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.0/MathJax.js?config=TeX-AMS-MML_HTMLorMML 
+

requirements.txt

@@ -0,0 +1,6 @@
+# Documentation static site generator & deployment tool
+mkdocs>=1.1.2
+
+# Add your custom theme if not inside a theme_dir
+# (https://github.com/mkdocs/mkdocs/wiki/MkDocs-Themes)
+#  mkdocs-material>=5.4.0

src/Makefile

@@ -2,7 +2,8 @@
 .DEFAULT_GOAL := all
 
 FC=mpif90
-FFLAGS=-Wall -mcmodel=large -O0 -g -fbacktrace -fcheck=all -ffpe-trap=invalid,zero,overflow,underflow,denormal 
+#FFLAGS=-Wall -mcmodel=large -O0 -g -fbacktrace -fcheck=all -ffpe-trap=invalid,zero,overflow,underflow,denormal 
+FFLAGS=-O3 -g
 
 OBJDIR=obj
 SRCS := $(wildcard *.f90)
@@ -32,7 +33,7 @@ $(OBJDIR)/%.o: %.f90
 
 
 deps:
-	@fortdepend -o Makefile.dep -i mpi -b obj -w
+	 @makedepf90 -b obj *.f90 > Makefile.dep
 
 #----------------- CLEAN UP -----------------#
 

src/Makefile.dep

@@ -11,9 +11,11 @@ obj/main :  \
 	obj/mode_calc.o \
 	obj/mode_convert.o \
 	obj/mode_create.o \
+	obj/mode_da.o \
 	obj/mode_merge.o \
 	obj/mode_metric.o \
 	obj/neighbors.o \
+	obj/opt_bubble.o \
 	obj/opt_deform.o \
 	obj/opt_delete.o \
 	obj/opt_disl.o \
@@ -39,8 +41,10 @@ obj/caller.o :  \
 	obj/mode_calc.o \
 	obj/mode_convert.o \
 	obj/mode_create.o \
+	obj/mode_da.o \
 	obj/mode_merge.o \
 	obj/mode_metric.o \
+	obj/opt_bubble.o \
 	obj/opt_deform.o \
 	obj/opt_delete.o \
 	obj/opt_disl.o \
@@ -51,6 +55,7 @@ obj/caller.o :  \
 	obj/parameters.o
 
 obj/elements.o :  \
+	obj/atoms.o \
 	obj/box.o \
 	obj/functions.o \
 	obj/parameters.o \
@@ -94,10 +99,17 @@ obj/mode_create.o :  \
 	obj/parameters.o \
 	obj/subroutines.o
 
+obj/mode_da.o :  \
+	obj/elements.o \
+	obj/io.o \
+	obj/neighbors.o \
+	obj/parameters.o
+
 obj/mode_merge.o :  \
 	obj/atoms.o \
 	obj/elements.o \
 	obj/io.o \
+	obj/neighbors.o \
 	obj/parameters.o \
 	obj/subroutines.o
 
@@ -108,11 +120,19 @@ obj/mode_metric.o :  \
 	obj/parameters.o
 
 obj/neighbors.o :  \
+	obj/box.o \
 	obj/elements.o \
 	obj/functions.o \
 	obj/parameters.o \
 	obj/subroutines.o
 
+obj/opt_bubble.o :  \
+	obj/atoms.o \
+	obj/box.o \
+	obj/elements.o \
+	obj/opt_group.o \
+	obj/parameters.o
+
 obj/opt_deform.o :  \
 	obj/box.o \
 	obj/elements.o \
@@ -132,9 +152,11 @@ obj/opt_disl.o :  \
 	obj/subroutines.o
 
 obj/opt_group.o :  \
+	obj/atoms.o \
 	obj/box.o \
 	obj/elements.o \
 	obj/parameters.o \
+	obj/sorts.o \
 	obj/subroutines.o
 
 obj/opt_orient.o :  \
@@ -164,4 +186,5 @@ obj/str.o :
 obj/subroutines.o :  \
 	obj/box.o \
 	obj/functions.o \
-	obj/parameters.o
+	obj/parameters.o \
+	obj/str.o

src/atoms.f90

@@ -1162,6 +1162,267 @@ END SELECT
 !
 END SUBROUTINE ATOMMASSSPECIES
 !
+
+subroutine realatomspecies(smass,species)
+!
+IMPLICIT NONE
+CHARACTER(LEN=2),INTENT(OUT):: species
+INTEGER:: mass
+REAL(dp),INTENT(IN):: smass
+
+if(mass_is_equal(smass,1.008d0)) then
+    species = 'H'
+else if(mass_is_equal(smass,2.014101777d0)) then
+      species = 'D'
+else if(mass_is_equal(smass,4.002602d0)) then 
+      species = 'He'
+else if(mass_is_equal(smass,6.94d0)) then
+      species = 'Li'
+else if(mass_is_equal(smass,9.012182d0)) then 
+      species = 'Be'
+else if(mass_is_equal(smass,10.81d0)) then
+      species = 'B'
+else if(mass_is_equal(smass,12.011d0)) then
+      species='C'
+else if(mass_is_equal(smass,14.007d0)) then
+      species='N'
+else if(mass_is_equal(smass,15.999d0)) then 
+      species='O'
+else if(mass_is_equal(smass,18.9984032d0)) then 
+      species='F'
+else if(mass_is_equal(smass,20.1797d0)) then 
+      species='Ne'
 !
+! n=3
+else if(mass_is_equal(smass,22.98976928d0)) then 
+      species = 'Na'
+else if(mass_is_equal(smass,24.305d0)) then 
+      species='Mg'
+else if(mass_is_equal(smass,26.9815386d0)) then 
+      species='Al'
+else if(mass_is_equal(smass,28.085d0)) then 
+      species='Si'
+else if(mass_is_equal(smass,30.973762d0)) then 
+      species='P'
+else if(mass_is_equal(smass,32.06d0)) then 
+      species='S'
+else if(mass_is_equal(smass,35.45d0)) then 
+      species='Cl'
+else if(mass_is_equal(smass,39.948d0)) then 
+      species='Ar'
+!
+! n=4
+else if(mass_is_equal(smass,39.0983d0)) then 
+      species='K'
+else if(mass_is_equal(smass,40.078d0)) then 
+      species='Ca'
+else if(mass_is_equal(smass,44.955912d0)) then 
+      species='Sc'
+else if(mass_is_equal(smass,47.867d0)) then 
+      species='Ti'
+else if(mass_is_equal(smass,50.9415d0)) then 
+      species='V'
+else if(mass_is_equal(smass,51.9961d0)) then 
+      species='Cr'
+else if(mass_is_equal(smass,54.938045d0)) then 
+      species='Mn'
+else if(mass_is_equal(smass,55.845d0)) then 
+      species='Fe'
+else if(mass_is_equal(smass,58.933195d0)) then 
+      species='Co'
+else if(mass_is_equal(smass,58.6934d0)) then 
+      species='Ni'
+else if(mass_is_equal(smass,63.546d0)) then 
+      species='Cu'
+else if(mass_is_equal(smass,65.38d0)) then 
+      species='Zn'
+else if(mass_is_equal(smass,69.723d0)) then 
+      species='Ga'
+else if(mass_is_equal(smass,72.63d0)) then 
+      species='Ge'
+else if(mass_is_equal(smass,74.9216d0)) then 
+      species='As'
+else if(mass_is_equal(smass,78.96d0)) then 
+      species='Se'
+else if(mass_is_equal(smass,79.904d0)) then 
+      species='Br'
+else if(mass_is_equal(smass,83.798d0)) then 
+      species='Kr'
+!
+! n=5
+else if(mass_is_equal(smass,85.4678d0)) then 
+      species='Rb'
+else if(mass_is_equal(smass,87.62d0)) then 
+      species='Sr'
+else if(mass_is_equal(smass,88.90585d0)) then 
+      species='Y'
+else if(mass_is_equal(smass,91.224d0)) then 
+      species='Zr'
+else if(mass_is_equal(smass,92.90638d0)) then 
+      species='Nb'
+else if(mass_is_equal(smass,95.96d0)) then 
+      species='Mo'
+else if(mass_is_equal(smass,98.906d0)) then 
+      species='Tc'
+else if(mass_is_equal(smass,101.07d0)) then 
+      species='Ru'
+else if(mass_is_equal(smass,102.90550d0)) then 
+      species='Rh'
+else if(mass_is_equal(smass,106.42d0)) then 
+      species='Pd'
+else if(mass_is_equal(smass,107.8682d0)) then 
+      species='Ag'
+else if(mass_is_equal(smass,112.411d0)) then 
+      species='Cd'
+else if(mass_is_equal(smass,114.818d0)) then 
+      species='In'
+else if(mass_is_equal(smass,118.71d0)) then 
+      species='Sn'
+else if(mass_is_equal(smass,121.76d0)) then 
+      species='Sb'
+else if(mass_is_equal(smass,127.60d0)) then 
+      species='Te'
+else if(mass_is_equal(smass,126.90447d0)) then 
+      species='I'
+else if(mass_is_equal(smass,131.293d0)) then 
+      species='Xe'
+!
+! n=6
+else if(mass_is_equal(smass,132.9054519d0)) then 
+      species='Cs'
+else if(mass_is_equal(smass,137.327d0)) then 
+      species='Ba'
+! Lanthanides
+else if(mass_is_equal(smass,138.90547d0)) then 
+      species='La'
+else if(mass_is_equal(smass,140.116d0)) then 
+      species='Ce'
+else if(mass_is_equal(smass,140.90765d0)) then 
+      species='Pr'
+else if(mass_is_equal(smass,144.242d0)) then 
+      species='Nd'
+else if(mass_is_equal(smass,144.91d0)) then 
+      species='Pm'
+else if(mass_is_equal(smass,150.36d0)) then 
+      species='Sm'
+else if(mass_is_equal(smass,151.964d0)) then 
+      species='Eu'
+else if(mass_is_equal(smass,157.25d0)) then 
+      species='Gd'
+else if(mass_is_equal(smass,158.92535d0)) then 
+      species='Tb'
+else if(mass_is_equal(smass,162.50d0)) then 
+      species='Dy'
+else if(mass_is_equal(smass,164.93032d0)) then 
+      species='Ho'
+else if(mass_is_equal(smass,167.259d0)) then 
+      species='Er'
+else if(mass_is_equal(smass,168.93421d0)) then 
+      species='Tm'
+else if(mass_is_equal(smass,173.054d0)) then 
+      species='Yb'
+else if(mass_is_equal(smass,174.9668d0)) then 
+      species='Lu'
+! End of Lanthanides
+else if(mass_is_equal(smass,178.49d0)) then 
+      species='Hf'
+else if(mass_is_equal(smass,180.94788d0)) then 
+      species='Ta'
+else if(mass_is_equal(smass,183.84d0)) then 
+      species='W'
+else if(mass_is_equal(smass,186.207d0)) then 
+      species='Re'
+else if(mass_is_equal(smass,190.23d0)) then 
+      species='Os'
+else if(mass_is_equal(smass,192.217d0)) then 
+      species='Ir'
+else if(mass_is_equal(smass,195.084d0)) then 
+      species='Pt'
+else if(mass_is_equal(smass,196.966569d0)) then 
+      species='Au'
+else if(mass_is_equal(smass,200.59d0)) then 
+      species='Hg'
+else if(mass_is_equal(smass,204.38d0)) then 
+      species='Tl'
+else if(mass_is_equal(smass,207.2d0)) then 
+      species='Pb'
+else if(mass_is_equal(smass,208.9804d0)) then 
+      species='Bi'
+else if(mass_is_equal(smass,209.98d0)) then 
+      species='Po'
+else if(mass_is_equal(smass,209.99d0)) then 
+      species='At'
+else if(mass_is_equal(smass,222.02d0)) then 
+      species='Rn'
+!
+! n=7
+else if(mass_is_equal(smass,233.d0)) then 
+      species='Fr'
+else if(mass_is_equal(smass,226.d0)) then 
+      species='Ra'
+! Actinides
+else if(mass_is_equal(smass,227.d0)) then 
+      species='Ac'
+else if(mass_is_equal(smass,232.03806d0)) then 
+      species='Th'
+else if(mass_is_equal(smass,231.03588d0)) then 
+      species='Pa'
+else if(mass_is_equal(smass,238.02891d0)) then 
+      species='U'
+else if(mass_is_equal(smass,237.d0)) then 
+      species='Np'
+else if(mass_is_equal(smass,244.d0)) then 
+      species='Pu'
+else if(mass_is_equal(smass,243.d0)) then 
+      species='Am'
+else if(mass_is_equal(smass,247.d0)) then 
+      species='Cm'
+else if(mass_is_equal(smass,247.d0)) then 
+      species='Bk'
+else if(mass_is_equal(smass,251.d0)) then 
+      species='Cf'
+else if(mass_is_equal(smass,252.d0)) then 
+      species='Es'
+else if(mass_is_equal(smass,257.d0)) then 
+      species='Fm'
+else if(mass_is_equal(smass,258.d0)) then 
+      species='Md'
+else if(mass_is_equal(smass,259.d0)) then 
+      species='No'
+else if(mass_is_equal(smass,262.d0)) then 
+      species='Lr'
+! End of actinides
+else if(mass_is_equal(smass,265.d0)) then 
+      species='Rf'
+else if(mass_is_equal(smass,268.d0)) then 
+      species='Db'
+else if(mass_is_equal(smass,271.d0)) then 
+      species='Sg'
+else if(mass_is_equal(smass,270.d0)) then 
+      species='Bh'
+else if(mass_is_equal(smass,277.d0)) then 
+      species='Hs'
+else if(mass_is_equal(smass,276.d0)) then 
+      species='Mt'
+else if(mass_is_equal(smass,281.d0)) then 
+      species='Ds'
+else if(mass_is_equal(smass,280.d0)) then 
+      species='Rg'
+else if(mass_is_equal(smass,285.17d0)) then 
+      species='Cn'
+else if(mass_is_equal(smass,284.d0)) then 
+      species='Uu'
+else if(mass_is_equal(smass,289.d0)) then 
+      species='Fl'
+else if(mass_is_equal(smass,288.d0)) then 
+      species='Mc'
+else if(mass_is_equal(smass,293.d0)) then 
+      species='Lv'
+else if(mass_is_equal(smass,294.d0)) then 
+      species='Ts'
+else if(mass_is_equal(smass,294.d0)) then 
+      species='Og'
 !
+end if
+end subroutine realatomspecies
 END MODULE atoms

src/box.f90

@@ -5,14 +5,11 @@ module box
     implicit none 
 
     real(kind=dp) :: box_bd(6) !Global box boundaries
+    real(kind=dp) :: box_ori(3,3) !Box orientation
     character(len=3) :: box_bc !Box boundary conditions (periodic or shrinkwrapped)
     logical :: bound_called 
     !The subbox variables contain values for each subbox, being the boxes read in through some
-    !command. Currently only mode_merge will require sub_boxes, for mode_create it will always 
-    !allocate to only 1 sub_box
-    integer :: sub_box_num = 0
-    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
+    !command. 
 
     !Below are some simulation parameters which may be adjusted if reading in restart files
     integer :: timestep=0
@@ -29,42 +26,6 @@ module box
         bound_called=.false.
     end subroutine box_init
 
-    subroutine alloc_sub_box(n)
-        !Allocate the sub_box variables
-
-        integer, intent(in) :: n
-
-        integer :: i
-
-        allocate(sub_box_ori(3,3,n), sub_box_bd(6,n))
-        do i = 1, n
-            sub_box_ori(:,:,i) = identity_mat(3)
-            sub_box_bd(:,i) = 0.0_dp
-        end do
-    end subroutine alloc_sub_box
-
-    subroutine grow_sub_box(n)
-        !Grows sub box arrays, this is only called when a new file is read in
-        integer, intent(in) :: n
-
-        integer, allocatable :: temp_array_bd(:,:,:), temp_file(:)
-        real(kind=dp), allocatable :: temp_ori(:,:,:), temp_bd(:,:)
-        !Allocate temporary arrays
-        allocate(temp_ori(3,3,sub_box_num+n),temp_bd(6,sub_box_num+n), &
-                 temp_array_bd(2,2,sub_box_num+n), temp_file(sub_box_num+n))
-
-        !Move allocation for all sub_box_arrays
-        temp_ori(:,:,1:sub_box_num) = sub_box_ori
-        temp_ori(:,:,sub_box_num+1:) = 0.0_dp
-        call move_alloc(temp_ori, sub_box_ori)
-
-        temp_bd(:, 1:sub_box_num) = sub_box_bd
-        temp_bd(:, sub_box_num+1:) = 0.0_dp
-        call move_alloc(temp_bd, sub_box_bd)
-
-        return
-    end subroutine grow_sub_box
-
     subroutine grow_box(temp_box_bd)
         !This function takes in a temporary box boundary and adjusts the overall box boundaries
         !to include  it
@@ -85,24 +46,6 @@ module box
         return
     end subroutine grow_box
 
-
-    subroutine in_sub_box(r, which_sub_box)
-        !This returns which sub_box a point is in. It returns the first sub_box with boundaries which 
-        !contain the point.
-        real(kind=dp), dimension(3), intent(in) :: r
-        integer, intent(out) :: which_sub_box
-
-        integer :: i 
-
-        do i = 1, sub_box_num
-            if( in_block_bd(r, sub_box_bd(:,i))) then 
-                which_sub_box = i
-                exit
-            end if
-        end do 
-    return 
-    end subroutine in_sub_box
-
     subroutine reset_box
         !This subroutine just resets the box boundary and position
         box_bc = "ppp"
@@ -115,4 +58,5 @@ module box
         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

src/caller.f90

@@ -1,6 +1,7 @@
 module caller
     !this module just calls modes and options
 
+    use mode_da
     use mode_create
     use mode_convert
     use mode_merge
@@ -14,6 +15,7 @@ module caller
     use opt_delete
     use opt_redef_box
     use opt_slip_plane
+    use opt_bubble
     use box
 
 
@@ -37,6 +39,8 @@ module caller
             call metric(arg_pos)
         case('--calc')
             call calc(arg_pos)
+        case('--da')
+            call da(arg_pos)
         case default
             print *, "Mode ", trim(adjustl(mode)), " currently not accepted. Please check documentation for ", &
                      "accepted modes and rerun."
@@ -56,8 +60,14 @@ module caller
             call group(arg_pos)
         case('-ow')
             arg_pos = arg_pos + 1
+        case('-novel')
+            arg_pos = arg_pos+1
         case('-wrap')
             arg_pos = arg_pos + 1
+        case('-norefine')
+            arg_pos = arg_pos + 1
+        case('-all_new')
+            arg_pos = arg_pos + 1
         case('-orient')
             call orient_opt(arg_pos)
         case('-unorient')
@@ -68,8 +78,8 @@ module caller
             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('-box_ori')
+            call set_box_ori(arg_pos)
         case('-deform')
             call deform(arg_pos)
         case('-delete')
@@ -82,6 +92,8 @@ module caller
             call redef_box(arg_pos)
         case('-slip_plane')
             call run_slip_plane(arg_pos)
+        case('-bubble')
+            call bubble(arg_pos)
         case default
             print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.'
             stop 3 

src/elements.f90

@@ -2,6 +2,7 @@ module elements
     
     !This module contains the elements data structures, structures needed for building regions
     !and operations that are done on elements
+    use atoms
     use parameters
     use functions
     use subroutines
@@ -11,10 +12,10 @@ 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(:), tag_ele(:) !Element size
+    integer, allocatable :: size_ele(:), lat_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(:,:,:)
+    real(kind=dp), allocatable :: force_node(:,:,:,:), virial_node(:,:,:,:), energy_node(:,:,:), vel_node(:,:,:,:)
 
     integer, save :: ele_num !Number of elements
     integer, save :: node_num !Total number of nodes
@@ -22,19 +23,21 @@ module elements
 
     !Data structure used to represent atoms 
     integer, allocatable ::  type_atom(:)!atom type 
-    integer, allocatable :: sbox_atom(:), tag_atom(:)
+    integer, allocatable :: 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(:)
+    real(kind=8), allocatable :: force_atom(:,:), virial_atom(:,:), energy_atom(:), vel_atom(:,:)
 
     !Mapping atom type to provided name
     character(len=2), dimension(10) :: type_to_name
     integer :: atom_types = 0
+    real(kind=dp) :: masses(10)
 
     !Variables for creating elements based on primitive cells
-    real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3), bcc_cell(3,8), bcc_mat(3,3), bcc_inv(3,3)
-    integer :: cubic_faces(4,6)
+    real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3), bcc_cell(3,8), bcc_mat(3,3), bcc_inv(3,3), &
+                     cube20(3,20)   
+    integer :: cubic_faces(4,6), oneonetwopairs(2,6)
 
     !Below are lattice type arrays which provide information on the general form of the elements.
     !We currently have a limit of 10 lattice types for simplicities sake but this can be easily increased.
@@ -43,7 +46,7 @@ module elements
     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 
+    !These variables contain informatione on the basis, for simplicities sake we limit 
     !the user to the definition of 10 lattice types with 10 basis atoms at each lattice point.
     !This can be easily increased with no change to efficiency
     integer :: max_basisnum, basisnum(10) !Max basis atom number, number of basis atoms in each lattice type 
@@ -72,13 +75,44 @@ module elements
                               0.0_dp, 1.0_dp, 1.0_dp /), &
                              shape(fcc_cell))
 
+        !Now initialize the 20 node cube element
+        cube20 = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, &
+                             1.0_dp, 0.0_dp, 0.0_dp, &
+                             1.0_dp, 1.0_dp, 0.0_dp, &
+                             0.0_dp, 1.0_dp, 0.0_dp, &
+                             0.0_dp, 0.0_dp, 1.0_dp, &
+                             1.0_dp, 0.0_dp, 1.0_dp, &
+                             1.0_dp, 1.0_dp, 1.0_dp, &
+                             0.0_dp, 1.0_dp, 1.0_dp, &
+                             0.5_dp, 0.0_dp, 0.0_dp, &
+                             1.0_dp, 0.5_dp, 0.0_dp, &
+                             0.5_dp, 1.0_dp, 0.0_dp, &
+                             0.0_dp, 0.5_dp, 0.0_dp, &
+                             0.0_dp, 0.0_dp, 0.5_dp, &
+                             1.0_dp, 0.0_dp, 0.5_dp, &
+                             1.0_dp, 1.0_dp, 0.5_dp, &
+                             0.0_dp, 1.0_dp, 0.5_dp, &
+                             0.5_dp, 0.0_dp, 1.0_dp, &
+                             1.0_dp, 0.5_dp, 1.0_dp, &
+                             0.5_dp, 1.0_dp, 1.0_dp, &
+                             0.0_dp, 0.5_dp, 1.0_dp /), &
+                             shape(cube20))
+
         !Now we create a list containing the list of vertices needed to describe the 6 cube faces
         cubic_faces(:,1) = (/ 1, 2, 3, 4 /) 
         cubic_faces(:,2) = (/ 1, 2, 6, 5 /)
         cubic_faces(:,3) = (/ 2, 3, 7, 6 /)
-        cubic_faces(:,4) = (/ 3, 4, 8, 7 /)
-        cubic_faces(:,5) = (/ 1, 4, 8, 5 /)
+        cubic_faces(:,4) = (/ 1, 4, 8, 5 /)
+        cubic_faces(:,5) = (/ 4, 3, 7, 8 /)
         cubic_faces(:,6) = (/ 5, 6, 7, 8 /)
+        
+        !Now mark which node pairs represent the 112 directions. This is used for the dislocation analysis.
+        oneonetwopairs(:,1) = (/ 1, 3 /)
+        oneonetwopairs(:,2) = (/ 1, 6 /)
+        oneonetwopairs(:,3) = (/ 2, 7 /)
+        oneonetwopairs(:,4) = (/ 1, 8 /)
+        oneonetwopairs(:,5) = (/ 4, 7 /)
+        oneonetwopairs(:,6) = (/ 5, 7 /)
 
         !!Now initialize the fcc primitive cell
         fcc_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, &
@@ -119,6 +153,7 @@ module elements
 
         max_basisnum = 0
         basisnum(:) = 0
+        basis_type(:,:) = 0
         ng_node(:) = 0
         node_num = 0
         node_atoms = 0
@@ -179,7 +214,7 @@ module elements
         
         !Allocate element arrays
         if(n > 0) then 
-            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), &
+            allocate(type_ele(n), tag_ele(n), size_ele(n), lat_ele(n), r_node(3,max_basisnum, max_ng_node,n), &
                        stat=allostat)
             if(allostat > 0) then 
                 print *, "Error allocating element arrays in elements.f90 because of: ", allostat
@@ -189,7 +224,7 @@ module elements
 
         if(m > 0) then 
             !Allocate atom arrays
-            allocate(type_atom(m), sbox_atom(m), tag_atom(m), r_atom(3,m), stat=allostat)
+            allocate(type_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
@@ -233,11 +268,6 @@ module elements
                 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)
@@ -265,11 +295,6 @@ module elements
                 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
@@ -280,11 +305,11 @@ module elements
         end if
     end subroutine
 
-    subroutine add_element(tag, type, size, lat, sbox, r)
+    subroutine add_element(tag, type, size, lat, r)
         !Subroutine which adds an element to the element arrays
-        integer, intent(in) :: size, lat, sbox, tag
+        integer, intent(in) :: size, lat, tag
         character(len=100), intent(in) :: type
-        real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node)
+        real(kind=dp), intent(in) :: r(:,:,:)
 
         integer :: newtag
 
@@ -304,15 +329,14 @@ module elements
         type_ele(ele_num) = type
         size_ele(ele_num) = size
         lat_ele(ele_num) = lat
-        sbox_ele(ele_num) = sbox
         r_node(:,:,:,ele_num) = r(:,:,:)
 
 
     end subroutine add_element
 
-    subroutine add_atom(tag, type, sbox, r)
+    subroutine add_atom(tag, type, r)
         !Subroutine which adds an atom to the atom arrays
-        integer, intent(in) :: type, sbox, tag
+        integer, intent(in) :: type,  tag
         real(kind=dp), intent(in), dimension(3) :: r
 
         integer :: newtag
@@ -325,29 +349,38 @@ module elements
         end if
         !Check to see if we need to grow the arrays
         call grow_ele_arrays(0,1)
-        tag_atom(atom_num) = tag
+        tag_atom(atom_num) = newtag
         type_atom(atom_num) = type
         r_atom(:,atom_num) = r(:)
-        sbox_atom(atom_num) = sbox
 
     end subroutine add_atom
 
-    subroutine add_atom_type(type, inttype)
+    subroutine add_atom_type(mass, inttype, force_new)
         !This subroutine adds a new atom type to the module list of atom types
-        character(len=2), intent(in) :: type
+        real(kind=dp), intent(in) :: mass
         integer, intent(out) :: inttype
 
+        logical, optional, intent(in) :: force_new
+
         integer :: i
-        logical :: exists
+        logical :: exists, force
         
+        if(present(force_new)) then 
+            force = force_new
+        else 
+            force = .false.
+        end if
+
         exists = .false.
-        do i=1, 10
-            if(type == type_to_name(i)) then 
-                exists = .true.
-                inttype = i
-                exit
-            end if
-        end do 
+        if(.not.force) then 
+            do i=1, 10
+                if(is_equal(mass, masses(i))) then 
+                    exists = .true.
+                    inttype = i
+                    exit
+                end if
+            end do 
+        end if
 
         if (exists.eqv..false.) then 
             atom_types = atom_types+1
@@ -355,7 +388,8 @@ module elements
                 print *, "Defined atom types are greater than 10 which is currently not supported."
                 stop 3
             end if
-            type_to_name(atom_types) = type
+            call atommassspecies(mass, type_to_name(atom_types))
+            masses(atom_types) = mass
             inttype = atom_types
         end if
         return
@@ -374,6 +408,8 @@ module elements
                 ng_node(i) = 8
             case('bcc')
                 ng_node(i) = 8
+            case('20fcc')
+                ng_node(i) = 20
             end select
 
             if(ng_node(i) > max_ng_node) max_ng_node = ng_node(i)
@@ -393,18 +429,18 @@ module elements
         !This subroutine returns the interpolated atoms from the elements. 
 
         !Arguments
-        character(len=100), intent(in) :: type !The type of element that it is
+        character(len=*), intent(in) :: type !The type of element that it is
         integer, intent(in) :: esize !The number of atoms per side
         integer, intent(in) :: lat_type !The integer lattice type of the element
-        real(kind=dp), dimension(3,max_basisnum, max_ng_node), intent(in) :: r_in !Nodal positions
-        integer, dimension(max_basisnum*max_esize**3), intent(out) :: type_interp !Interpolated atomic positions
+        real(kind=dp), dimension(:,:,:), intent(in) :: r_in !Nodal positions
+        integer, dimension(:), 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
+        real(kind = dp), optional, intent(in) :: eng(:,:), f(:,:,:), &
+                                                 v(:,:,:)
+        real(kind=dp), dimension(:,:), optional,intent(out) :: data_interp !Interpolated atomic positions
 
         !Internal variables
-        integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp
+        integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp, adjust
         real(kind=dp) :: a_shape(max_ng_node)
         real(kind=dp) :: t, s, r
         
@@ -426,8 +462,36 @@ module elements
             lat_type_temp = lat_type
         end select
 
-        select case(trim(adjustl(type)))
+        select case(type)
         case('fcc','bcc')
+            !Now loop over all the possible sites
+            do it = 1, esize
+                t=-1.0_dp +(it-1)*(2.0_dp/(esize-1))   
+                do is =1, esize
+                    s=-1.0_dp +(is-1)*(2.0_dp/(esize-1))   
+                    do ir = 1, esize
+                        r=-1.0_dp +(ir-1)*(2.0_dp/(esize-1))   
+                        call rhombshape(r,s,t,a_shape)
+
+                        do ibasis = 1, bnum
+                            ia = ia + 1
+                            do inod = 1, 8
+                                type_interp(ia) = basis_type(ibasis,lat_type_temp)
+                                r_interp(:,ia) = r_interp(:,ia) + a_shape(inod) * r_in(:,ibasis,inod)
+
+                                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
+
+                    end do
+                end do
+            end do 
+        case('20fcc')
             !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)
@@ -435,11 +499,11 @@ module elements
                     s = (1.0_dp*(is-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
                     do ir = 1, esize
                         r = (1.0_dp*(ir-1) - (esize-1)/2)/(1.0_dp*(esize-1)/2)
-                        call rhombshape(r,s,t,a_shape)
+                        call quad_rhombshape(r,s,t,a_shape)
 
                         do ibasis = 1, bnum
                             ia = ia + 1
-                            do inod = 1, 8
+                            do inod = 1, 20
                                 type_interp(ia) = basis_type(ibasis,lat_type_temp)
                                 r_interp(:,ia) = r_interp(:,ia) + a_shape(inod) * r_in(:,ibasis,inod)
 
@@ -464,6 +528,68 @@ module elements
         return
     end subroutine interpolate_atoms
 
+    subroutine interpolate_vel(type, esize, lat_type, vel_in, vel_interp)
+        !This subroutine returns the interpolated atoms from the elements. 
+
+        !Arguments
+        character(len=*), intent(in) :: type !The type of element that it is
+        integer, intent(in) :: esize !The number of atoms per side
+        integer, intent(in) :: lat_type !The integer lattice type of the element
+        real(kind=dp), dimension(:,:,:), intent(in) :: vel_in !Nodal positions
+        real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: vel_interp !Interpolated atomic positions
+
+        !Internal variables
+        integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp, adjust
+        real(kind=dp) :: a_shape(max_ng_node)
+        real(kind=dp) :: t, s, r
+        
+        !Initialize some variables
+        vel_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
+        !the basis is 0,0,0, and the type doesn't matter
+
+        select case(lat_type)
+        case(0)
+            bnum=1
+            lat_type_temp = 1
+        case default
+            bnum = basisnum(lat_type)
+            lat_type_temp = lat_type
+        end select
+
+        select case(type)
+        case('fcc','bcc')
+            !Now loop over all the possible sites
+            do it = 1, esize
+                t=-1.0_dp +(it-1)*(2.0_dp/(esize-1))   
+                do is =1, esize
+                    s=-1.0_dp +(is-1)*(2.0_dp/(esize-1))   
+                    do ir = 1, esize
+                        r=-1.0_dp +(ir-1)*(2.0_dp/(esize-1))   
+                        call rhombshape(r,s,t,a_shape)
+
+                        do ibasis = 1, bnum
+                            ia = ia + 1
+                            do inod = 1, 8
+                                vel_interp(:,ia) = vel_interp(:,ia) + a_shape(inod) * vel_in(:,ibasis,inod)
+
+                            end do
+                        end do
+
+                    end do
+                end do
+            end do 
+        end select
+        if (ia /= bnum*esize**3) then 
+            print *, "Incorrect interpolation"
+            stop 3
+        end if
+        return
+
+    end subroutine interpolate_vel
+
     subroutine rhombshape(r,s,t, shape_fun)
         !Shape function for rhombohedral elements
         real(kind=8), intent(in) :: r, s, t
@@ -482,6 +608,42 @@ module elements
         return
     end subroutine rhombshape
 
+    subroutine quad_rhombshape(r,s,t, shape_fun)
+        !Shape function for 20 node quadratic rhombohedral elements
+        real(kind=8), intent(in) :: r, s, t
+        real(kind=8), intent(out) :: shape_fun(:)
+
+        !Corner nodes
+        shape_fun(1) = (1.0_dp-r)*(1.0_dp-s)*(1.0_dp-t)*(-r-s-t-2)/8.0_dp
+        shape_fun(2) = (1.0_dp+r)*(1.0_dp-s)*(1.0_dp-t)*(r-s-t-2)/8.0_dp
+        shape_fun(3) = (1.0_dp+r)*(1.0_dp+s)*(1.0_dp-t)*(r+s-t-2)/8.0_dp
+        shape_fun(4) = (1.0_dp-r)*(1.0_dp+s)*(1.0_dp-t)*(-r+s-t-2)/8.0_dp
+        shape_fun(5) = (1.0_dp-r)*(1.0_dp-s)*(1.0_dp+t)*(-r-s+t-2)/8.0_dp
+        shape_fun(6) = (1.0_dp+r)*(1.0_dp-s)*(1.0_dp+t)*(r-s+t-2)/8.0_dp
+        shape_fun(7) = (1.0_dp+r)*(1.0_dp+s)*(1.0_dp+t)*(r+s+t-2)/8.0_dp
+        shape_fun(8) = (1.0_dp-r)*(1.0_dp+s)*(1.0_dp+t)*(-r+s+t-2)/8.0_dp
+
+        !Side nodes, first node r is zero
+        shape_fun(9) = (1-r*r)*(1-s)*(1-t)/4.0_dp
+        shape_fun(11) = (1-r*r)*(1+s)*(1-t)/4.0_dp
+        shape_fun(17) = (1-r*r)*(1-s)*(1+t)/4.0_dp
+        shape_fun(19) = (1-r*r)*(1+s)*(1+t)/4.0_dp
+
+        !node s is zero
+        shape_fun(10) = (1+r)*(1-s*s)*(1-t)/4.0_dp
+        shape_fun(12) = (1-r)*(1-s*s)*(1-t)/4.0_dp
+        shape_fun(18) = (1+r)*(1-s*s)*(1+t)/4.0_dp
+        shape_fun(20) = (1-r)*(1-s*s)*(1+t)/4.0_dp
+
+        !node t is zero
+        shape_fun(13) = (1-r)*(1-s)*(1-t*t)/4.0_dp
+        shape_fun(14) = (1+r)*(1-s)*(1-t*t)/4.0_dp
+        shape_fun(15) = (1+r)*(1+s)*(1-t*t)/4.0_dp
+        shape_fun(16) = (1-r)*(1+s)*(1-t*t)/4.0_dp
+
+
+    end subroutine quad_rhombshape
+
     subroutine delete_atoms(num, in_index)
         !This subroutine deletes atoms from the atom arrays 
         integer, intent(in) :: num
@@ -533,19 +695,18 @@ module elements
         !accidentally overwrite values which need to be deleted
         do i = num, 1, -1
             node_num = node_num - ng_node(lat_ele(sorted_index(i)))
+            node_atoms = node_atoms - ng_node(lat_ele(sorted_index(i)))*basisnum(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
                 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
@@ -555,12 +716,47 @@ module elements
 
     subroutine wrap_atoms
         !This subroutine wraps atoms back into the simulation cell if they have exited for any reason
+
         integer :: i
         do i = 1, atom_num
             call apply_periodic(r_atom(:,i))
         end do 
+
     end subroutine wrap_atoms
 
+    subroutine wrap_elements
+        integer :: i, inod, ibasis, j, node_in_bd(3,max_ng_node)
+        real(kind =dp) :: box_len(3)
+
+        do j = 1, 3
+            box_len(j) = box_bd(2*j) - box_bd(2*j-1)
+        end do
+        print *, box_bd
+        do i = 1, ele_num
+            node_in_bd = 0
+            do inod = 1, ng_node(lat_ele(i))
+                do ibasis = 1, basisnum(lat_ele(i))
+                    do j = 1, 3
+                        if(r_node(j,ibasis,inod,i) < box_bd(2*j-1)) then 
+                            node_in_bd(j,inod) = 1
+                        else if(r_node(j,ibasis,inod,i) > box_bd(2*j)) then 
+                            node_in_bd(j,inod) = -1
+                        end if
+                    end do
+                end do
+            end do
+
+            do j = 1, 3
+                if(all(node_in_bd(j,:) == 1)) then 
+                    r_node(j,:,:,i) = r_node(j,:,:,i) + box_len(j)
+                else if(all(node_in_bd(j,:) == -1)) then 
+                    r_node(j,:,:,i) = r_node(j,:,:,i) - box_len(j)
+                end if
+            end do
+        end do
+    end subroutine wrap_elements
+
+
     subroutine def_new_box
         !This subroutine calculates new box boundaries based on minimum and maximum nodal/atomic positions
         integer :: i, j, inod, ibasis
@@ -613,6 +809,9 @@ do i = 1, atom_num
         else if(trim(adjustl(pos_string)) == '-inf') then
             pos=box_bd(2*i-1)
 
+        else if(trim(adjustl(pos_string)) == 'mid') then 
+            pos=(box_bd(2*i)+box_bd(2*i-1))/2
+
         else if (trim(adjustl(pos_string)) == 'rand') then 
             call random_number(rand)
             pos = (box_bd(2*i)-box_bd(2*i-1))*rand + box_bd(2*i-1)
@@ -636,7 +835,6 @@ do i = 1, atom_num
                 end do 
             end if
 
-            read(cone, *) face
             if ((face < 1).or.(face > 6)) stop "Current face number must be 1,2,3,4,5,6. Please check documentation" 
             !Now get the position
             call offset_pos(ele, face, rand_ele_pos)
@@ -659,6 +857,15 @@ do i = 1, atom_num
             end if
             pos = box_bd(2*i) - pos
 
+        else if ((index(pos_string,'-') > 0).and.(index(pos_string,'mid')>0)) then 
+            !Now extract the number we are reducing from midinity
+            if(index(pos_string,'mid') < index(pos_string,'-')) then 
+                read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos
+            else
+                read(pos_string(1:index(pos_string,'-')-1), *, iostat=iospara) pos
+            end if
+            pos = (box_bd(2*i)+box_bd(2*i-1))/2.0_dp - pos
+
         else if ((index(pos_string,'+') > 0).and.(index(pos_string,'inf')>0)) then 
             !Now extract the number we are reducing from infinity
            if(index(pos_string,'inf') < index(pos_string,'+')) then 
@@ -668,6 +875,15 @@ do i = 1, atom_num
             end if
             pos = box_bd(2*i-1) + pos
 
+        else if ((index(pos_string,'+') > 0).and.(index(pos_string,'mid')>0)) then 
+            !Now extract the number we are reducing from midinity
+           if(index(pos_string,'mid') < index(pos_string,'+')) then 
+                read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos
+            else
+                read(pos_string(1:index(pos_string,'+')-1), *, iostat=iospara) pos
+            end if
+            pos = (box_bd(2*i)+box_bd(2*i-1))/2.0_dp + pos
+
         else if ((index(pos_string,'*') > 0).and.(index(pos_string,'inf')>0)) then 
             !Now extract the number we are reducing from infinity
            if(index(pos_string,'inf') < index(pos_string,'*')) then 
@@ -721,11 +937,11 @@ do i = 1, atom_num
 
     case('fcc')
         !First we have to extract the element lattice parameter
-        call matrix_inverse(sub_box_ori(:,:,sbox_ele(ie)),3,ori_inv)
+        call matrix_inverse(box_ori(:,:),3,ori_inv)
         r_cubic_node = r_node(:,1,:,ie)
         r_cubic_node = matmul(fcc_inv,matmul(ori_inv,r_cubic_node))
         lp = (maxval(r_cubic_node(1,:))-minval(r_cubic_node(1,:)))/(esize-1)
-        pos = matmul(sub_box_ori(:,:,sbox_ele(ie)),matmul(fcc_mat,pos))*lp
+        pos = matmul(box_ori,matmul(fcc_mat,pos))*lp
         pos = pos + r_node(:,1, 1, ie)
 
     case default
@@ -830,14 +1046,71 @@ do i = 1, atom_num
 
     end subroutine 
 
+    subroutine alloc_vel_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(vel_node)) then 
+                deallocate(vel_node)
+            end if 
+            allocate(vel_node(3,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(vel_atom)) then 
+                deallocate(vel_atom)
+            end if 
+            allocate(vel_atom(3,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 alloc_vel_arrays
+
+
     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)
+        real(kind=dp), intent(in) :: eng, force(:), virial(:)
+
+        integer :: size_atom_array
+        real(kind=dp), allocatable :: tmp_eng(:), tmp_force(:,:), tmp_virial(:,:)
+
+        size_atom_array=size(tag_atom)
+        if(ia > atom_num) then 
+            stop "ia in add_atom_dat shouldn't be larger than atom_num"
+        else if(ia > size(energy_atom)) then
+            allocate(tmp_eng(size(energy_atom)+1024))
+            allocate(tmp_force(3, size(energy_atom)+1024)) 
+            allocate(tmp_virial(6, size(energy_atom)+1024))
+
+            tmp_eng=0.0_dp
+            tmp_eng(1:size(energy_atom)) = energy_atom
+            call move_alloc(tmp_eng, energy_atom)
+
+            tmp_force=0.0_dp
+            tmp_force(:, 1:size(force_atom, 2)) = force_atom
+            call move_alloc(tmp_force, force_atom)
+
+            tmp_virial=0.0_dp
+            tmp_virial(:, 1:size(virial_atom,2)) = virial_atom
+            call move_alloc(tmp_virial, virial_atom)
+        end if
 
         energy_atom(ia) = eng
         force_atom(:,ia) = force(:) 
         virial_atom(:,ia) = virial(:) 
+
         vflag = .true.
         return
     end subroutine add_atom_data
@@ -848,10 +1121,12 @@ do i = 1, atom_num
         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)
+        integer :: inod
         vflag = .true.
         energy_node(:,:,ie)   = eng
         force_node(:,:,:,ie)  = force
         virial_node(:,:,:,ie) = virial
+        
         return
     end subroutine add_element_data
 
@@ -862,4 +1137,27 @@ do i = 1, atom_num
         node_num = 0
     end subroutine reset_data
 
+    function ele_in_bounds(bd, etype, esize, lat_type, r_in )
+        real(kind=dp), intent(in) :: bd(6), r_in(3, max_basisnum, max_ng_node)
+        integer, intent(in) ::  lat_type, esize
+        character(len=*), intent(in) :: etype
+
+        integer :: i
+        real(kind=dp) :: rinterp(3, max_basisnum*max_esize**3) 
+        integer :: tinterp(max_basisnum*max_esize**3) 
+        logical :: ele_in_bounds
+
+        ele_in_bounds=.false.
+        call interpolate_atoms(etype, esize, lat_type, r_in, tinterp, rinterp) 
+
+        do i=1, esize**3
+            if(in_block_bd(rinterp(:,i), bd))then 
+                ele_in_bounds = .true.
+                exit
+            end if
+        end do
+
+        return
+    end function ele_in_bounds
+
 end module elements

src/functions.f90

@@ -168,11 +168,11 @@ END FUNCTION StrDnCase
 
         do i =1 ,3
             !Check upper bound
-            if(v(i) > (box_bd(2*i)+10.0_dp**(-10)) ) then 
+            if(v(i) > (box_bd(2*i))) then 
                 in_block_bd =.false.
                 exit
             !Check lower bound
-            else if (v(i) < (box_bd(2*i-1)-10.0_dp**(-10))) then 
+            else if (v(i) < (box_bd(2*i-1))+1d-6) then 
                 in_block_bd = .false.
                 exit
             end if
@@ -254,6 +254,19 @@ END FUNCTION StrDnCase
     return
     end function is_equal
 
+    function mass_is_equal(A, B)
+        !This function checks if too numbers are equal within a tolerance
+        real(kind=dp), intent(in) :: A, B
+        logical :: mass_is_equal
+
+        if((A>(B - 10.0_dp**(-2))).and.(A < (B+10.0_dp**(-2)))) then 
+            mass_is_equal = .true.
+        else
+            mass_is_equal = .false.
+        end if
+    return
+    end function mass_is_equal
+
     pure function unitvec(n,vec)
         integer, intent(in) :: n
         real(kind=dp), intent(in) :: vec(n)

src/main.f90

@@ -38,9 +38,11 @@ program main
     !Call initialization functions
     call lattice_init
     call box_init
-    call random_seed
+    call init_random_seed
     force_overwrite=.false.
+    norefine=.false.
     wrap_flag = .false.
+    all_new = .false.
 
     end_mode_arg = 0
     
@@ -61,11 +63,16 @@ program main
         !This lets us know if we need to wrap atomic positions back into the cell
         case('-wrap')
             wrap_flag=.true.
+
+        case('-norefine')
+            norefine=.true.
         
         case('-set_cac')
             call set_cac(i)
         case('-set_types')
             call set_types(i)
+        case('-all_new')
+            all_new=.true.
         end select 
     end do 
     !Determine if a mode is being used and what it is. The first argument has to be the mode
@@ -106,14 +113,17 @@ program main
     end do
 
     !If wrap flag was passed then call the wrap atoms command
-    if(wrap_flag) call wrap_atoms
+    if(wrap_flag) then 
+        call wrap_atoms
+        call wrap_elements
+    end if
 
     !If we called the boundary command then we adjust the box bounds
     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 ((trim(adjustl(mode)) /= "--metric").and.(trim(adjustl(mode)) /= "--calc"))then    
+    if ((trim(adjustl(mode)) /= "--metric").and.(trim(adjustl(mode)) /= "--calc").and.(trim(adjustl(mode)) /= "--da"))then    
         if ((outfilenum == 0)) then
             argument = 'none'
             call get_out_file(argument)
@@ -121,4 +131,5 @@ program main
         call write_out
     end if
 
+    return
 end program main

src/mode_calc.f90

@@ -25,6 +25,9 @@ module mode_calc
         case('tot_virial')
             allocate(calculated(6))
             call calc_tot_virial
+        case('tot_energy')
+            allocate(calculated(1))
+            call calc_tot_energy
         case default
             print *, trim(adjustl(calc_opt)), " is not accepted as a calc option in mode_calc"
             stop 3
@@ -67,23 +70,19 @@ module mode_calc
         calculated = 0
         do i = 1, atom_num
             do j = 1, 6
-                calculated(j) = calculated(j) + virial_atom(j, i)
+                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)))
+                        calculated(j) = calculated(j) - virial_node(j,ibasis,inod,i)*(size_ele(i)**3)/ng_node(lat_ele(i))
                     end do
-                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
@@ -91,5 +90,26 @@ module mode_calc
 
         print *, "Total virial is calculated as : (v11, v22, v33, v32, v31, v21)"
         print *, calculated
-    end subroutine 
+    end subroutine calc_tot_virial
+
+    subroutine calc_tot_energy
+        integer :: i, inod, ibasis, j
+        calculated=0
+        !Sum atom energies
+        do i = 1, atom_num
+            calculated(1) = calculated(1) + energy_atom(i)
+        end do
+
+        !Sum element energies
+        do i =1, ele_num
+            do inod=1, ng_node(lat_ele(i))
+                do ibasis=1, basisnum(lat_ele(i))
+                    calculated(1) = calculated(1) + energy_node(ibasis, inod, i)*(size_ele(i)**3)/ng_node(lat_ele(i))
+                end do
+            end do
+        end do
+
+        print *, 'Total energy in eV is:'
+        print *, calculated(1)
+    end subroutine calc_tot_energy
 end module mode_calc

src/mode_create.f90

@@ -12,10 +12,10 @@ module mode_create
 
     character(len=100) :: name, element_type    
     real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), &
-                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3)
+                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3),  basis_pos(3,10)
     integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), &
-               basis_pos(3,10), esize_nums, esize_index(10)
-    logical :: dup_flag, dim_flag, efill
+               esize_nums, esize_index(10)
+    logical :: dup_flag, dim_flag, efill, crossb(3)
 
     real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:)
     integer, allocatable :: elat(:)
@@ -67,7 +67,7 @@ module mode_create
             end do
             !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')
+            case('fcc', '20fcc')
                 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')
@@ -79,10 +79,10 @@ module mode_create
             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)
+                box_len(i) = box_bd(2*i) - box_bd(2*i-1)
             end do
 
 
-
         else if(dim_flag) then 
             !As a note everything is defined so that the lattice parameter is multiplied in at the end 
             !so we have to divide all the real Angstroms units by the lattice parameter
@@ -93,7 +93,7 @@ module mode_create
             end do
             !Now get the rotated box vertex positions in lattice space. Should be integer units
             select case(trim(adjustl(element_type)))
-            case('fcc')
+            case('fcc', '20fcc')
                 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
@@ -122,7 +122,7 @@ module mode_create
                 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(0,element_type, esize, 1, 1, r_node_temp)
+            call add_element(0,element_type, esize, 1, r_node_temp)
         end if
         
         !If we passed the dup_flag or dim_flag then we have to convert the lattice points and add them to the atom/element arrays
@@ -131,9 +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)
+                call build_with_rhomb(box_lat_vert, fcc_mat, 8, fcc_inv)
             case('bcc')
-                call build_with_rhomb(box_lat_vert, bcc_mat)
+                call build_with_rhomb(box_lat_vert, bcc_mat, 8, bcc_inv)
+            case('20fcc')
+                call build_with_rhomb(box_lat_vert, fcc_mat, 20, fcc_inv)
             case default
                 print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ",&
                          "element type"
@@ -150,7 +152,7 @@ module mode_create
             if(lat_atom_num > 0) then 
                 do i = 1, lat_atom_num
                     do ibasis = 1, basisnum(1)
-                        call add_atom(0,basis_type(ibasis, 1), 1, (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis))
+                        call add_atom(0,basis_type(ibasis, 1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis))
                     end do
                 end do
                 deallocate(r_atom_lat)
@@ -164,16 +166,16 @@ module mode_create
                         end do
                     end do 
 
-                    call add_element(0,element_type, elat(i), 1, 1, r_node_temp)
+                    call add_element(0,element_type, elat(i), 1, r_node_temp)
                 end do 
             end if
         end if
 
-        !The last thing we do is setup the sub_box_boundaries
-        call alloc_sub_box(1)
-        sub_box_num = 1
-        sub_box_ori(:,:,1) = orient
-        sub_box_bd(:,1) = box_bd
+        !If any elements are fully outside the box then wrap them back in
+        if (any(crossb)) then 
+            call wrap_elements
+        end if
+        box_ori = orient
     end subroutine create
     !This subroutine parses the command and pulls out information needed for mode_create
     subroutine parse_command(arg_pos)
@@ -181,9 +183,10 @@ module mode_create
         integer, intent(out) :: arg_pos
         integer ::  ori_pos, i, j, arglen, stat
         character(len=100) :: textholder
-        character(len=20) :: orient_string
+        character(len=100) :: orient_string
         character(len=2) :: btype
-        logical :: isortho, isrighthanded
+        logical :: isortho, isrighthanded, bool
+        real(kind=dp) :: mass
 
 
         !Pull out all required positional arguments
@@ -201,6 +204,7 @@ module mode_create
         call get_command_argument(5, textholder, arglen)
         if(arglen==0) STOP "Esize missing in mode create"
         read(textholder, *, iostat=stat) esize
+        max_esize = esize
         if(stat > 0) STOP "Error reading esize"
 
         arg_pos = 6 
@@ -248,6 +252,12 @@ module mode_create
                     read(textholder, *) origin(i)
                     arg_pos = arg_pos + 1
                 end do
+            case('crossb')
+                do i = 1, 3
+                    call get_command_argument(arg_pos, textholder)
+                    read(textholder, *) crossb(i) 
+                    arg_pos = arg_pos + 1
+                end do
             case('basis')
                 call get_command_argument(arg_pos, textholder)
                 read(textholder, *) basisnum(1)
@@ -256,7 +266,8 @@ module mode_create
                 do i = 1, max_basisnum
                     call get_command_argument(arg_pos, btype)
                     arg_pos = arg_pos+1
-                    call add_atom_type(btype, basis_type(i,1))
+                    call atommass(btype, mass)
+                    call add_atom_type(mass, basis_type(i,1))
                     do j = 1, 3
                         call get_command_argument(arg_pos, textholder)
                         read(textholder, *) basis_pos(j,i)
@@ -273,14 +284,16 @@ module mode_create
             end select
         end do
 
+
         !Calculate the lattice periodicity length in lattice units
         do i = 1, 3
             lattice_space(i) = norm2(orient(i,:))
         end do
-
+        !Now normalize the orientation matrix
+        orient =  matrix_normal(orient,3)
         !Check special periodicity relations
         select case(trim(adjustl(element_type)))
-        case('fcc')
+        case('fcc', '20fcc')
             do i = 1,3
                 !Check if one of the directions is 110
                 if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(orient(i,3),0.0_dp))).or.&
@@ -307,8 +320,6 @@ module mode_create
                 end if
             end do
         end select 
-        !Now normalize the orientation matrix
-        orient =  matrix_normal(orient,3)
         !Now check these to make sure they are right handed and orthogonal
         call check_right_ortho(orient, isortho, isrighthanded)
         if (.not.isortho) then 
@@ -325,25 +336,28 @@ module mode_create
         if (basisnum(1) == 0) then
             max_basisnum = 1
             basisnum(1) = 1
-            call add_atom_type(name, basis_type(1,1)) !If basis command not defined then we use name as the atom_name
+            call atommass(name, mass)
+            call add_atom_type(mass, basis_type(1,1)) !If basis command not defined then we use name as the atom_name
             basis_pos(:,1) = 0.0_dp
         end if
-        !
+        
     end subroutine
     
-    subroutine build_with_rhomb(box_in_lat, transform_matrix)
+    subroutine build_with_rhomb(box_in_lat, transform_matrix, nn, transform_inverse)
         !This subroutine returns all the lattice points in the box in r_lat
 
         !Inputs
         integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space
         real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space
+        real(kind=dp), dimension(3,3), intent(in) :: transform_inverse !The inverse transform
+        integer, intent(in) :: nn
         !Internal variables
-        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)
+        integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz,  numlatpoints, ele(3,nn), rzero(3), efill_size, &
+                   vlat(3), temp_lat(3,nn), m, n, o, j, k, nump_ele, efill_temp_lat(3,nn), filzero(3), &
+                   bd_ele_lat(6), efill_ele(3,nn), ebd(6), shift_vec(3), type_interp(max_basisnum*max_esize**3)
+        real(kind=dp) :: v(3), temp_nodes(3,1,nn), r(3), centroid_bd(6), vreal(3), r_interp(3, max_basisnum*max_esize**3)
         logical, allocatable :: lat_points(:,:,:)
-        logical :: node_in_bd(8), add, lat_points_ele(esize,esize,esize), intersect_bd(3)
+        logical :: node_in_bd(nn), add, lat_points_ele(esize,esize,esize), intersect_bd(3)
 
         !Do some value initialization
         max_esize = esize
@@ -360,6 +374,8 @@ module mode_create
             numlatpoints = numlatpoints*(bd_in_lat(2*i)-bd_in_lat(2*i-1))
         end do
 
+        if(numlatpoints < 0) stop "number of lattice points is negative, this can occur if the model is too big"
+
         !Allocate the correct lat variable
         select case(esize)
         !Atomistics
@@ -393,7 +409,12 @@ module mode_create
             !If we are working with elements we have to use more complex code
 
             !Initialize finite element 
-            ele(:,:) = (esize-1) * cubic_cell(:,:)
+            select case(element_type)
+            case('fcc','bcc')
+                ele(:,:) = (esize-1) * cubic_cell(:,:)
+            case('20fcc')
+                ele(:,:) = (esize-1)*cube20
+            end select
             
             !Make a 3 dimensional array of lattice points. This array is indexed by the integer lattice position. 
             !A value of true means that the coordinate is a lattice point which is within the box_bd
@@ -414,6 +435,7 @@ module mode_create
                 end do
             end do
 
+
             !Now we redefine bd_in_lat The first 3 indices contains limits for the lat_points array
             bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10
             bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10
@@ -423,7 +445,7 @@ module mode_create
                 do iy = 1, bd_in_array(2)
                     do ix = 1, bd_in_array(1)
                         node_in_bd(:) = .false.
-                        do inod = 1, 8
+                        do inod = 1, nn
                             vlat = ele(:,inod) + (/ ix, iy, iz /)
 
                             !Check to see if the node resides at a position containing a lattice point within the box
@@ -445,16 +467,22 @@ module mode_create
             !Now build the finite element region
             lat_ele_num = 0
             lat_atom_num = 0
-            allocate(r_lat(3,8,numlatpoints/esize), elat(numlatpoints/esize))
+            allocate(r_lat(3,nn,numlatpoints/esize), elat(numlatpoints/esize))
 
-            ele(:,:) = (esize-1) * cubic_cell(:,:)
+            select case(element_type)
+            case('fcc','bcc')
+                ele(:,:) = (esize-1) * cubic_cell(:,:)
+            case('20fcc')
+                ele(:,:) = (esize-1)*cube20
+            end select
+ 
             !Redefined the second 3 indices as the number of elements that fit in the bounds
             do i = 1, 3
                 bd_in_array(3+i) = bd_in_array(i)/esize
             end do
 
             !Now start the element at rzero
-            do inod=1, 8
+            do inod=1, nn
                 ele(:,inod) = ele(:,inod) + rzero
             end do
             do iz = -bd_in_array(6), bd_in_array(6)
@@ -463,7 +491,7 @@ module mode_create
                         node_in_bd(:) = .false.
                         temp_nodes(:,:,:) = 0.0_dp
                         temp_lat(:,:) = 0
-                        do inod = 1, 8
+                        do inod = 1, nn
                             vlat= ele(:,inod) + (/ ix*(esize), iy*(esize), iz*(esize) /)
                             !Transform point back to real space for easier checking
                             ! v = matmul(orient, matmul(transform_matrix,v))
@@ -479,6 +507,31 @@ module mode_create
                             !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.
+                            else if(any(crossb)) then 
+                                vreal=0
+                                do i = 1, 3
+                                    if(crossb(i)) then 
+                                        if(temp_nodes(i,1,inod) < box_bd(2*i-1)) then 
+                                            vreal(i) = temp_nodes(i,1,inod)+box_len(i)
+                                        else if(temp_nodes(i,1,inod) > box_bd(2*i)) then
+                                            vreal(i) = temp_nodes(i,1,inod)-box_len(i)
+                                        else 
+                                            vreal(i) = temp_nodes(i,1,inod)
+                                        end if
+                                    else
+                                        vreal(i) = temp_nodes(i,1,inod)
+                                    end if
+                                end do
+                                v = matmul(transform_inverse, matmul(orient_inv, vreal))
+                                do i = 1, 3
+                                    vlat(i) = nint(v(i) - bd_in_lat(2*i-1)+5)
+                                end do
+                                if(any(vlat > shape(lat_points)).or.any(vlat < 1)) then 
+                                    continue 
+                                !If within array boundaries check to see if it is a lattice point
+                                else if(lat_points(vlat(1),vlat(2),vlat(3))) then 
+                                    node_in_bd(inod) = .true.
+                                end if
                             end if
                         end do
 
@@ -488,14 +541,50 @@ module mode_create
                             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,:))
-                                    do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:))
-                                        lat_points(m,n,o) = .false.
-                                    end do 
+                            if(any(crossb)) then 
+                                call interpolate_atoms('fcc', esize, 0, temp_nodes, type_interp, r_interp)
+                                j= 0
+                                 do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))
+                                    do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:))
+                                        do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:))
+                                            j=j+1
+                                            do i = 1, 3
+                                                if(crossb(i)) then 
+                                                    if(r_interp(i,j) < box_bd(2*i-1)) then 
+                                                        vreal(i) = r_interp(i,j)+box_len(i)
+                                                    else if(r_interp(i,j) > box_bd(2*i))  then
+                                                        vreal(i) = r_interp(i,j)-box_len(i)
+                                                    else 
+                                                        vreal(i) = r_interp(i,j)
+                                                    end if
+                                                else
+                                                    vreal(i) = r_interp(i,j)
+                                                end if
+                                            end do
+                                            v = matmul(transform_inverse, matmul(orient_inv, vreal))
+                                            do i = 1, 3
+                                                vlat(i) = nint(v(i) - bd_in_lat(2*i-1)+5)
+                                            end do
+
+                                            if(lat_points(vlat(1), vlat(2), vlat(3))) then  
+                                                lat_points(vlat(1), vlat(2), vlat(3)) = .false.
+                                            else
+                                                print *, "Lat points should be true not false"
+                                                stop 2
+                                            end if
+                                        end do 
+                                    end do
+                                end do                               
+                            else
+                                !Now set all the lattice points contained within an element to false
+                                do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))
+                                    do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:))
+                                        do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:))
+                                            lat_points(m,n,o) = .false.
+                                        end do 
+                                    end do
                                 end do
-                            end do
+                            end if
 
                         !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
@@ -515,8 +604,8 @@ module mode_create
                                 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), &
+                            lat_points_ele(1:(bd_ele_lat(2)-bd_ele_lat(1)+1),1:(bd_ele_lat(4)-bd_ele_lat(3)+1),&
+                            1:(bd_ele_lat(6)-bd_ele_lat(5))+1)= 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
@@ -537,7 +626,7 @@ module mode_create
                                     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
+                                                do inod = 1,nn
                                                     vlat = efill_ele(:,inod) + (/ i, j, k /)
                                                     if (.not.lat_points_ele(vlat(1),vlat(2),vlat(3))) cycle xe
                                                      do o = 1,3

src/mode_da.f90

@@ -0,0 +1,312 @@
+module mode_da
+    !This mode is used to calculate the dislocation analysis algorithm for both CG and atomistic regions
+
+    use parameters
+    use io
+    use elements
+    use neighbors
+    
+    implicit none 
+
+    integer :: num_d_points
+    integer, allocatable :: is_slipped(:,:,:) 
+    real(kind=dp), allocatable :: disl_line(:,:), disl_data(:)
+
+    private :: write_xyz
+    public 
+    contains
+
+    subroutine da(arg_pos)
+        !This is the main calling subroutine for the dislocation analysis code
+        integer, intent(out) :: arg_pos
+        real(kind=dp), dimension(6) :: temp_box_bd
+
+        integer :: ppos
+        character(len=100) :: outfile
+
+        !Now call the dislocation analysis code for the coarse-grained elements
+        call parse_command(arg_pos)
+
+        call read_in(1, (/0.0_dp,0.0_dp,0.0_dp/), temp_box_bd)
+
+        !Now allocate the necessary variables
+        if(allocated(is_slipped)) deallocate(is_slipped)
+        allocate(is_slipped(max_basisnum, max_ng_node, ele_num))
+
+        if(allocated(disl_line)) deallocate(disl_line)
+        num_d_points=0
+        allocate(disl_line(3,1024))
+        allocate(disl_data(1024))
+
+        call cga_da
+
+        !Now create the output file num and write out to xyz format
+        ppos = scan(trim(infiles(1)),".", BACK= .true.)
+        if ( ppos > 0 ) then 
+            outfile = 'da_'//infiles(1)(1:ppos)//'vtk'
+        else 
+            outfile = 'da_'//infiles(1)//'.vtk'
+        end if
+        call write_da_vtk(outfile)        
+
+        ppos = scan(trim(infiles(1)),".", BACK= .true.)
+        if ( ppos > 0 ) then 
+            outfile = 'da_'//infiles(1)(1:ppos)//'xyz'
+        else 
+            outfile = 'da_'//infiles(1)//'.xyz'
+        end if
+        call write_da_xyz(outfile)
+        
+        return
+    end subroutine da
+
+
+    subroutine parse_command(arg_pos)
+        !This subroutine parses the arguments for mode command
+        integer, intent(out)  :: arg_pos
+        integer :: i, arglen
+        character(len = 100) :: textholder
+        logical :: file_exists
+
+        !Read the input file 
+        call get_command_argument(2,textholder, arglen)
+        if (arglen == 0) stop "Missing file for dislocation analysis"
+        call get_in_file(textholder)
+
+        arg_pos = 3
+        return
+    end subroutine parse_command
+
+    subroutine cga_da
+
+        integer :: i, j, k, l, ibasis, nei, close_neighbors(2,4,6), far_neighbor(4,6), minindex, minindices(2), linevec(3)
+
+        integer :: face_types(ele_num*6), face_ele(6*ele_num), diff_pair(2,6)
+        real(kind = dp) :: face_centroids(3, ele_num*6), r(3), rc, vnode(3, max_basisnum, 4), vnorm, &
+                           max_node_dist, ndiff, rmax(3), rmin(3), rnorm, v1(3), v2(3), v1norm, v2norm, theta1, theta2, &
+                            diff_mag(6), diffvec(3)
+        real(kind=dp), allocatable :: disl_line_temp(:,:), disl_data_temp(:)
+        logical :: is_edge
+
+        !Initialize variables
+        l = 0
+        max_node_dist = 0
+
+        !Now save the close and far neighbors of the nodes. This is done to attempt to figure out the character of the dislocation
+        !by comparing how the burgers vector is distributed over the face.
+        do j = 1, 6
+            close_neighbors(1, 1, j) = cubic_faces(4, j)
+            close_neighbors(2, 1, j) = cubic_faces(2, j)
+            far_neighbor(1,j) = cubic_faces(3,j)
+
+            close_neighbors(1, 2, j) = cubic_faces(1, j)
+            close_neighbors(2, 2, j) = cubic_faces(3, j)
+            far_neighbor(2,j) = cubic_faces(4,j)
+
+            close_neighbors(1, 3, j) = cubic_faces(2, j)
+            close_neighbors(2, 3, j) = cubic_faces(4, j)
+            far_neighbor(3,j) = cubic_faces(1,j)
+
+            close_neighbors(1, 4, j) = cubic_faces(3, j)
+            close_neighbors(2, 4, j) = cubic_faces(1, j)
+            far_neighbor(4,j) = cubic_faces(2,j)
+        end do
+
+        !First calculate all of the face centroids
+        do i = 1, ele_num
+            do j = 1, 6
+                r(:) = 0.0_dp
+                rmax= -Huge(1.0)
+                rmin= Huge(1.0)
+                do k = 1, 4
+                    do ibasis = 1, basisnum(lat_ele(i))
+                        r = r + r_node(:, ibasis, cubic_faces(k,j), i) 
+                        rmax= max(rmax, r_node(:, ibasis, cubic_faces(k,j), i))
+                        rmin= min(rmin, r_node(:, ibasis, cubic_faces(k,j), i))
+                    end do
+                end do  
+                
+                rnorm=norm2(rmax-rmin)
+                max_node_dist = max(max_node_dist, rnorm)
+                r = r/(basisnum(lat_ele(i))*4)
+
+                !add the face centroids, the type, and map the elements faces to the face arrays
+                l = l + 1
+                face_centroids(:, l) = r 
+                face_types(l) = j
+                face_ele(l) = i
+
+
+            end do
+        end do
+
+        !Now set the cut off distance as half the distance from opposite nodes in the largest element
+        rc = 0.5*max_node_dist
+        call calc_NN(l, face_centroids(:,1:l), rc)
+
+        !Now loop overall the faces and make sure the nearest neighbor is the opposite type. If it isn't than we dscard it 
+        is_slipped = 0
+        do i = 1, l
+            nei = nn(i)             
+            
+            !Skip if it's 0 or the closest face belongs to the same element
+            if ((nei == 0).or.(face_ele(i) == face_ele(nei))) cycle
+            
+
+            !Check the face types, the way that the faces are laid out in the cubic_faces array face 1's opposite is 6 and 
+            ! face 2's opposite is 5 and etc
+            vnode = 0
+            if(face_types(i) == (7-face_types(nei))) then 
+                vnorm = 0
+                do j = 1, 4
+                    do ibasis = 1, basisnum(lat_ele(face_ele(i)))
+                        !Compute the vectors between all nodes at the face.
+                        vnode(:,ibasis,j) = r_node(:,ibasis, cubic_faces(j,face_types(i)), face_ele(i)) - &
+                                       r_node(:,ibasis, cubic_faces(j,face_types(nei)), face_ele(nei))
+                    end do
+                end do
+
+                !Now calculate the maximum difference between nodes
+                vnorm = 0
+                do j = 1, 4
+                    do k=j,4
+                        v1=vnode(:,1,j) - vnode(:,1,k)
+                        v1norm=norm2(v1)
+                        if (vnorm < v1norm) then
+                            vnorm = v1norm
+                            diffvec= v1
+                        end if 
+                    end do
+                end do
+                !Now calculate the difference between the largest norm and the smallest norm, if it's larger than 0.5 then mark it
+                !as slipped. This value probably can be converted to a variable value that depends on the current lattice parameter
+                !I think 0.5 works ok though.
+                if (vnorm>0.5_dp) then 
+                    print *, "Element number ", face_ele(i), " is dislocated along face ", face_types(i), &
+                             " with neighbor ", face_ele(nei), " with max displacement of ", vnorm
+                    l=0
+                    do j = 1, 4
+                        is_slipped(:, cubic_faces(j,face_types(i)), face_ele(i)) = 1
+
+                        !This portion of the code is used to determine what the character of the dislocation most likely is
+                        !effectively how this works is we look for the two nodes with the most similar burgers vector.
+                        !The vector between these two nodes is close to the line direction and as a result we can probably estimate
+                        !if it's close to edge, screw, or if it's pretty fairly mixed.
+                        do k = j+1, 4
+                            l=l+1
+                            diff_mag(l) = norm2(vnode(:, 1, j)-vnode(:,1,k))
+                            diff_pair(1,l)=cubic_faces(j, face_types(i))
+                            diff_pair(2,l)=cubic_faces(k, face_types(i))
+                        end do
+                    end do
+                    minindex = minloc(diff_mag,1)
+                    !Now figure out if the min difference between nodes is associated with a 112 direction    
+                    is_edge = .false.
+                    do j = 1,6
+                        if((diff_pair(1,minindex) == oneonetwopairs(1,j)).and.(diff_pair(2,minindex)==oneonetwopairs(2,j))) then
+                            is_edge=.true.
+                            exit
+                        else if((diff_pair(2,minindex)==oneonetwopairs(1,j)).and.(diff_pair(1,minindex)==oneonetwopairs(2,j)))then
+                            is_edge=.true.
+                            exit
+                        end if
+                    end do
+
+                    if(is_edge) then 
+                        print *, 'Dislocation has primarily edge character'
+                    else
+                        print *, 'Dislocation has primarily screw character'
+                    end if
+
+                    if( i < nei) then 
+                        num_d_points=num_d_points+2
+                        if(num_d_points > size(disl_line,2)) then 
+                            allocate(disl_line_temp(3,size(disl_line,2)+1024))
+                            disl_line_temp=0.0_dp
+                            disl_line_temp(:, 1:size(disl_line,2))=disl_line
+                            call move_alloc(disl_line_temp, disl_line)
+                        end if
+                        if(num_d_points/2 > size(disl_data)) then 
+                            allocate(disl_data_temp(size(disl_data)+1024))
+                            disl_data_temp=0
+                            disl_data_temp(1:size(disl_data)) = disl_data
+                            call move_alloc(disl_data_temp, disl_data)
+                        end if
+                        print *, num_d_points/2, vnorm
+                        disl_data(num_d_points/2)=vnorm
+                        disl_line(:,num_d_points-1) =r_node(:, 1, diff_pair(1,minindex), face_ele(i))
+                        disl_line(:,num_d_points) =r_node(:, 1, diff_pair(2,minindex), face_ele(i))
+                    end if
+                end if
+            end if
+        end do
+
+    end subroutine cga_da
+
+    subroutine write_da_vtk(outfile)
+        character (len=*), intent(in) :: outfile
+
+        integer :: i 
+        open(unit=11, file=outfile, action='write', status='replace', position='rewind')
+
+        write(11, '(a)') '# vtk DataFile Version 4.0.1'
+        write(11, '(a)') 'Dislocation line analyzed via cacmb dislocation analysis'
+        write(11, '(a)') 'ASCII'
+        write(11, '(a)') 'DATASET UNSTRUCTURED_GRID'
+        write(11, *) 'POINTS', num_d_points, 'float'
+       
+        !Write the points to the file
+        do i = 1, num_d_points
+            write(11, '(3f23.15)') disl_line(:, i)
+        end do
+
+        write(11,*) 'CELLS', num_d_points/2, num_d_points/2*3
+        do i=1, num_d_points, 2
+            write(11,'(3i16)') 2, i-1, i
+        end do
+        
+        write(11,*) 'CELL_TYPES', num_d_points/2
+        do i=1, num_d_points, 2
+            write(11,'(i16)') 3
+        end do
+
+        write(11,*) 'CELL_DATA', num_d_points/2
+        write(11, '(a)') 'SCALARS burgers_vec, float'
+        write(11, '(a)') 'LOOKUP_TABLE default'
+        do i=1, num_d_points, 2
+            write(11, '(f23.15)') disl_data((i+1)/2)
+        end do
+
+        close(11)
+
+    end subroutine write_da_vtk
+
+    subroutine write_da_xyz(outfile)
+        !This subroutine write the element positions to a .xyz file and marks whether they are slipped or not 
+        character(len=*), intent(in) :: outfile
+        integer :: i, ibasis, inod, outn
+        
+        open(unit=11, file=outfile, action='write', status='replace', position='rewind')
+
+        !Write number of node_atoms
+        write(11, '(i16)') node_atoms
+        !Write comment
+        write(11, '(a)') "is_slipped x y z"
+        !Write nodal positions
+        outn = 0
+        do i = 1, ele_num
+            do inod = 1, ng_node(lat_ele(i))
+                do ibasis = 1, basisnum(lat_ele(i))
+                        write(11, '(1i16, 3f23.15)') is_slipped(ibasis,inod,i), 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
+        
+        close(11)
+    end subroutine write_da_xyz
+end module mode_da

src/mode_merge.f90

@@ -6,18 +6,23 @@ module mode_merge
     use io
     use subroutines
     use elements
+    use neighbors
+
+    implicit none
 
     character(len=4) :: dim
     integer :: in_num, new_starts(2)
-    real(kind=dp) :: shift_vec(3)
-    logical :: shift_flag
+    real(kind=dp) :: shift_vec(3), replace_vec(3) 
+    character(len=100) :: replace_str(3)
+    logical :: shift_flag, replace_flag
+    real(kind=dp), private, save :: rc_off 
 
     public 
     contains
     subroutine merge(arg_pos)
         
         integer, intent(out) :: arg_pos 
-        integer :: i
+        integer :: i, j 
         real(kind=dp) :: displace(3), temp_box_bd(6)
 
         print *, '-----------------------Mode Merge---------------------------'
@@ -55,7 +60,27 @@ module mode_merge
                 call read_in(i, displace, temp_box_bd)
             end if 
 
+
             if(shift_flag) call shift(new_starts, i)
+
+            if(replace_flag.and.(i>1)) then 
+                !Parse the replace vector
+                do j = 1, 3
+                    call parse_pos(j, replace_str(j), replace_vec(j))
+                end do
+                call replace(new_starts, temp_box_bd)
+            end if
+
+        end do
+
+
+        !Now reset tags
+        do i = 1, atom_num
+            tag_atom(i) = i
+        end do
+        
+        do i = 1, ele_num
+            tag_ele(i) = i
         end do
 
         return
@@ -113,6 +138,17 @@ module mode_merge
                     if (arglen==0) stop "Missing vector component for shift command"
                     read(textholder, *) shift_vec(i)
                 end do
+            case('replace')
+                replace_flag = .true.
+                do i = 1,3
+                    arg_pos = arg_pos + 1
+                    call get_command_argument(arg_pos, replace_str(i), arglen)
+                    if (arglen==0) stop "Missing vector component for shift command"
+                end do
+                arg_pos = arg_pos+1
+                call get_command_argument(arg_pos, textholder, arglen)
+                read(textholder,*) rc_off
+
             case default
                 !If it isn't an available option to mode merge then we just exit
                 exit 
@@ -152,7 +188,7 @@ module mode_merge
         end do
 
         !If we don't include the wrap command then we have to increase the size of the box
-        if(.not.(wrap_flag)) then 
+        if(.not.(wrap_flag).and..not.(trim(adjustl(dim))=='none')) then 
             do i = 1,3
                 if (alldims(i:i) /= trim(adjustl(dim))) then 
                     if (current_shift(i) < -lim_zero) then 
@@ -167,4 +203,147 @@ module mode_merge
         end if
 
     end subroutine shift
+
+    subroutine replace(array_start, rbox_bd)
+        integer, intent(in) :: array_start(2)
+        real(kind = dp), intent(in) :: rbox_bd(6)
+
+        integer ::  ibasis, inod, del_num, del_index(atom_num), nump_ele, interp_start
+        integer ::  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, vlat(3), 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 :: in_bd, lat_points(max_esize, max_esize, max_esize)
+        real(kind=dp) :: del_bd(6)
+        integer :: i, c(3), ci, cj, ck, num_nei, nei, delete_num
+
+        !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(:,:,:,:)
+
+        !First apply the replace vec to all new nodes and elements
+        do i = array_start(1), atom_num
+            r_atom(:,i) = r_atom(:, i) + replace_vec
+        end do  
+
+        do i = array_start(2), ele_num
+            do inod = 1, ng_node(lat_ele(i))
+                do ibasis=1, basisnum(lat_ele(i))
+                    r_node(:, ibasis,inod, i) = r_node(:, ibasis,inod, i) + replace_vec 
+                end do
+            end do
+
+        end do
+
+        !Calculate new boundary
+        do i = 1, 6
+            del_bd(i) = rbox_bd(i) + replace_vec((i-1)/2 + 1)
+        end do
+        
+        del_num = 0
+        del_index=0
+        interp_start = atom_num +1
+        !Now loop over all old elements, 
+        do ie = 1, array_start(2)-1
+            !If any element points are within the boundary then we run the refine code
+            if(ele_in_bounds(del_bd, type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie))) then 
+                added_points=0
+                del_num = del_num + 1
+                del_index(del_num) = ie
+
+                !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 add the leftover lattice points as atoms, only if they aren't within the new boundaries
+                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))
+                                    added_points=added_points + 1
+                                    call add_atom(0, basis_type(ibasis,lat_ele(ie)), ratom(:,ibasis))
+                                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 if
+        end do
+
+        !Once all atoms are added we delete all of the elements
+        call delete_elements(del_num, del_index)
+
+        !Now delete overlapping atoms
+        allocate(which_cell(3,atom_num))
+
+        !First pass the atom list and atom num to the algorithm which builds the cell list
+        print *, rc_off
+        call build_cell_list(atom_num, r_atom, 4*rc_off, cell_num, num_in_cell, cell_list, which_cell)
+
+        !Now loop over every atom and figure out if it has neighbors within the rc_off
+        del_num = 0
+        atom_loop: do i = 1, atom_num
+            
+            !c is the position of the cell that the atom belongs to 
+            c = which_cell(:,i)
+
+            !Check to make sure it hasn't already been deleted
+            if(all(c /= 0)) then 
+                !Now loop over all neighboring cells
+                do ci = -1, 1, 1
+                    do cj = -1, 1, 1
+                        do ck = -1, 1, 1  
+
+                            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).and.(nei/= 0)) then 
+
+                                    !Now check to see if it is in the cutoff radius, if it is add it to the delete code
+                                    if (norm2(r_atom(:,nei)-r_atom(:,i)) < rc_off) then 
+                                        
+                                        del_num = del_num + 1
+
+                                        !Make sure to delete the older value
+                                        if( (i < array_start(1)).or.(i > interp_start)) then 
+                                            del_index(del_num) = i
+                                            which_cell(:,i) = 0
+                                            cycle atom_loop
+                                        else
+                                            del_index(del_num) = nei
+                                            which_cell(:,nei) = 0 
+                                            cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci) = 0
+                                        end if
+                                    end if
+                                end if
+                            end do 
+                        end do
+                    end do 
+                end do 
+            end if
+
+        end do atom_loop
+
+        print *, "Replace command deletes ", del_num, " atoms"
+        !Now delete all the atoms
+        call delete_atoms(del_num, del_index(1:del_num))
+
+
+        return
+    end subroutine replace
+
 end module mode_merge

src/mode_metric.f90

@@ -10,7 +10,7 @@ module mode_metric
     
     integer :: nfiles
     character(len=100) :: metric_type
-    real(kind=dp) :: rc_off
+    real(kind=dp) :: rc_off, old_box_len(3)
 
     !Save reference positions
     integer :: np, npreal, nmet
@@ -50,23 +50,31 @@ module mode_metric
         end select
 
         !Now set the reference positions
+        print *, "Combining 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)
+        print *, "Getting initial neighbor list"
+        call calc_neighbor(rc_off, r_zero, np)
 
+        !Save box lengths for applying periodic boundaries 
+        do i=1,3
+            old_box_len(i) = box_bd(2*i) - box_bd(2*i-1)
+        end do
         !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
+            print *, "Calculating metric for file number ", i
             call calc_metric
             !Now create the output file num and write out to xyz format
             ppos = scan(trim(infiles(i)),".", BACK= .true.)
@@ -124,8 +132,11 @@ module mode_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)
+                         U(3,3), R(3,3), Rskew(3,3), oldrij(3), box_len(3)
 
+        do i = 1, 3
+            box_len(i) = box_bd(2*i) - box_bd(2*i-1)
+        end do
         !Loop over all points
         do ip = 1, np
             eta(:,:) = 0.0_dp
@@ -134,8 +145,17 @@ module mode_metric
             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)
+                rij = r_curr(:,nei) - r_curr(:,ip)  +periodvec(:,jp,ip)*box_len 
+                if (norm2(rij) > 5*rc_off) then 
+                    do j = 1, 3
+                        if (r_curr(j,nei) > r_curr(j,ip)) then 
+                            rij(j) = rij(j) - box_len(j)
+                        else if(r_curr(j,nei) < r_curr(j,ip)) then 
+                            rij(j) = rij(j) + box_len(j)
+                        end if
+                    end do
+                end if
+                oldrij = r_zero(:,nei) - r_zero(:,ip) + periodvec(:,jp,ip)*old_box_len 
 
                 !Calculate  eta and omega
                 do i = 1,3
@@ -190,9 +210,10 @@ module mode_metric
         integer :: i, inod, ibasis
 
         npoints=0
-
+        rout = -huge(1.0_dp)
         if(atom_num > 0) then 
             do i = 1, atom_num
+                call apply_periodic(r_atom(:,i))
                 rout(:,tag_atom(i))  = r_atom(:,i)
                 npoints = npoints + 1
             end do
@@ -202,8 +223,11 @@ module mode_metric
             do i = 1, ele_num
                 do inod = 1, ng_node(lat_ele(i))
                     do ibasis = 1, basisnum(lat_ele(i))
+                        call apply_periodic(r_node(:,ibasis,inod,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
@@ -223,14 +247,14 @@ module mode_metric
         
         select case(metric_type)
         case('def_grad')
-            write(11,*) "type element x y z F11 F12 F13 F21 F22 F23 F31 F32 F33"
+            write(11,*) "id 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)"
+            write(11,*) "id 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))
+                write(11,*) tag_atom(i), type_atom(i), 0, r_atom(:,i), met(:,tag_atom(i))
             end do
         end if
 
@@ -238,7 +262,8 @@ module mode_metric
             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), &
+                        write(11,*) atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis, &
+                                    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

src/neighbors.f90

@@ -4,10 +4,12 @@ module neighbors
     use elements
     use subroutines
     use functions
+    use box
 
-    integer, allocatable :: nei_list(:,:), nei_num(:)
+    integer, allocatable :: nei_list(:,:), nei_num(:), nn(:), periodvec(:,:,:)
     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)
@@ -47,22 +49,25 @@ module neighbors
         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
+        num_in_cell = 0
+        
         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
+            if(r_list(1,i) < box_bd(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
+                c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off)) + 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(1:cell_lim, :, :, :) = cell_list(1:cell_lim, :, :, :)
                 resize_cell_list(cell_lim+1:, :, :, :) = 0
                 call move_alloc(resize_cell_list, cell_list)
+                cell_lim = cell_lim + 10
             end if
 
             cell_list(num_in_cell(c(1),c(2),c(3)),c(1),c(2),c(3)) = i
@@ -77,44 +82,70 @@ module neighbors
         
         real(kind=dp), intent(in) :: rc_off
         integer, intent(in) :: n
-        real(kind=dp), dimension(3,n) :: r_list
+        real(kind=dp), dimension(:, :), intent(in) :: r_list
 
-        integer :: i, c(3), ci, cj, ck, num_nei, nei
+        integer :: i, j, c(3),cn(3), ci, cj, ck, num_nei, nei, v(3), period_dir(3)
         !Variables for cell list code
         integer, dimension(3) ::cell_num
         integer, allocatable :: num_in_cell(:,:,:), cell_list(:,:,:,:)
         integer :: which_cell(3,n)
+        real(kind=dp) :: r(3), box_len(3)
+        logical :: period_bd(3)
         
-        !First reallocate the neighbor list codes
+        !First reallocate the neighbor list variables
         if (allocated(nei_list)) then
-            deallocate(nei_list,nei_num)
+            deallocate(nei_list,nei_num, periodvec)
         end if
 
-        allocate(nei_list(100,n),nei_num(n))
+        allocate(nei_list(100,n),nei_num(n), periodvec(3,100,n))
+        nei_list=0
+        periodvec=0
+        nei_num=0
 
         !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)
 
+        do i=1, 3
+            if (box_bc(i:i) == 'p') then 
+                period_bd(i) = .true.
+            else
+                period_bd(i)=.false.
+            end if
+            box_len(i) = box_bd(2*i) - box_bd(2*i-1)
+        end do
         !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
+            if(r_list(1,i) < box_bd(1)) cycle
 
-            !c is the position of the cell that the point
+            !c is the position of the cell that the point belongs to 
             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
+                        !First try to apply periodic boundaries 
+                        v=(/ ck, cj, ci /)
+                        cn=0
+                        period_dir=0
+                        do j=1, 3
+                            if ((c(j) + v(j) == 0).and.period_bd(j)) then 
+                                cn(j) = cell_num(j)
+                                period_dir(j) = 1
+                            else if ((c(j) + v(j) > cell_num(j)).and.period_bd(j)) then 
+                                cn(j) = 1
+                                period_dir(j) = -1
+
+                            else if ((c(j)+v(j) >= 1) .and. (c(j)+v(j) <= cell_num(j))) then
+                                cn(j) = c(j) + v(j)
+                            end if
+                        end do
                         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)
+                        do num_nei = 1, num_in_cell(cn(1),cn(2),cn(3))
+                            nei = cell_list(num_nei,cn(1),cn(2),cn(3))
 
                             !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
@@ -122,10 +153,13 @@ module neighbors
 
                                 !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 
+                        
+                                r = r_list(:,nei) + period_dir*box_len
+                                if (norm2(r-r_list(:,i)) < rc_off) then 
                                     
                                     nei_num(i) = nei_num(i) + 1
                                     nei_list(nei_num(i), i) = nei
+                                    periodvec(:,nei_num(i),i) = period_dir
 
                                 end if
                             end if
@@ -139,4 +173,66 @@ module neighbors
         return
     end subroutine calc_neighbor
 
+    subroutine calc_NN(n,points, rc_off)
+
+        integer, intent(in) :: n
+        real(kind=dp), intent(in) :: points(3,n)
+        real(kind=dp), intent(in) :: rc_off
+
+        integer :: i, c(3), ci, cj, ck, nei
+        !cell arrays
+        integer, dimension(3) ::cell_num
+        integer, allocatable :: num_in_cell(:,:,:), cell_list(:,:,:,:)
+        integer :: which_cell(3,n)
+        real(kind = dp) :: rmin
+
+        !First reallocate the neighbor list codes
+        if (allocated(nn)) then
+            deallocate(nn)
+        end if
+
+        allocate(nn(n))
+        nn=0
+
+        !First build the cell lists
+        call build_cell_list(n, points, rc_off, cell_num, num_in_cell, cell_list, which_cell)
+        
+        pointloop: do i = 1, n
+            !First check to see if the point is a filler point, if so then skip it
+            if(points(1,i) < -Huge(-1.0_dp)+1) cycle
+
+            !Get the positon of the cell 
+            c = which_cell(:,i)
+
+            !Initialize the min vec
+            rmin=rc_off
+
+            !loop over all neighboring cells
+            do ci = -1, 1, 1
+                do cj = -1, 1, 1
+                    do ck = -1, 1, 1
+                        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 
+                            if((nei /= i)) then 
+                                !If it's the minimum position than we add it to the nearest neighbor list and updat e the min vec
+                                
+                                if (norm2(points(:,nei)-points(:,i)) < rmin) then 
+                                    rmin = norm2(points(:, nei) - points(:,i))
+                                    nn(i)=(nei)
+                                end if
+                            end if
+                        end do 
+                    end do
+                end do
+            end do
+        end do pointloop
+
+        return
+    end subroutine calc_NN
+
 end module neighbors

src/opt_bubble.f90

@@ -0,0 +1,146 @@
+module opt_bubble
+    !This module contains the bubble option which can be used to create voids with specific pressures of certain atoms
+
+    use atoms
+    use parameters
+    use elements 
+    use box
+    use opt_group
+
+    implicit none
+
+    real(kind=dp), private :: br, brat, c(3)
+    character(len=2), private :: species
+
+    public 
+    contains
+    subroutine bubble(arg_pos)
+        integer, intent(inout) :: arg_pos
+
+        integer :: new_type, n, j, i, atom_num_pre
+        real(kind = dp)  :: p(3), rand, factor, per, vol, mass
+
+        print *, '------------------------------------------------------------'
+        print *, 'Option Bubble'
+        print *, '------------------------------------------------------------'
+        !First we parse the bubble command
+        call parse_bubble(arg_pos)
+
+        !Now we use the existing group code to delete a sphere which represents the bubble
+        centroid=c
+        radius = br
+        type='all'
+        gshape='sphere'
+        group_nodes = .true.
+        group_atom_types=0
+        call get_group
+        call refine_group
+        call get_group
+        atom_num_pre= atom_num
+        call delete_group
+
+        !Now we create a new atom type with the desired species
+        call atommass(species, mass)
+        call add_atom_type(mass, new_type)
+
+        !Now we calculate the number of atoms we need to add for the desired pressure
+        !print *, "Creating a bubble with pressure", bp, " at temperature ", bt, " with radius ", br
+        !
+        !factor = 1.0e24/6.02214e23
+        !if (bp <= 10.0) then
+        !    per=factor*(3.29674113+4.51777872*bp**(-0.80473167))
+        !else if (bp .le. 20.0) then
+        !    per=factor*(4.73419689-0.072919483*bp)
+        !else
+        !    per=factor*(4.73419689-0.072919483*bp)
+        !    print *, 'warning: pressure is too high'
+        !    print *, 'equation of state is only valid for < 20 GPa'
+        !endif 
+        !vol = 4.0*pi/3.0*br**3.0
+
+        n = brat*(atom_num_pre-atom_num)
+        print *, "Adding ", n, " atoms of species ", species
+
+        !Now add n atoms randomly within the sphere 
+        do i = 1, n
+            do while(.true.)
+                do j = 1, 3
+                    call random_number(rand)
+                    p(j) = rand*(2*br) + c(j)-br
+                end do
+                if (norm2(p-c) < br) exit
+            end do
+            
+            call add_atom(0, new_type, p)
+        end do
+
+    end subroutine bubble 
+
+    subroutine parse_bubble(arg_pos) 
+        integer, intent(inout) :: arg_pos
+
+        integer :: i, arglen
+        real(kind=dp) :: mass
+        character(len=100) :: tmptxt
+
+        
+        !First read in the bubble centroid
+        do i = 1, 3
+            arg_pos = arg_pos + 1 
+            call get_command_argument(arg_pos, tmptxt, arglen)
+            print *,  trim(adjustl(tmptxt))
+            call parse_pos(i, tmptxt, c(i))
+        end do
+
+        !Now the bubble radius
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, tmptxt, arglen)
+        print *,  trim(adjustl(tmptxt))
+        if(arglen == 0) stop "Missing bubble radius"
+        read(tmptxt, *) br
+
+        !Now bubble ratio
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, tmptxt, arglen)
+        print *,  trim(adjustl(tmptxt))
+        if(arglen == 0) stop "Missing bubble ratio"
+        read(tmptxt, *) brat
+        
+        !Now the bubble species
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, species, arglen)
+        print *,  species
+        if(arglen == 0) stop "Missing bubble species"
+
+
+        !OPtional arguments
+        do while(.true.)
+            if(arg_pos > command_argument_count()) exit
+            arg_pos=arg_pos+1
+
+            call get_command_argument(arg_pos, tmptxt)
+            tmptxt=adjustl(tmptxt)
+            select case(trim(tmptxt))
+
+            case('excludetypes')
+                arg_pos=arg_pos + 1
+                call get_command_argument(arg_pos, tmptxt, arglen)
+                if(arglen==0) stop "Missing number of atoms to exclude"
+                read(tmptxt, *) exclude_num
+                do i=1,exclude_num
+                    arg_pos = arg_pos + 1
+                    call get_command_argument(arg_pos, tmptxt, arglen)
+                    if(arglen==0) stop "Missing exclude atom types"
+                    call atommass(tmptxt, mass)
+                    call add_atom_type(mass, exclude_types(i))
+                end do
+
+            case default
+                exit
+            end select
+        end do
+        
+        return
+    end subroutine parse_bubble
+end module opt_bubble
+

src/opt_deform.f90

@@ -10,6 +10,7 @@ module opt_deform
 
     real(kind=dp), save :: applied_strain
     integer, save :: sdim
+    logical :: percent 
 
     public 
     contains
@@ -21,9 +22,11 @@ module opt_deform
 
         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)
+        real(kind=dp) :: frac_atom(atom_num), frac_node(max_basisnum, max_ng_node, ele_num), blen
+
 
         !initialize some variables
+        percent=.false.
         dims(1) = 'x'
         dims(2) = 'y'
         dims(3) = 'z'
@@ -44,7 +47,13 @@ module opt_deform
         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
+        if(percent) then 
+            blen = box_bd(2*sdim) - box_bd(2*sdim-1) 
+            box_bd(2*sdim) = box_bd(2*sdim-1) + blen*applied_strain
+        else
+            box_bd(2*sdim) = box_bd(2*sdim) + applied_strain
+        end if
+
         print *, "New box bounds are ", box_bd(2*sdim-1:2*sdim)
 
         !Now reassign the positions
@@ -90,6 +99,25 @@ module opt_deform
         if (arg_len == 0) stop "Missing strain in deform command"
         read(textholder, *) applied_strain
 
+        !Now parse the additional options which may be present
+        do while(.true.)
+            if(arg_pos > command_argument_count()) exit
+            !Pull out the next argument which should either be a keyword or an option
+            arg_pos=arg_pos+1
+            call get_command_argument(arg_pos, textholder)
+            textholder=adjustl(textholder)
+
+            !Choose what to based on what the option string is 
+            select case(trim(textholder))
+            case('percent')
+                percent=.true. 
+            case default
+                arg_pos=arg_pos-1
+                !if it isn't an available option to opt_group then we just exit
+                exit 
+            end select
+        end do
+
         arg_pos = arg_pos + 1
 
     end subroutine parse_deform

src/opt_delete.f90

@@ -8,6 +8,7 @@ module opt_delete
     implicit none 
 
     real(kind=dp) :: rc_off
+    logical :: first
 
     public  
     contains
@@ -34,6 +35,7 @@ module opt_delete
         character(len=100) :: textholder
         arg_pos = arg_pos + 1
 
+        first = .true.
         call get_command_argument(arg_pos, textholder, arg_len)
         if(arg_len==0) stop "Missing argument to delete command"
 
@@ -49,6 +51,24 @@ module opt_delete
         end select 
 
         arg_pos = arg_pos + 1
+
+        do while(.true.)
+            if(arg_pos > command_argument_count()) exit
+            !Pull out the next argument which should either be a keyword or an option
+            arg_pos=arg_pos+1
+            call get_command_argument(arg_pos, textholder)
+            textholder=adjustl(textholder)
+
+            select case(trim(textholder))
+            case('first')
+                first=.true.
+            case('last')
+                first=.false.
+            case default
+                !if it isn't an available option to opt_group then we just exit
+                exit 
+            end select
+        end do
     end subroutine parse_delete
 
     subroutine delete_overlap
@@ -62,14 +82,16 @@ module opt_delete
         integer, dimension(3) :: cell_num
         integer, allocatable :: num_in_cell(:,:,:), which_cell(:,:)
         integer, allocatable :: cell_list(:,:,:,:)
+        logical :: deleted_atom(atom_num)
         
         allocate(which_cell(3,atom_num))
 
         !First pass the atom list and atom num to the algorithm which builds the cell list
-        call build_cell_list(atom_num, r_atom, rc_off, cell_num, num_in_cell, cell_list, which_cell)
+        call build_cell_list(atom_num, r_atom, 5*rc_off, cell_num, num_in_cell, cell_list, which_cell)
 
         !Now loop over every atom and figure out if it has neighbors within the rc_off
         delete_num = 0
+        deleted_atom=.false.
         atom_loop: do i = 1, atom_num
             
             !c is the position of the cell that the atom belongs to 
@@ -93,22 +115,39 @@ module opt_delete
 
                                 !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).and.(nei/= 0)) then 
+                                if((nei /= i).and.(.not.deleted_atom(nei))) then 
 
                                     !Now check to see if it is in the cutoff radius, if it is add it to the delete code
                                     if (norm2(r_atom(:,nei)-r_atom(:,i)) < rc_off) then 
                                         
                                         delete_num = delete_num + 1
-                                        for_delete(delete_num) = max(i,nei)
+                                        if(first) then
+                                            for_delete(delete_num) = min(i,nei)
+                                            deleted_atom(min(i,nei)) = .true.
+                                        else
+                                            for_delete(delete_num) = max(i,nei)
+                                            deleted_atom(max(i,nei)) = .true.
+                                        end if
 
                                         !Now zero out the larger index
-                                        if(i > nei) then 
-                                            which_cell(:,i) = 0
-                                            cycle atom_loop
+                                        if(first) then 
+                                            if(i < nei) then 
+                                                which_cell(:,i) = 0
+                                                cycle atom_loop
+                                            else
+                                                which_cell(:,nei) = 0 
+                                                cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci) = 0
+                                            end if
                                         else
-                                            which_cell(:,nei) = 0 
-                                            cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci) = 0
+                                            if(i > nei) then 
+                                                which_cell(:,i) = 0
+                                                cycle atom_loop
+                                            else
+                                                which_cell(:,nei) = 0 
+                                                cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci) = 0
+                                            end if
                                         end if
+
                                     end if
                                 end if
                             end do 

src/opt_disl.f90

@@ -30,7 +30,7 @@ module opt_disl
         integer, intent(inout) :: arg_pos 
 
         print *, '--------------------Option Dislocation-----------------------'
-
+        
         select case(trim(adjustl(option)))
         case('-disl')
             call parse_disl(arg_pos)
@@ -113,7 +113,7 @@ module opt_disl
     subroutine disl
         !This subroutine creates the actual dislocation
 
-        integer :: i, sub_box, inod, ibasis, ix, iy, iz
+        integer :: i, 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)
 
@@ -240,7 +240,7 @@ module opt_disl
     subroutine dislgen
         !This subroutine creates the actual dislocation
 
-        integer :: i, sub_box, inod, ibasis
+        integer :: i, inod, ibasis
         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)
 
@@ -250,14 +250,6 @@ module opt_disl
         be = sin(char_angle*pi/180.0_dp)*b
         bs = cos(char_angle*pi/180.0_dp)*b
 
-        !Figure out which sub box you are in so you can access it's orientation, this code will not work
-        !with overlapping sub_boxes
-        do i = 1, sub_box_num
-            if(in_block_bd(centroid, sub_box_bd(:,i))) then 
-                sub_box = i
-                exit
-            end if
-        end do
 
         !Construct the slip system orientation matrix in an unrotated system 
         slipx = cross_product(slip_plane, line)
@@ -267,8 +259,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)
+        disp_transform = matmul(box_ori(:,:), ss_inv)
         call matrix_inverse(disp_transform,3,inv_transform)
 
         if(atom_num > 0) then 
@@ -381,18 +372,6 @@ module opt_disl
 
         arg_pos = arg_pos + 1
 
-        !Now check to make sure that the dimension selected is actually a 1 1 1 direction.
-        ! call in_sub_box(centroid, sbox)
-        ! if(.not.((abs(sub_box_ori(loop_normal,1,sbox)) == abs(sub_box_ori(loop_normal,2,sbox))).and. &
-        !          (abs(sub_box_ori(loop_normal,2,sbox)) == abs(sub_box_ori(loop_normal,3,sbox))).and. &
-        !          (abs(sub_box_ori(loop_normal,3,sbox)) == abs(sub_box_ori(loop_normal,1,sbox))))) then 
-
-        !     print *, "The selected dimension ", loop_normal, " for sub_box ", sbox, " is ", &
-        !              sub_box_ori(loop_normal,:,sbox), " which is not in the (111) family of planes"
-        !     STOP 3 
-        ! end if
-
-
     end subroutine
 
     !Code for the creation of dislocation loops is based on functions from atomsk
@@ -541,7 +520,7 @@ module opt_disl
         end do
         return
 
-    end subroutine
+    end subroutine disloop
 
     !********************************************************
     ! SOLIDANGLE
@@ -652,18 +631,6 @@ module opt_disl
 
         arg_pos=arg_pos+1
 
-        !Now check to make sure that the dimension selected is actually a 1 1 1 direction.
-        ! call in_sub_box(centroid, sbox)
-        ! if(.not.((abs(sub_box_ori(loop_normal,1,sbox)) == abs(sub_box_ori(loop_normal,2,sbox))).and. &
-        !          (abs(sub_box_ori(loop_normal,2,sbox)) == abs(sub_box_ori(loop_normal,3,sbox))).and. &
-        !          (abs(sub_box_ori(loop_normal,3,sbox)) == abs(sub_box_ori(loop_normal,1,sbox))))) then 
-
-        !     print *, "The selected dimension ", loop_normal, " for sub_box ", sbox, " is ", &
-        !              sub_box_ori(loop_normal,:,sbox), " which is not in the (111) family of planes"
-        !     STOP 3 
-        ! end if
-
-
     end subroutine parse_vacancydisloop
 
 

src/opt_orient.f90

@@ -10,7 +10,7 @@ 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(:,:,:)
+    real(kind=dp), dimension(3,3) :: orig_box_ori
     character(len=3) :: orig_box_bc 
 
     public 
@@ -22,7 +22,7 @@ module opt_orient
 
         integer :: i, j, k, ibasis, inod
         logical :: isortho, isrighthanded, matching
-        real(kind=dp) :: inv_sub_box_ori(3,3,sub_box_num)
+        real(kind=dp) :: inv_box_ori(3,3)
 
 
         !First parse the orient command
@@ -31,36 +31,25 @@ module opt_orient
         !Now rotate the basis. To do this we transform the basis to [100] [010] [001] and then 
         !transform to user specified basis.
 
-        !Find all inverse orientation matrices for all sub_boxes
-        do i = 1, sub_box_num
-            call matrix_inverse(sub_box_ori(:,:,i), 3, inv_sub_box_ori(:,:,i))
-        end do
+        call matrix_inverse(box_ori(:,:), 3, inv_box_ori(:,:))
 
         !Now transform all atoms
         do i = 1, atom_num
-            r_atom(:,i) = matmul(new_orient,matmul(inv_sub_box_ori(:,:,sbox_atom(i)),r_atom(:,i)))
+            r_atom(:,i) = matmul(new_orient,matmul(inv_box_ori(:,:),r_atom(:,i)))
         end do 
 
         !Now transform all elements
         do i = 1, ele_num
             do inod =1, ng_node(lat_ele(i))
                 do ibasis = 1, basisnum(lat_ele(i))
-                    r_node(:,ibasis,inod,i) = matmul(new_orient,matmul(inv_sub_box_ori(:,:,sbox_ele(i)),r_node(:,ibasis,inod,i)))
+                    r_node(:,ibasis,inod,i) = matmul(new_orient,matmul(inv_box_ori(:,:),r_node(:,ibasis,inod,i)))
                 end do
             end do  
         end do
 
-        !Now save the original sub_box_ori and overwrite them
-        if(allocated(orig_sub_box_ori)) deallocate(orig_sub_box_ori)
-
-        allocate(orig_sub_box_ori(3,3,sub_box_num))
-        orig_sub_box_ori = sub_box_ori
-
-        !Now overwrite the orientations
-        do i = 1, sub_box_num
-            sub_box_ori(:,:,i) = new_orient
-        end do 
-
+        !Now save the original box_ori 
+        orig_box_ori = box_ori
+        box_ori = new_orient
         !Save original box boundaries
         orig_box_bd = box_bd
 
@@ -70,18 +59,16 @@ module opt_orient
         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
+            do k = 1, 3
+                if(.not.is_equal(orig_box_ori(i,k), new_orient(i,k))) then 
+                    matching = .false.
+                    exit
+                end if
+            end do 
             if(.not.matching) box_bc(i:i)='s'
         end do
 
-            call def_new_box
+        call def_new_box
     end subroutine orient_opt
 
     subroutine parse_orient(arg_pos)
@@ -121,24 +108,24 @@ module opt_orient
         real(kind=dp) :: inv_ori(3,3)
 
         !Now rotate the basis. To do this we transform the basis to [100] [010] [001] and then 
-        !transform to the original sbox_ele
+        !transform to the original box_ori
 
         !Find the inverse for the new orientation matrix
         call matrix_inverse(new_orient, 3, inv_ori)
 
-        !Recover original sub_box_ori
-        sub_box_ori = orig_sub_box_ori 
+        !Recover original box_ori and box_bd
+        box_ori = orig_box_ori 
 
         !Now transform all atoms
         do i = 1, atom_num
-            r_atom(:,i) = matmul(sub_box_ori(:,:,sbox_atom(i)),matmul(inv_ori(:,:),r_atom(:,i)))
+            r_atom(:,i) = matmul(box_ori(:,:),matmul(inv_ori(:,:),r_atom(:,i)))
         end do 
 
         !Now transform all elements
         do i = 1, ele_num
             do inod =1, ng_node(lat_ele(i))
                 do ibasis = 1, basisnum(lat_ele(i))
-                    r_node(:,ibasis,inod,i) = matmul(sub_box_ori(:,:,sbox_ele(i)),matmul(inv_ori,r_node(:,ibasis,inod,i)))
+                    r_node(:,ibasis,inod,i) = matmul(box_ori(:,:),matmul(inv_ori,r_node(:,ibasis,inod,i)))
                 end do
             end do  
         end do
@@ -148,29 +135,4 @@ module opt_orient
         box_bc = orig_box_bc
     end subroutine unorient
 
-    subroutine sbox_ori(arg_pos)
-        integer, intent(inout) :: arg_pos
-
-        integer :: i, sbox_in, arg_len
-        real(kind = dp) :: new_orient(3,3)
-        character(len=100) :: textholder 
-        character(len=20) :: ori_string
-
-        arg_pos = arg_pos + 1
-        call get_command_argument(arg_pos, textholder,arg_len)
-        if (arg_len== 0) stop 'Missing sbox in sbox_ori command'
-        read(textholder,*) sbox_in
-
-        do i = 1, 3
-            arg_pos = arg_pos + 1
-            call get_command_argument(arg_pos, ori_string, arg_len)
-            if (arg_len == 0) print *, "Missing orientation vector in -orient option"
-            call parse_ori_vec(ori_string, new_orient(i,:))
-            new_orient(i,:) = new_orient(i,:) / norm2(new_orient(i,:))
-        end do
-
-        sub_box_ori(:,:,sbox_in) = new_orient
-        arg_pos = arg_pos + 1
-
-    end subroutine sbox_ori
 end module opt_orient

src/opt_redef_box.f90

@@ -3,6 +3,7 @@ module opt_redef_box
     use box
     use elements
     use subroutines
+
     implicit none
     
     character(len=3) :: redef_dim, new_bc
@@ -13,7 +14,8 @@ module opt_redef_box
     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)
+        integer :: i, j, inod, ibasis, iatom, delete_list(atom_num), delete_num, type_interp(max_basisnum*max_esize**3), &
+                   new_snum 
         real(kind=dp):: r_interp(3, max_basisnum*max_esize**3)
         logical :: node_out(8)
 
@@ -64,7 +66,7 @@ module opt_redef_box
                 !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))
+                        call add_atom(0,type_interp(iatom),  r_interp(:,iatom))
                     end if
                 end do 
             end if
@@ -128,4 +130,26 @@ module opt_redef_box
         arg_pos = arg_pos + 1
     end subroutine parse_redef_box
 
+    subroutine set_box_ori(arg_pos)
+
+        integer, intent(inout) :: arg_pos 
+
+        integer :: arg_len, i
+        real(kind=dp) :: new_orient(3,3)
+        character(len=100) :: textholder
+        character(len=20) :: ori_string
+
+        do i = 1, 3
+            arg_pos = arg_pos + 1
+            call get_command_argument(arg_pos, ori_string, arg_len)
+            if (arg_len == 0) print *, "Missing orientation vector in -orient option"
+            call parse_ori_vec(ori_string, new_orient(i,:))
+            new_orient(i,:) = new_orient(i,:) / norm2(new_orient(i,:))
+        end do
+
+        box_ori(:,:) = new_orient
+        arg_pos = arg_pos + 1
+
+    end subroutine set_box_ori
+
 end module opt_redef_box

src/opt_slip_plane.f90

@@ -17,11 +17,13 @@ module opt_slip_plane
         !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 :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3), ibasis,&
+                    starting_anum
                 
         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
+                         maxp, minp, vel_interp(3, max_basisnum*max_esize**3)
+        real(kind=dp), allocatable :: vel_tmp(:,:)
 
         integer :: type_interp(max_basisnum*max_esize**3)
         logical :: lat_points(max_esize,max_esize, max_esize)
@@ -58,11 +60,25 @@ module opt_slip_plane
                 
                 !If we aren't efilling then just refine the element
                 if(.not.efill) then 
+                    starting_anum=atom_num
                     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))
+                        call add_atom(0, type_interp(ia), r_interp(:,ia))
                     end do
+
+                    if(allocated(vel_atom)) then  
+                        call interpolate_vel(type_ele(ie), size_ele(ie), lat_ele(ie), vel_node(:,:,:,ie), vel_interp)
+                        if(size(vel_atom,2) < size(type_atom)) then 
+                            allocate(vel_tmp(3, size(type_atom)))
+                            vel_tmp=0.0_dp
+                            vel_tmp(:, 1:size(vel_atom,2)) = vel_atom
+                            call move_alloc(vel_tmp, vel_atom)
+                        end if
+                        do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3           
+                            vel_atom(:, starting_anum+ia) = vel_interp(:, ia)
+                        end do
+                    end if
                 !If we are efilling then the code is slightly more complex
                 else
                     !First populate the lat points array 
@@ -93,7 +109,7 @@ module opt_slip_plane
                                         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) 
+                                            call add_element(0, type_ele(ie), esize, lat_ele(ie),  rfill) 
                                             new_ele_num = new_ele_num + 1
                                         end if
                                     end do latloop
@@ -110,7 +126,7 @@ module opt_slip_plane
                                     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))
+                                        call add_atom(0, basis_type(ibasis,lat_ele(ie)), ratom(:,ibasis))
                                     end do
                                 end if
                             end do

src/parameters.f90

@@ -13,6 +13,7 @@ module parameters
     
     !Numeric constants
     real(kind=dp), parameter :: pi = 3.14159265358979323846_dp
+    real(kind=dp), parameter :: pvt2n =  1362.626470955822
     
     !Mode type that is being called 
     character(len=100) :: mode

src/subroutines.f90

@@ -2,6 +2,8 @@ module subroutines
     use parameters
     use functions
     use box
+    use str
+
     implicit none
 
     integer :: allostat, deallostat
@@ -149,29 +151,43 @@ module subroutines
 
     subroutine parse_ori_vec(ori_string, ori_vec)
         !This subroutine parses a string to vector in the format [ijk]
-        character(len=20), intent(in) :: ori_string
+        character(*), intent(in) :: ori_string
         real(kind=dp), dimension(3), intent(out) :: ori_vec
+        real(kind=dp) :: prefactor
+        integer :: start, fin
 
         integer :: i, ori_pos, stat
 
-        ori_pos=2
-        do i = 1,3
-            do while(ori_string(ori_pos:ori_pos)==' ')    
-                ori_pos=ori_pos+1
-            end do
-            if (ori_string(ori_pos:ori_pos) == '-') then 
-                ori_pos = ori_pos + 1
-                read(ori_string(ori_pos:ori_pos), *, iostat=stat) ori_vec(i)
-                if (stat>0) STOP "Error reading orientation value"
-                ori_vec(i) = -ori_vec(i)
-                ori_pos = ori_pos + 1
-            else
-                read(ori_string(ori_pos:ori_pos), *, iostat=stat) ori_vec(i)
-                if(stat>0) STOP "Error reading orientation value"
-                ori_pos=ori_pos + 1
+        do i = 1, len(ori_string)
+            if(ori_string(i:i) == "[") then 
+                start=i+1
+            else if(ori_string(i:i) == "]") then 
+                fin = i-1
             end if
         end do
-
+        
+        ori_pos=start
+
+        if (tok_count(ori_string(start:fin)) == 3) then 
+            read(ori_string(start:fin),*) ori_vec(:)
+        else
+            do i = 1,3
+                do while(ori_string(ori_pos:ori_pos)==' ')    
+                    ori_pos=ori_pos+1
+                end do
+                if (ori_string(ori_pos:ori_pos) == '-') then 
+                    ori_pos = ori_pos + 1
+                    read(ori_string(ori_pos:ori_pos), *, iostat=stat) ori_vec(i)
+                    if (stat>0) STOP "Error reading orientation value"
+                    ori_vec(i) = -ori_vec(i)
+                    ori_pos = ori_pos + 1
+                else
+                    read(ori_string(ori_pos:ori_pos), *, iostat=stat) ori_vec(i)
+                    if(stat>0) STOP "Error reading orientation value"
+                    ori_pos=ori_pos + 1
+                end if
+            end do
+        end if
         return
     end subroutine parse_ori_vec
     
@@ -241,5 +257,52 @@ module subroutines
         return
     end subroutine check_right_ortho
 
+    subroutine init_random_seed()
+        implicit none
+        integer, allocatable :: seed(:)
+        integer :: i, n, un, istat, dt(8), pid, t(2), s
+        integer(8) :: count, tms
+      
+        call random_seed(size = n)
+        allocate(seed(n))
+        ! First try if the OS provides a random number generator
+        open(newunit=un, file="/dev/urandom", access="stream", &
+             form="unformatted", action="read", status="old", iostat=istat)
+        if (istat == 0) then
+           read(un) seed
+           close(un)
+        else
+           ! Fallback to XOR:ing the current time and pid. The PID is
+           ! useful in case one launches multiple instances of the same
+           ! program in parallel.
+           call system_clock(count)
+           if (count /= 0) then
+              t = transfer(count, t)
+           else
+              call date_and_time(values=dt)
+              tms = (dt(1) - 1970) * 365_8 * 24 * 60 * 60 * 1000 &
+                   + dt(2) * 31_8 * 24 * 60 * 60 * 1000 &
+                   + dt(3) * 24 * 60 * 60 * 60 * 1000 &
+                   + dt(5) * 60 * 60 * 1000 &
+                   + dt(6) * 60 * 1000 + dt(7) * 1000 &
+                   + dt(8)
+              t = transfer(tms, t)
+           end if
+           s = ieor(t(1), t(2))
+           pid = getpid() + 1099279 ! Add a prime
+           s = ieor(s, pid)
+           if (n >= 3) then
+              seed(1) = t(1) + 36269
+              seed(2) = t(2) + 72551
+              seed(3) = pid
+              if (n > 3) then
+                 seed(4:) = s + 37 * (/ (i, i = 0, n - 4) /)
+              end if
+           else
+              seed = s + 37 * (/ (i, i = 0, n - 1 ) /)
+           end if
+        end if
+        call random_seed(put=seed)
+      end subroutine init_random_seed
 
 end module subroutines