Properties of Selenium under Tensile Stretch: A Molecular Dynamics Simulation

Philipp Baumgart, Katharina Berlin, Angelika Bogatyrev, John Kappelmaier, Ida Kemmer, Roland Rest, Vincent Salz, Jente Schienke, Emely Wiesner, Christina Oligschleger^*

Department of Applied Science, Bonn-Rhine-Sieg University of Applied Sciences, Rheinbach, Germany.

DOI: 10.4236/msce.2025.133003                

Abstract

Selenium is the third element of the group of chalcogens and is especially used in electrotechnical applications. Since many manufacturing processes use thin films of the material, knowledge of the mechanical properties is important for the detailed process layout. During the studies, the mechanical properties of trigonal selenium were determined. Therefore simulations of a uniaxial tensile test with a constant strain rate of 4.2⋅10^–6/fs using molecular dynamic methods were performed. Finally a great dependency of the mechanical properties on temperature, direction and strain was observed. In particular, anisotropy in fracture behavior can be demonstrated. This can be explained by the anisotropy of trigonal selenium. In addition, a density distribution was calculated, showing that the mechanical properties are influenced by diffusion processes occurring at increased temperatures. The changes in dynamic processes point to a ductile-brittle transition.

Keywords

Selenium, Tensile Test, Molecular Dynamics

1. Introduction

Selenium is the third element of the chalcogen group (6th main group in the periodic table). Chemically, it behaves similarly to sulfur. Thus, many compounds exist where both elements are used analogously. Physically selenium is often used as a semiconductor where it behaves as a p-type, like tellurium—the fourth element of this group. Additionally, it displays both photovoltaic action and photoconductivity [1].

There are at least eight different selenium allotropes claimed to exist in literature. This encompasses three monoclinic (α, β, γ) polymorphisms, one hexagonal (trigonal) crystalline, and two amorphous forms [2]. The structure and crystallinity of the element have major effects on both chemical and physical as well as the mechanical properties of the material and therefore impact the possible usage.

Presently, selenium is usually extracted from the anode slime resulting from copper refining, usually containing between 5% and 25% of the element. Additionally, it may also be recovered from past applications, such as xerographic photocopier drums. The material has many possible uses. Semiconductors in rectifiers and xerographic photoreceptors for photovoltaic cells and formerly photocopiers, pigments for the decoloration and coloration of some glasses and glazes yielding various colors like black, bronze, pink, purple, and maroon depending on the other pigments present, and metallurgical use as an alloying element [1].

The mechanical properties of selenium are most relevant for its application as a semiconductor due to them usually being made from multiple thin films [1]. With changing thermal and mechanical loads, these layers can crack and spall off due to the difference in their physical and mechanical properties like the density, thermal expansion and Young’s modulus [3]. Thus, it is important to simulate these properties and the possible behavior of these materials to better adjust the layer design, e.g., by including bond layers. Even though the possibility of running in situ tensile tests in the nanoscale exists for transmission electron microscopy, these are more focused on the spatial resolution of stress in the material. This serves to give possible reasons for changing physical and mechanical properties on the nanoscale but does not evaluate the fundamental one—the structure change on an atomic scale—in detail. However, it is the basis for the calculation of the strain, which may be obtained via electron diffraction by relating the present lattice parameters to the unstressed lattice [4].

For both the visualization and the more detailed evaluation of the occurring phenomena on an atomic scale, molecular dynamics (MD) can serve as an important tool.

2. Materials and Methods

2.1. Molecular Dynamics

In molecular dynamics the behaviour of molecules and atoms is described by simulation. As a basis for such simulations, statistical physical equations are used. The dynamics are mainly described using the second Newton’s law shown in Equation (1) [5].

$F_i=-{\nabla }_{i}U=m_i\cdot {\ddot{r}}_i$ (1)

Where $F_i$ equals the force (which corresponds to the negative gradient of the potential energy $U$ ), $m_i$ stands for the mass, and ${\ddot{r}}_i=a_i$ is the acceleration with regard to an atom of interest. Forces in between the molecules are depicted with “orce-fields” (potentials), these fields can be based on different empirically determined potential functions. These functions describe electrostatic, covalent, and van-der-Waals interactions. There is a huge variety of models which can be applied, depending on the molecules to be observed [5].

To correctly determine the position and acceleration of the molecules, the equation of motion is integrated over short periods of time typically in the range of femtoseconds. As a result, these functions are set up over time, and the information gained consists of: the position, angle, orientation and velocity of the molecule. What is not calculated in classical molecular dynamics are quantum mechanical effects. The accuracy and precision of the results are highly dependent on the used force field [5].

2.2. Potential

Interatomic potentials are used to balance the attractive and repulsive forces between neighboring atoms within a system. Thus ensuring that the system’s overall potential energy reaches its lowest possible value. These potentials are established by fitting functions to static properties derived from experimental studies or first-principle data. Experimental data often include measurements like cohesive energy, modulus values, lattice parameters, thermal expansion coefficients and vibrational data. With advances in in computational power and efficiency, significant progress has been made in developing various interatomic potentials for metals, alloys, semiconductors, ceramics, intermetallic compounds and organic materials. Generally, interatomic potentials are categorised into four types: pair potential, embedded atom method, Tersoff potential, and ionic solid potential [6].

For this work, a three-body potential was used, because short-ranged pair potentials are not usable for modeling covalently bonded systems [7] [8]. Thus, a three-body term was adopted which is a critical component for accurately modeling the interactions within selenium, particularly in both its molecular and crystalline forms. The three-body potential used in this study accounts for interactions among triplets of atoms, incorporating both radial distances and angular dependencies. The potential is structured to ensure that it can effectively replicate the geometries and binding energies of small selenium clusters and extend this accuracy to bulk selenium in its various phases with particular focus on structural, dynamic and thermodynamic properties of crystalline modifications [9] [10]. The potential energy $U$ is expressed as a sum of two-body and three-body terms, and represented in equation (2).

$U={\displaystyle \sum }_{i<j}{V}_2\left(r_{ij}\right)+{\displaystyle \sum }_{i<j<k}h\left(r_{ij},r_{kj},{\Theta }_{ijk}\right)+\text{cyclic}\text{permutations}$ (2)

Here, $V_2\left(r_{ij}\right)$ represents the two-body potential, which depends on the distance $r_{ij}$ between the two atoms $i$ and $j$ . The three-body interaction term $h\left(r_{ij},r_{jk},{\Theta }_{ijk}\right)$ is a function of the two distances $r_{ij}$ and $r_{jk}$ , and the angle ${\Theta }_{ijk}$ formed by these bonds at the atom $j$ .

2.3. Details of Simulation

In the following subsection the experimental setup and the observables are going to be explained.

2.3.1. Experimental Setup

In this experimental setup for the analysis of trigonal selenium, five temperatures are analysed under tensile testing: 100 K, 200 K, room temperature 298 K (approx. 300 K), 400 K, and 500 K. Before carrying out the tensile tests the sample comprising $N=384$ atoms was heated up to the desired temperature and afterward equilibrated. To provide stable temperature conditions a canonical ensemble was used, thus the number of atoms and the volume were kept constant and the temperature was kept constant and equilibrated via a velocity scaling scheme which ensures deviations of about 5%. In order to reduce surface effects we applied periodic bondary conditions. After this procedure, the actual tensile tests were carried out. The tensile tests were uniaxial, strain-controlled and conducted in all three directions $x,y$ and $z$ in individual analysis. The strain of the respective side length started with 1%, followed by 5% after that the strain was increased in 5% increments until the sample ruptured. All tensile tests are carried out a constant strain rate of 2.0⋅10⁻⁶/MDS≡4.219⋅10⁻⁶/fs, which is a typical value used in MD simulations [11]. When the target volume is reached with the desired tensile strain, the system is brought back into equilibrium as a canonical ensemble, i.e. we equilibrate the structures for additional 5000 to 100,000 MDS depending on the applied tensile strain.

2.3.2. Observables

For the evaluation the following different observables were considered with the overarching goal of analysing and describing the energy and structural changes of the trigonal selenium in short-range, micro-, and macroscopic scales during the tensile tests.

First, the (potential and total) energy of the system and the displacements of the atoms were observed during the MD simulations and recorded over the time (in molecular dynamic time steps-MDS). From this the elongation at break can be seen in a drastic change of energy, since bonds are broken due to separation of atoms. Additionally, the structures were monitored and visualized using VMD [12]. This allows conclusions to be drawn about the spatial arrangement of the atoms [13].

Changes within the short- and medium-range order of the structure are revealed by the pair distribution function $g\left(r\right)$ .

$g\left(r\right)=\langle \frac{n\left(r\right)}{4\pi r^2{\rho }_N\text{Δ}r}\rangle$ (3)

The pair distribution is averaged over all reference particles in the structure (ensemble mean), in which $n\left(r\right)$ is the number of particles at a distance $r$ within a spherical shell with surface area $4\pi r^2$ around the reference atom. $\text{Δ}r$ is the thickness of the spherical shell and ${\rho }_N$ is the particle number density. The time dependence of the structural changes were determined by evaluation of the van Hove correlation-function. In general the van Hove correlation-function describes the probability of an atom being within of a radius $r$ of a target atom, which gets analysed to assess changes in the short-range order. The self-consistent-part van Hove correlation-function describes the movement of a single atom in relation to its surroundings, providing insights into self-diffusion and movement patterns [14].

$G\left(r,t\right)=\frac{1}{N}\langle {\displaystyle \sum _{i,j=1}^N{\displaystyle \int \delta \left(r+r_i\left(0\right)-r^{\prime }\right)\delta \left(r^{\prime }-r_j\left(t\right)\right)\text{d}{r}^{\prime }}}\rangle$ (4)

To acquire information about material compression and elongation, the mean density distribution was analysed. The average density $\rho$ and its standard deviation $s_{\rho }$ were observed over time, allowing overall density changes to be recognised. In order to obtain more detailed data about the density distribution in the material on a smaller scale and thus identify local inhomogeneities, the densities of smaller subsystems and their distribution were also recorded. Finally, in order to further analyze material changes in terms of density, such as phase transition processes, the “ideal” density, at which the material is regularly distributed throughout the volume, is calculated and compared with the measured densities.

3. Results

During the tensile testing simulation the total energy [eV] is plotted as a function of time [MDS]. The different temperatures 100 K to 500 K are plotted individually, each with all three Cartesian directions ($x,y,z$ ). The applied tensile strains are 1%, 5%, 10%, 15% and 20% of the sample length. Different tendencies can be observed by comparing the obtained data.

3.1. Direction Dependency

The trigonal selenium shows different structural appearance depending on the regarded direction, see Figure 1.

(a) (b)

Figure 1. Structure of trigonal Se, parallel (left) and perpendicular to the chains (right), visualized using Molden [15].

This leads to a direction dependency in tensile strength which can be evaluated in Figures 2-4, in which the total energies are plotted vs. time. The direction-dependent energy curves are displayed in separate diagrams for the simulated temperatures from 100 K (top) to 500 K (bottom). Figure 2 represents the results at a strain of 10%, Figure 3 depicts our findings at a strain of 5%, and in Figure 4 the energy curves are shown at a strain of 1%, resp.: A general tendency can be observed for the x- and y-direction since they exhibit a similar behaviour for lower strains (1% and 5%). For higher strains the total energy drops more rapidly (breakage in system) for y- than for x-direction. Regardless of this the drop of total energy in z-direction appears first and thus after the shortest period of time. After the drop a steady state in energy is reached. According to that the total energy is lower for the y-direction and highest for the z-direction.

3.2. Temperature Dependency

Changing the temperature in a system can lead to great changes in the structural appearance [16] [17] and strength of that system. In Figure 5, the total energies are plotted versus time at an applied strain of 10%. The temperature-dependent energy curves are given in three diagrams: on top of the figure the results are presented from the MD-runs by applying the tensile stress in x-direction, in the middle we present energy curves from simulations using the strain in y-direction, and on bottom the respective results are visualized from application a strain in z-direction.

As it can be seen in Figure 5 there are some overall tendencies for increasing temperature. First, it suggests that the noise increases, the total energy is thus not as steady as for lower temperatures. Also, the drop in the total energy (breakage in system) appears earlier, especially for the z-direction. In general the total energy in the system decreases with increasing temperature. Additionally, the plots for the different directions seem to anneal for higher temperatures which can be observed even better by comparing Figures 2-4.

Similarly to the total energy we also investigated the behaviour of the internal stress/pressure which is built up in our simulations. Figure 6 shows the internal pressure over time at different temperatures for all three spacial directions for the strain of 10%.

Internal pressure is negative because it represents the system’s response to the applied (external) expansion of the structure. The overall shape is the same for all three directions. There is a rise in stress (more negative) first followed by a drop (less negative) and a steady state close to 0 GPa. The individual peak shifts to higher times and lower stress with decreasing temperature. The only difference can be seen in x-direction where the plots of the temperatures below 400 K display a less high drop back than the other two temperatures.

3.3. Strain Dependency

The strain dependency gets illustrated in Figure 7. Therefore the total energy is

Figure 2. Depiction of direction dependency ($\epsilon =10\%$ ); Abscissa: Time in MDS, Ordinate: Total energy in electron volt [eV].

Figure 3. Depiction of direction dependency ($\epsilon =5\%$ ); Abscissa: Time in MDS, Ordinate: Total energy in electron volt [eV].

Figure 4. Depiction of direction dependency ($\epsilon =1\%$ ); Abscissa: Time in MDS, Ordinate: Total energy in electron volt [eV].

plotted against the time for all considered strains at a temperature of 200 K. The diagramme on top depicts the influence of strains in x-direction, the figure in the middle shows the results for strains aplied in y-direction, and on the bottom we present the results from the simulations with strain in z-direction. First, for the x- and y-direction the initial energy drop can be observed after a energetic maximum at ca. 33,000 MDS during the run with an application of an tensile strain with 10% elongation. In the cases when the elongations of 10%, 15% and 20% are applied, significant differences can be mentioned regarding the curve progression. After the initial energy drop, whereby in this case a direction dependence can be observed between x- and y-direction. In x-direction the total energy further increases with increasing strain, for the strains from 5% to 20% the final energies range from −2.38 eV to −2.35 eV.

Whereas in y-direction the curve progression initially decreases to a certain value and remains constant regardless of the strain (the final energies are around −2.38 eV). In contrast to the results gained in x- and y-direction the simulations in z-direction reveal a different behaviour. The energetic drop occurs even after a maximum strain of 5% is reached. This indicates structurals relaxation, but the final energy, with a value of −2.375 eV, is higher than that achieved by tensile strain in the x- or y-direction. For the tensile strains of 10%, 15% and 20% applied in z-direction the energy curves peak at slightly less than 30,000 MDS and each exhibit several relaxation levels with final energies ranging from −2.365 to −2.36 eV.

Figure 5. Depiction of temperature dependency ($\epsilon =10\%$ ); Abscissa: Time in MDS, Ordinate: Total energy in electron volt [eV].

Figure 6. Depiction of temperature dependency ($\epsilon =10\%$ ); Abscissa: Time in MDS, Ordinate: Pressure [GPa].

3.4. Density Distribution

As already mentioned a further evaluation of the obtained data was performed to determine the densities $\rho$ and the standard deviations of the densities $s_{\rho }$ . By taking the different strains and directions into account, the temperature dependence of the densities and their standard deviations are given in Table 1. Compared to an ideal density ${\rho }_{\text{ideal}}$ -under the assumption that the atoms are regularly spreaded in the volume- the densities in the simulated structures are reduced what is typical for less ordered systems. At all simulated temperatures, the densities for the different

Figure 7. Depiction of strain dependency (as an example at $T=200\text{K}$ ); Abscissa: Time in MDS, Ordinate: Total energy in electron volt [eV].

Table 1. Taking into account the different strains and directions, we show the temperature dependency of both the densities $\rho$ and standard devations $s_{\rho }$ , resp. In the last column the ideal density for equally distributed Se-atoms are given. All densities are in g/cm³.

Temperature	x-direction	y-direction	z-direction	Ideal density
strain	$\rho$ $s_{\rho }$	$\rho$ $s_{\rho }$	$\rho$ $s_{\rho }$	${\rho }_{\text{ideal}}$
K
0.01	4.62500 0.76547	4.60156 0.80127	4.63281 0.80157	4.80217
0.05	4.35938 0.70417	4.39062 0.85453	4.36719 0.70844	4.61923
0.10	4.12500 0.91430	4.15625 1.11058	4.13281 1.09773	4.40927
K
0.01	4.60156 0.75612	4.61719 0.88274	4.60156 0.80127	4.80217
0.05	4.35938 0.79779	4.37500 0.71261	4.38281 1.07071	4.61923
0.10	4.14062 0.94567	4.14062 1.30943	4.13281 1.14987	4.40927
0.15	4.01562 0.60576	3.98438 1.36350	4.00000 1.08253	4.21756
0.20	3.83594 0.78588	3.89844 1.40640	3.81250 1.19406	4.04183
K
0.01	4.60938 0.84071	4.63281 0.74079	4.62500 0.82443	4.80217
0.05	4.36719 0.88495	4.36719 0.89810	4.36719 1.02027	4.61923
0.10	4.12500 0.91430	4.15625 1.33134	4.11719 1.03736	4.40927
0.15	3.95312 1.08512	3.96094 1.27569	3.99219 1.14732	4.21756
0.20	3.82812 1.30270	3.86719 1.21913	3.81250 1.25779	4.04183
K
0.01	4.60156 0.82528	4.61719 0.89155	4.61719 0.86486	4.80217
0.05	4.36719 0.80157	4.37500 0.87500	4.40625 0.84722	4.61923
0.10	4.12500 0.92702	4.13281 1.08700	4.10938 0.95798	4.40927
0.15	3.96094 1.18026	3.97656 1.30931	4.00000 1.17925	4.21756
K
0.05	4.37500 1.10043	4.37500 0.98821	4.39844 0.80127	4.61923
0.10	4.08594 1.14412	4.12500 1.12152	4.11719 1.00289	4.40927
0.15	3.96094 0.89401	3.97656 1.04743	3.99219 0.89483	4.21756
0.20	3.78125 1.12803	3.84375 1.16885	3.82812 1.20292	4.04183

tensile strains have similar values with respect to the applied directions. Compared to the ideal densities, the calculated densities decrease to a similar extent with increasing tensile elongation.

Furthermore the density distribution can be examined with the van Hove correlation-function, which is visualised in Figure 8 and Figure 9. To do this, the atomic distance is plotted over time; the height of the peaks is colour-coded. We present the results of the simulations at two different temperatures, in this case 100 K and 500 K, each applying a maximum strain of 10% in x-direction. The lowest band (at approximately 2.3 Å) belongs to the nearest neighbour distance. This distance is maintained even at higher temperatures and tensile strains and indicates a pronounced short-range order. Additional bands indicate larger but regular distances between the atoms, this can be seen particularly at 100 K. With increasing temperature and time a general widening of the bands can be observed. At 500 K we find broad bands, indicating a larger spread in atomic distances. This result is in agreement with the observed rather large density deviations from Table 1.

Figure 8. Van Hove distribution at 100 K in x-direction by $\epsilon =10\%$ .

Figure 9. Van Hove distribution at 500 K in x-direction by $\epsilon =10\%$ .

4. Discussion

As described in chapter 3.1 there is a great dependency on the direction of the applied tensile strain. The first aspect that needs to be addressed is the anisotropy of the trigonal selenium. Therefore the general structure of the material is shown in Figure 1. In the cartesian z-direction the mechanical properties are dominated by the covalent bonds of the atomic chains, while the intermolecular bonds, in particular van der Waals interactions, determine the mechanical properties in the x- and y-direction. Therefore, higher tensile strengths as well as Young’s modulus result along the atomic chains due to the higher bonding energies of the covalent bonds in comparison to the intermolecular interactions between the molecular chains [18].

Also higher values regarding the total energy can be observed in z-direction until the breakage of the material compared to the other spatial directions. Here, the amount of total energy increases due to the application of force, which favors an increasing amount of total energy. Whereby also a strong dependence on temperature needs to be mentioned, which is explained in more detail below. Additionally the drop in the plot of the y- and x-direction can only be observed for a strain of 10% and furthermore the drop in the z-direction appears earlier in time. This reveals a higher pronounced brittleness of the material in z-direction.

The anisotropy is also visible in the fracture behaviour which is demonstrated in Figure 10. In x-direction corner cracks can be observed whereas plane cracks can be observed in y-direction. This can be related to the angle the force is applied on the selenium lattice. In y-direction, the force is applied perpendicular to the lattice planes, which causes fracture of whole plane. In x-direction, the force is applied in an angle of approximately 125˚ which causes diagonal fracture. In this small sized simulation, this results in the disruption of a corner. In contrast inz-direction no characteristic fracture behaviour becomes noticeable due to the fracture of covalent bonded atoms. This behaviour was already observed in literature where also a irregular fracture morphology was observed by applying uni-axial force in chain direction [18].

Figure 10. Tensile test simulation for all directions at $\epsilon =10\%$ . Timesteps for taking the images are stated in MDS at 100 K.

To further analyse the temperature dependencies the observations of chapter 3.2 need to be regarded. The value of the total energy increases with increasing temperature whereas the time for the initial breakage of the material occurs earlier. This would indicate an increasing embrittlement of the material as the temperature rises. Besides at 500 K no further effects regarding the anisotropy become evident. In addition, a stronger noise of the energy curves becomes apparent with increasing temperature. These aspects can be addressed to diffusion processes of the material. The diffusion gets apparent by the van Hove distribution visualised by Figure 8 and Figure 9. At 100 K narrow bands are visible, indicating an ordered structure in which the atoms oscillate around their equilibrium positions. On the other hand, at 500 K, broader bands can be observed, which are still clearly separated at the beginning of the tensile test (i.e. at early time of the simulation), but are no longer clearly separated and overlap in course of the tensile elongation and the relaxation phase, i.e. at later times of the MD simulation. The probability that the atoms are at a certain distance from each other is therefore widened, which implies a favoured diffusion. In addition, the clearly separated bands of the material at 100 K indicate the intact crystalline structure. On closer inspection, shifts in the bands associated with chain spacing and longer atomic distances (i.e. distances between 3 and 5 Å) occur at around 33000 MDS. This can be attributed to the fracture of the material. At a temperature of 500 K this cannot be observed also regarding the fact, that no clear fracture is visible at 500 K. However, the broadness of the structure bands exhibit a loss in crystallinity. As explained above, Figure 6 shows the relaxation behaviour. While at lower temperatures the selenium reveals brittle behaviour displayed by the abrupt drops in stress when cracks occur. At higher temperatures those abrupt drops in stress don’t occur, since the cracking is circumvented by diffusion of atoms. As the cracks in the x-axis usually are only partial, the system can retain some of the stress beyond the pull phase. While at lower temperatures this stress is retained, at higher temperatures diffusion allows for an reduction of the stress. All plots indicate that with increasing temperature the behaviour from brittle to ductile changes gradually.

Analysing the findings of chapter 3.3 also the strain dependency needs to be discussed. The lowest strain of 1% is not able to cause a break in any direction. Likewise no breakage can be observed for strains of 5% in x- and y-direction but in z-direction, whilst a breakage becomes visible for all directions for higher strains than 5%. The time of fracture is reached contemporaneously in all directions. The curve progression in x- and y-direction appears to be direction dependent and not strain dependent since the curves show a similar behaviour. Whereas in x-direction the increase of the total energy needs to be addressed after the initial energy drop. These differences could be related to the fracture morphology previously mentioned. In z-direction the strain dependence becomes more apparent after the initial energy drop. The phenomena can be explained by a relaxation processes occurring after the pulling. The lattice doesn’t fracture completely so there is still an amount of stress, it can bear. This stress is reduced after the phase of pulling while holding the load constant.

Further analysis of the density distribution 3.4 as well as visualization of the structures (see Figure 10) lead to the conclusion that under certain conditions partially amorphous structures can be formed from crystals. Since the density of amorphous Se ranges from 4.25 to 4.45 g/cm³, these values can easily be achieved by applying tensile strains [19]. Especially tensile strains in z-direction, i.e. parallel to the covalently bonded atomic chains, cause regions with amorphous structures. We can therefore conclude that the application of certain tensile strains is suitable for generating partially amorphous structures without going through the melting phase. Different degrees of amorphization have an impact on e.g. vibrational properties [16] [20]. What is not investigated here is the influence of these partly crystalline, partly amorphous structures on dynamic, thermodynamic and electronic properties. These types of structures could lead to interesting applications and are the subject of further research [21] [22].

5. Conclusions

The results of the simulation enable a deeper understanding of the anisotropic material behavior of selenium. These could be described using a uniaxial tensile test. Besides the anisotropic behaviour also dependencies regarding the temperature could be analysed to further understand the diffusion-processes of crystalline selenium with increasing temperature. To a limited extent it was also possible to verify the results with findings from the literature.

The angle between the direction of load and the lattice planes is the most important factor for the fracture behaviour and the bearable forces. The highest tensile strength and Young’s modulus are observed in the z-direction due to the covalent forces acting in this direction. In x-direction, the simulation shows slightly higher strength than in y-direction because the shift-fracture of the lattice planes allows a small uptake of energy while the catastrophic fracture of the whole plane does not.

With increasing temperature, the strength of the material decreases, especially in the z-direction. The tensile test behaviour for the different load directions converges at 500 K while it can be clearly distinguished at 100 K. Also diffusion processes may take place, causing deformation without an increase of stress acting on the system at higher temperatures.

The strain applied has no influence on the strength of the material. But in the y-direction a higher strain does not lead to further stress build-up, while in z- and especially in the x-direction, further stress is generated as the strain increases. This is caused by the difference between the complete fracture of the whole lattice plane in the y-direction compared to the partial fracture in x- and z-direction.

Future laboratory tests would be recommended to verify the simulation data. In addition, varying loads to further investigate the material behaviour, such as pressure loads, would be of interest for further investigations. Additionally, comparisons with other configurations of the element selenium could also lead to meaningful results.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Butterman, W. and Brown, Jr., R.D. (2004) Selenium Mineral Commodity Profiles. US Geological Survey. http://pubs.usgs.gov/of/2003/of03-018/of03-018.pdf
[2]	Greenwood, N. and Earnshaw, A. (1984) Chemistry of the Elements. Elsevier Science and Technology. https://www.sciencedirect.com/book/9780750633659/chemistry-of-the-elements
[3]	Abadias, G., Chason, E., Keckes, J., Sebastiani, M., Thompson, G.B., Barthel, E., et al. (2018) Review Article: Stress in Thin Films and Coatings: Current Status, Challenges, and Prospects. Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films, 36, Article 020801. https://doi.org/10.1116/1.5011790
[4]	Hÿtch, M.J. and Minor, A.M. (2014) Observing and Measuring Strain in Nanostructures and Devices with Transmission Electron Microscopy. MRS Bulletin, 39, 138-146. https://doi.org/10.1557/mrs.2014.4
[5]	Haberlandt, R., Fritzsche, S., Peinel, G. and Heinzinger, K. (2012) Molekulardynamik. Springer Vieweg.
[6]	Pal, S. and Reddy, K.V. (2024) Atomistic Simulation. In: Molecular Dynamics for Materials Modeling, CRC Press, 1-22. https://doi.org/10.1201/9781003323495-1
[7]	Stillinger, F.H. and Weber, T.A. (1985) Computer Simulation of Local Order in Condensed Phases of Silicon. Physical Review B, 31, 5262-5271. https://doi.org/10.1103/physrevb.31.5262
[8]	Stillinger, F.H., Weber, T.A. and LaViolette, R.A. (1986) Chemical Reactions in Liquids: Molecular Dynamics Simulation for Sulfur. The Journal of Chemical Physics, 85, 6460-6469. https://doi.org/10.1063/1.451426
[9]	Oligschleger, C., Jones, R.O., Reimann, S.M. and Schober, H.R. (1996) Model Interatomic Potential for Simulations in Selenium. Physical Review B, 53, 6165-6173. https://doi.org/10.1103/physrevb.53.6165
[10]	Oligschleger, C., Jones, R.O., Reimann, S.M. and Schober, H.R. (2020) Erratum: Model Interatomic Potential for Simulations in Selenium [Phys. Rev. B 53, 6165(1996)]. Physical Review B, 102, Article 99901. https://doi.org/10.1103/physrevb.102.099901
[11]	Koh, S.J.A. and Lee, H.P. (2006) Molecular Dynamics Simulation of Size and Strain Rate Dependent Mechanical Response of FCC Metallic Nanowires. Nanotechnology, 17, 3451-3467. https://doi.org/10.1088/0957-4484/17/14/018
[12]	Humphrey, W., Dalke, A. and Schulten, K. (1996) VMD: Visual Molecular Dynamics. Journal of Molecular Graphics, 14, 33-38. https://doi.org/10.1016/0263-7855(96)00018-5
[13]	Frenkel, D. and Smit, B. (2002) Monte Carlo Simulations. In: Understanding Molecular Simulation, Elsevier, 23-61. https://doi.org/10.1016/b978-012267351-1/50005-5
[14]	Hopkins, P., Fortini, A., Archer, A.J. and Schmidt, M. (2010) The Van Hove Distribution Function for Brownian Hard Spheres: Dynamical Test Particle Theory and Computer Simulations for Bulk Dynamics. The Journal of Chemical Physics, 133, Article 224505. https://doi.org/10.1063/1.3511719
[15]	Schaftenaar, G., Vlieg, E. and Vriend, G. (2017) Molden 2.0: Quantum Chemistry Meets Proteins. Journal of Computer-Aided Molecular Design, 31, 789-800. https://doi.org/10.1007/s10822-017-0042-5
[16]	Bermejo, F.J., García-Hernández, M., Mompeán, F.J., MacMorrow, D. and Martinez, J.L. (1995) Nature of the First Diffraction Peak in Glassy Selenium. Physical Review B, 51, 11932-11935. https://doi.org/10.1103/physrevb.51.11932
[17]	Kirchhoff, F., Kresse, G. and Gillan, M.J. (1998) Structure and Dynamics of Liquid Selenium. Physical Review B, 57, 10482-10495. https://doi.org/10.1103/physrevb.57.10482
[18]	Qin, J., Sui, C., Qin, Z., Wu, J., Guo, H., Zhen, L., et al. (2021) Mechanical Anisotropy in Two-Dimensional Selenium Atomic Layers. Nano Letters, 21, 8043-8050. https://doi.org/10.1021/acs.nanolett.1c02294
[19]	Reyes-Retana, J.A. and Valladares, A.A. (2010) Structural Properties of Amorphous Selenium: An AB Initio Molecular-Dynamics Simulation. Computational Materials Science, 47, 934-939. https://doi.org/10.1016/j.commatsci.2009.11.026
[20]	Oligschleger, C., Facius, C., Kutz, H., Langen, C., Thumm, M., von Brühl, S., et al. (2009) Molecular Dynamics Simulation of Structural and Dynamic Properties of Selenium Structures with Different Degrees of Amorphization. Journal of Physics: Condensed Matter, 21, Article 405402. https://doi.org/10.1088/0953-8984/21/40/405402
[21]	Ksiażek, K., Wacke, S., Górecki, T. and Górecki, C. (2007) Effect of the Milling Conditions on the Degree of Amorphization of Selenium by Milling in a Planetary Ball Mill. Journal of Physics: Conference Series, 79, Article 012037. https://doi.org/10.1088/1742-6596/79/1/012037
[22]	Hegedüs, J. and Kugler, S. (2005) Growth of Amorphous Selenium Thin Films: Classical versus Quantum Mechanical Molecular Dynamics Simulation. Journal of Physics: Condensed Matter, 17, 6459-6468. https://doi.org/10.1088/0953-8984/17/41/016