# Bending Deformation Driven by Molecular Rotation

Pedro A. Santos-Florez  Department of Physics and Astronomy, University of
Nevada, Las Vegas, NV 89154, USA  Shinnosuke Hattori  Advanced Research
Laboratory, Technology Infrastructure Center, Technology Platform, Sony Group
Corporation, 4-14-1 Asahi-cho, Atsugi-shi 243-0014, Japan  Qiang Zhu
<EMAIL_ADDRESS>Department of Physics and Astronomy, University of Nevada, Las
Vegas, NV 89154, USA  Department of Mechanical Engineering and Engineering
Science, University of North Carolina at Charlotte, Charlotte, NC 28223, USA

###### Abstract

In recent years, certain molecular crystals have been reported to possess
surprising flexibility by undergoing significant elastic or plastic
deformation in response to mechanical loads. However, despite this
experimental evidence, there currently exists no atomistic mechanism to
explain the physical origin of this phenomenon from numerical simulations. In
this study, we investigate the mechanical behavior of three naphthalene
diimide derivatives, which serve as representative examples, using direct
molecular dynamics simulations. Our simulation trajectory analysis suggests
that molecular rotational freedom is the key factor in determining a crystal’s
mechanical response, ranging from brittle fracture to elastic or plastic
deformation under mechanical bending. Additionally, we propose a rotation-
dependent potential energy surface as a means to classify organic materials’
mechanical responses and identify new candidates for future investigation.

##  I INTRODUCTION

While most molecular crystals are brittle, there exists a class of compliant
organic crystals that can bend under a large mechanical stress loading Naumov
_et al._ (2015); Saha  _et al._ (2018). Since early 2000, a growing number of
mechanically flexible organic crystals have been reported experimentally Reddy
_et al._ (2006); Takamizawa and Miyamoto (2014); Ghosh and Reddy (2012); Panda
_et al._ (2015); Krishna  _et al._ (2016); Yadav  _et al._ (2019); Raju  _et
al._ (2018); Mishra and Sun (2020); Zhang  _et al._ (2021). In general, the
mechanical response of an organic solid depends on both the molecular
substance and crystal packing. A remarkable example is shown in Fig. 1. Three
crystals, made of similar molecules from naphthalene diimide derivatives, were
found to exhibit distinct responses from brittle fracture to compliant
deformation with either reversible (elastic) or irreversible (plastic)
characteristic Devarapalli  _et al._ (2019). The flexible nature of organic
materials is vital for a variety of applications, e.g., high-performance
modular solar cells Root  _et al._ (2017), actuators Li  _et al._ (2019),
photochemistry Mutai  _et al._ (2021), fluorescence Di  _et al._ (2022);
Hayashi and Koizumi (2016), electronics Wang  _et al._ (2019); Samanta  _et
al._ (2023), optics Liu  _et al._ (2018), as well as drug tabulation Sun
(2017).

Figure 1:  The simulated bending on three different materials based on
naphthalene diimide derivatives. (a) brittle Pr (50.3$\times$7.0$\times$6.8
nm3), (b) elastic Et (50.7$\times$6.4$\times$6.6 nm3) and (c) elastic/plastic
Me (50.2$\times$6.4$\times$6.9 nm3). These three crystals consist of very
similar molecules that differ only in the side groups. In the left panel, the
initial and finally deformed configurations are colored by the accumulated
molecular rotation ($\alpha$) along the $x$-axis. The corresponding molecules
and the definition of rotation angles are shown in the right panel.

In the recent years, various computational techniques have been introduced to
characterize the mechanical properties of molecular systems Devarapalli  _et
al._ (2019); Wang and Sun (2018, 2019, 2020); Ootani and Kubo (2022);
Matveychuk  _et al._ (2022). They include topological analysis, elastic
properties calculation Wang and Sun (2019), and shear/tensile simulations
Devarapalli  _et al._ (2019); Ootani and Kubo (2022). These techniques are
successful in identifying brittle materials. Within an interlocked environment
in crystal packing, molecular motions in those materials are largely
restricted, resulting a brittleness under bending Wang and Sun (2019). On the
other hand, some materials are featured by a strong anisotropy with plausible
slip planes Wang and Sun (2019); Bryant  _et al._ (2018). Therefore, these
materials become compliant over a broad range of applied stress along some
crystallographic directions. However, all available techniques fail to explain
the difference between the elastic and plastic materials. While there have
been plenty of studies on the bending of metals Zhu  _et al._ (2012); Zhang
_et al._ (2014); Nöhring  _et al._ (2016); Zhuo and Beom (2018); Katakam and
Yedla (2020); He and Ma (2022), to our knowledge, no attempts have been made
to directly simulate the bending of organic materials at the atomistic level.

Among the compliant crystals, ductile materials are often favored in
engineering applications Sun (2017). Hence, researchers attempted to use the
well established dislocation theory to explain the observed plasticity on
organic materials Reddy  _et al._ (2006); Saha  _et al._ (2018). Similar to
the plastic deformation in ductile metals, it was proposed that mechanical
shearing can occur via the slippage of dislocated molecular layers on the
molecular crystals with a layered packing Reddy  _et al._ (2010); Wang and Sun
(2018). Using these slip planes, a bending model was proposed to explain the
underlying mechanism Reddy  _et al._ (2005a). Although the dislocation is not
uncommon in molecular crystals Mathew  _et al._ (2013); Mathew and Picu
(2013); Olson  _et al._ (2018), there has been no direct experimental evidence
to support that dislocation pre-exists or appears in the organic crystals
under bending. Furthermore, this mechanism fails to explain the observed
large-scale elastic deformation. In fact, two crystals in Fig. 1b-c can
undergo either elastic or plastic deformation despite the apparent similarity.
Clearly, our current understanding on molecular bending remains limited.

In this work, we present our efforts in questing the molecular bending
mechanism through atomistic simulation. To achieve this goal, we started by
developing a simulation protocol that can directly model the bending of
organic crystals at the atomic level. Next, the simulation results were
carefully analyzed to classify and understand the atomistic mechanism of
materials-dependent deformations from brittle fracture to elastic or plastic
deformation. Furthermore, we demonstrated that molecular rotational freedom is
the key factor in determining a crystal’s mechanical response. Finally, we
introduced a rotation-dependent potential energy surface as a means to
classify organic materials’ mechanical responses and identify new candidates
for future screening of new mechanically flexible organic crystals.

##  II Methods

###  II.1 Crystal Structures of Three Systems

In this study, we focused on three systems consisting of naphthalene diimide
derivatives as shown in Fig. 1. The three molecules share the same backbone
while differing only in the side chains. The brittle crystal consists of the
molecules with the propyl group, featured by the orthorhombic space group
$Pbca$ with one molecule in the asymmetric unit. On the other hand, the
elastic/plastic crystals have the ethyl/methyl groups, both adopting the
monoclinic space group $P2_{1}/c$ with half a molecule in the asymmetric unit.
For convenience, we follow the previous literature Devarapalli  _et al._
(2019) to name these systems according to their molecular functional groups
(i.e., Pr, Et, Me). In all three cases, the weak interaction plane formed by
alkyl groups is (001). In Fig. 2, each molecule in the unit cell is colored by
the alignment along the $y$-axis. Clearly, the overall molecular packing in
the brittle-Pr crystal are more complex. Since there exist eight different
types of molecular alignments due to the $mmm$ symmetry operations, the Pr
crystal has molecules aligned in different ways within the same (001) layer.
On the contrary, there are only two types of molecular alignments in the Et/Me
crystals. And the (001) layer in Et/Me crystals has all molecules aligned in
the same direction. Table LABEL:table1 summarizes the crystallographic
information of three molecular crystals.

Table 1: The crystallographic information of three molecular crystals. Among them, Pr denotes the brittle crystal with the CSD refcode of DAHLOQ; Et is the elastic crystal with the CSD refcode of BIYRIM01; and Me is the plastic crystal with the CSD refcode of DAHMUX. The column of size list the number of molecules in each model. System |  Symmetry |  Size |  $a$ (Å) |  $b$ (Å) |  $c$ (Å) |  $\beta$ (∘)
---|---|---|---|---|---|---
Pr | $Pbca$ | 8 | 6.96 | 17.24 | 27.58 | 90.0
Et | $P2_{1}/c$ | 2 | 4.84 | 7.74 | 18.32 | 90.1
Me | $P2_{1}/c$ | 2 | 4.62 | 8.02 | 17.02 | 94.0
Figure 2:  The crystal structures of (a) Pr, (b) Et (c) Me systems.

###  II.2 Atomistic Modeling of Bending

To directly simulate the bending of organic crystal, we employed a three-point
bending model within a partial periodic boundary condition based on the LAMMPS
package Plimpton (1995) at room temperature. In our calculation, we performed
non-equilibrium molecular dynamics (MD) simulation by applying the indentation
on the molecular slab model (see Fig. 3). Both $x$ and $y$-axes are under the
constraint of periodic boundary conditions, while the $c$-axis is not
periodic. We rotated the crystal structures with the matrix of [[0,0,1],
[0,-1,0], [1,0,0]], and then built the super cell slab models with sufficient
vacuum separation. The slab correction was applied to remove the slab-slab
interactions from the periodic images. Due to the non-triclinic box
restriction on the computation of slab correction, the $\beta$ angles for the
slabs of Et and Me were to be set to 90∘, which are slightly different from
the ideal values. However, this compromise should not change the results
largely.

Figure 3:  The schematic setup of a bending simulation model.

Along the non-periodic $z$-axis, a cylindrical indenter with the radius of 30
Å is applied on top of the slab center in the unit cell. To mimic two other
contacting points in the three-points bending simulation, the last layer of
molecules in the bottom region were frozen in the entire simulation. In
addition, the first columns of molecules on both left and right sides of the
unit cell are defined as the border. The rest atoms not belonging to the
frozen and border groups are set to the mobile group that can move freely. To
ensure a sufficient heat bath, we first performed Langevin thermostat
Schneider and Stoll (1978) on both mobile and border groups, followed by a
second thermal equilibration on only the border atoms. The fully equilibrated
sample was used to perform three-points bending simulation with only the
border atoms being under the Langevin thermostat to mimic the external
temperature reservoir. Upon bending, an indenter was used to push into the
simulation slab in a flow with the rate of 10 m/s. When the system reaches the
maximum indentation depth, the indenter was kept for 300-500 ps to allow the
system achieves thermal equilibrium. Afterwards, the indenter will move upward
with the previous rate to mimic the release of indenter process. To check the
dependence of indentation rate, we also varied the rates from 2-50 m/s. We
found that these rates roughly led to similar results. However, a rate faster
than 200 m/s may result in nonphysical phase transition for the Me sample. It
is also possible that the change of indenter shape, size and temperature may
change the results significantly. These factors will be the subject of future
work.

###  II.3 Force Field Choices and Benchmark

To reliably simulate the deformation of organic crystals at the atomistic
level, it is crucial to choose an accurate interatomic force field model. In
this work, we developed a computational pipeline to automate the generation of
molecular force fields from AmberTools20 Case  _et al._ (2020), based on the
General Amber Force Field (GAFF) Wang  _et al._ (2004) framework with atomic
charges using semi-empirical AM1-BCC method Jakalian  _et al._ (2000). To
confirm that the simulation results are not due to the artifact of force field
choices, we repeated the simulations using the OpenFF-toolkit Wagner  _et al._
(2023) with different parametrization protocol Boothroyd  _et al._ (2023);
Shirts  _et al._ (2017), as well as Density Functional based Tight Binding
(DFTB) Hourahine  _et al._ (2020). The OpenFF-toolkitWagner  _et al._ (2023)
was employed to generate the OpenFF model by assigning atom types based on
direct chemical perception, utilizing an atom-by-atom assignment approach
through the use of SMIRKS (SMILES Arbitrary Target Specification) patterns.
The OpenFF 2.0.0 (Sage) force field Boothroyd  _et al._ (2023) was adopted for
atom typing, and the ParmEd package Shirts  _et al._ (2017) was employed for
input file format conversion, ensuring compatibility with various molecular
simulation engines and the accurate representation of molecular topologies. In
the DFTB model, we used the DFTB+ code Hourahine  _et al._ (2020) with the
inclusion of van de Waals dispersion based on the Tkatchenko-Scheffler method
Tkatchenko and Scheffler (2009).

Table LABEL:table2 lists the computed equilibrium cell parameters with
different types of force fields at both zero and room temperatures. Clearly,
the GAFF model, as well as other models, yield similar results that are
comparable with the experimental values Devarapalli  _et al._ (2019). For
convenience, we will mainly employ the GAFF model in our following
simulations.

Table 2: The comparison of different force field models in describing the equilibrium cell parameters.  System |  Cell |  Experiment |  CVFF |  GAFF |  OpenFF |  DFTB
---|---|---|---|---|---|---
| (Å) | 300 K | 300 K | 300 K | 300 K | 0 K
|  | Ref. Devarapalli  _et al._ (2019) | Ref. Devarapalli  _et al._ (2019) | this work | this work | this work
Pr | $a$ | 6.96 | 7.54 | 7.30 | 7.35 | 6.69
$b$ | 17.24 | 16.99 | 17.41 | 17.62 | 17.08
$c$ | 27.58 | 28.14 | 27.09 | 27.20 | 27.74
Et | $a$ | 4.84 | 5.02 | 5.07 | 4.93 | 4.53
$b$ | 7.74 | 7.66 | 7.79 | 7.86 | 7.81
$c$ | 18.32 | 19.88 | 19.07 | 19.05 | 18.52
Me | $a$ | 4.62 | 4.60 | 4.58 | 4.50 | 4.29
$b$ | 8.02 | 7.87 | 8.28 | 8.19 | 8.02
$c$ | 17.02 | 18.66 | 18.40 | 17.82 | 16.69

##  III Computational Results and Discussions

To make a fair comparison, we set up all model sizes close to 50.0 $\times$
7.0 $\times$ 7.0 nm3 as summarized in Table LABEL:table3. For each system, we
also added the vacuum of 120 Å to allow the materials bend sufficiently. In
addition, we considered two kinds of Me models, including (i) the supercell
after the isobaric-isothermal (NPT) equilibration; and (ii) the supercell with
the experimental cell parameters. Although these two configurations only
differ slightly, it has been found they led to different elastic/plastic
deformation processes in the subsequent bending simulation. All supercell slab
models were then used to perform the three-points bending simulation as
illustrated in Fig. 3 with an indentation rate of 10 m/s under 300 K. For each
system, we ran the indentation simulation for multiple times to determine the
maximum indentation depth (5-20 nm) that leads to the formation of crack.
Before the maximum indentation depth is reached, we also continued the
simulation by releasing the indenter with the same rate to check if the
deformation process is reversible.

Table 3: The details of slab models used in the bending simulation. The column of size lists the number of molecules in each model. System |  Deformation |  Supercell | Size |  $a$ (Å) |  $b$ (Å) |  $c$ (Å)
---|---|---|---|---|---|---
Pr | brittle | $18\times 4\times 5$ | 5760 | 503.2 | 69.9 | 70.6
Et | elastic | $27\times 4\times 5$ | 6480 | 508.5 | 63.6 | 74.7
Me | elastic | $29\times 8\times 15$ | 6960 | 501.6 | 65.2 | 86.5
Me | plastic | $30\times 8\times 15$ | 7200 | 510.6 | 64.2 | 85.1

###  III.1 Direct Identification of Deformation Characteristics

Fig. 4 summarizes the simulated evolution of potential energy as a function of
indentation depth for all three materials. Encouragingly, our calculations
produced a sequence of deformations (including brittle fracture, elastic
deformation and plastic bending) that are similar to the previous experimental
observations Devarapalli  _et al._ (2019). First, Pr is clearly brittle as
evidenced by the abrupt drop of energy in Fig. 4a, which is consistent to the
appearance of crack pattern in Fig. 1a when the indenter reaches 3.5 nm. On
the other hand, Et is more complaint with a maximum indentation of 6.2 nm.
Applying further loading would lead to the formation of crack as well. If we
release the indentation before Et reaches 6.2 nm, the model roughly returned
to the original state. Therefore, this deformation is elastic. Interestingly,
Me can survive under more than 10 nm indentation without breaking under two
different setups. For the slab after a full NPT equilibration, it bends
elastically, as evidenced by the reversible energy versus indentation depth
relation (denoted as Me-elastic in Fig. 4b). When the slab has a small strain
in the initial configuration (see Table LABEL:table3), its energy curves upon
loading and unloading are no longer reversible. Compared to Me-elastic, this
sample achieves a lower energy when it approaches the maximum indentation
depth upon loading. When the indentation is released, it does not return to
the original state, but maintains a relatively higher energy. Therefore, the
whole deformation process is irreversible and plastic. The sample will be
called Me-plastic from now on.

Figure 4:  The evolution of average molecular potential energy as a function
of indentation depth upon (a) loading and (b) unloading.

To our knowledge, all previous computational studies were limited to indirect
simulations of tensile and shear tests Devarapalli  _et al._ (2019); Wang and
Sun (2019); Ootani and Kubo (2022); Matveychuk  _et al._ (2022). Here, our
calculations provide the first direct atomistic modeling on the experimentally
observed bending deformations. Compared to the simulation results, the elastic
and plastic samples are found to bend more significantly in real experiments
Devarapalli  _et al._ (2019). This is because that the material along $x$-axis
under the actual bending test can shrink to release the tensile stress.
However, our simulation model still obeys the periodic boundary condition
along the $x$-axis. Hence we expect that the degree of bending in simulation
is underestimated as compared to the real situation.

Figure 5:  The simulated distribution of accumulated rotational angles (with
respect to the initial configurations) for all materials upon the bending
loads (3.6 nm for Pr, 6.2 nm for Et and 10.1 nm for Me). Figure 6:  The list
of representative snapshots from the simulations of (a) Pr-brittle and (b) Et-
elastic deformations.

###  III.2 Atomistic motions upon the deformation

While analyzing their dynamic trajectories, we observed that molecules rotate
strongly upon bending. Fig. 1 defines the alignments ($\alpha,\beta,\gamma$)
for each molecule that can rotate along the $x,y,z$ axes in the Cartesian
coordinates. The distributions of molecular rotations under bending are shown
in Fig. 5. Given that indentation direction acts on the $z$-axis and the setup
of three bending points aligns along the $x$-axis, we expect that the rotation
along $y$ axis ($\beta$) is the primary motion under the loading. Indeed, Fig.
5 reveals that the rotation in $\beta$ is more pronounced that other
directions.

Figure 7:  The list of representative snapshots from the simulation of Me-
plastic deformation.

To understand the role of molecular rotation in the whole deformation process,
we plotted a few representative structures from the MD trajectory for each
system in Fig. 6 and Fig. 7 and analyzed their patterns as follows.

Pr-brittle. Upon deformation, we found that the sample continuously to bend
from 0 to 2.5 nm (the first row of Fig. 6a) and 3.5 nm (the second row of Fig.
6a). The Pr molecules barely rotate around $x$ and $z$ axis. However, the
rotation on $y$-axis is more pronounced and it symmetrically distributed
around the central indenter. When the indentation depth exceeds 4.2 nm (the
last row of Fig. 6a), the lower surface cracks due to a large tensile stress.

Et-elastic. Upto the indentation depth of 4.0 nm (the first row of Fig. 6b),
the Et molecules barely rotate around the $x$ and $z$ axis, while the rotation
on $y$-axes ($\beta$) is more pronounced and it symmetrically distributed
around the central indenter. However, it is clear that the molecules around
the center of $y$-axis do not rotate. Upon further indentation at 5.0 nm (the
second row of Fig. 6b) and 6.2 nm (the last row of Fig. 6b), the molecules at
the center of lower surface undergo a large rotation around the $x$ and $z$
due to a large tensile stress, but do not rotate around $y$. This suggests
that molecules upon tension prefer a rotation on $\alpha$ and $\gamma$, rather
than the rotation around $\beta$ due to the anisotropic behavior of its
potential energy landscape. Since the rotations are symmetrically distributed
around the indenter, it is an elastic deformation. When the indentation is
released, the process is supposed to be reversible. Me-Elastic sample
undergoes very similar processes except that the critical indetentation depth
(10.2 nm) is larger.

Me-plastic. At 5.5 nm, we found that the Me molecules near the indenter (first
row of Fig. 7) have alternative changes of $\alpha$ and $\gamma$ angles, which
is similar to that in Fig. 6b. In addition, these molecules have non-symmetric
distribution of $\beta$ angles, which signals a phase transition triggered by
the large compressive stress in the upper surface due to bending. This domain
of new phases, consisting of realigned molecules (denoted as the red dotted
eclipse), can easily slip along its interface with the parent domain. Upon
indentation, the molecules in the secondary domain do not gain enough momentum
to go downward as compared to other molecules. Therefore, the relative
slipping direction of the secondary domain is upward and we observed the
appearance of a bump near the indenter tip (second row of Fig. 7) at 6.7 nm.
As the tip continues to go down, the secondary domain keeps climbing up until
the bump reaches its maximum. In the mean time, the molecules at the center
bottom region are nearly flattened, which can trigger another phase transition
to form a new phase domain. Upon further compression, the flattened molecules
at the center bottom region create much empty space along the $z$-axis. Thus,
the secondary domain slips down to push the neighboring molecules down to fill
the empty space (third row of Fig. 7) at 9.5 nm. When the indentation is
released, the process is supposed to be irreversible at low temperature since
triggering the back transformation requires some energy barrier. Therefore, it
is a plastic deformation.

Clearly, such a plastic deformation process is driven by the molecular
rotation, which is different from metal’s plastic deformation that requires
the migration of dislocations Li  _et al._ (2018); Zhang  _et al._ (2014);
Nöhring  _et al._ (2016); Zhuo and Beom (2018); Katakam and Yedla (2020); He
and Ma (2022). In several recent experimental studies, it has been proposed
that molecular rotation may play a central role to generate a crystal twining
Takamizawa  _et al._ (2018) or phase transition Takamizawa and Miyamoto
(2014); Karothu  _et al._ (2016) which leads to plastic deformability. Our
simulation on Me-plastic revealed a similar atomistic picture except that its
new domain size is much smaller. Due to molecular rotation, some Me molecules
near the indenter form a new phase. The newly formed secondary phase can
freely slide along the interface to adjust the local stress. In the early
stage, the upward movement of re-aligned molecules results in a bump shape
near the indenter (instead of two bumps being symmetrically aligned near the
indenter). Such an asymmetric bump has actually been found in the bending
experiment Devarapalli  _et al._ (2019), which may provide another evidence to
support our modeling results. Given that most of the previous bending
experiments did not report the finding of new domains, it is likely that only
very small domains of rotated molecules can be formed due to energetic reasons
under the plastic bending deformation. In this case, a reliable atomstic
modelling is needed to capture such subtle details. Furthermore, if the
temperature is sufficiently high to cross the phase transition barrier, the
process may become reversible, similar to the previously reported superelastic
or shape-memory phenomenon Takamizawa and Miyamoto (2014); Karothu  _et al._
(2016); Takamizawa  _et al._ (2018).

###  III.3 Rotation-dependent Energy Map

So far, we have established the relation between molecular rotation and the
observed mechanical bending flexibility of organic crystals. However, we are
still unclear why some materials are more compliant than others and why we
observed two different deformation behaviors on the Me crystal with slightly
different initial configurations. To quest their physical origins, it is
necessary to examine the potential energy surface (PES) with respect to the
molecular rotations.

Figure 8:  The computed rotation-dependent energy map for Me-plastic from the
GAFF model and its physical interpretation. Figure 9:  The potential energy
surface from the GAFF model as a function of molecular rotation for three
crystals with different mechanical responses: (a) Pr-brittle, (b) Et-elastic,
(c) Me-elastic, and (d) Me-plastic deformations. The dashed lines in each
subplot denote direction of symmetric bending. The while region in (a) denotes
the rotations leading to energy exceeding 10000 kcal/mol.

To compute the rotation-dependent energy map, we started with the perfect
crystal structures and tracked the energy changes while systematically
rotating two groups of symmetrically-related molecules (colored in red and
blue in Fig. 2) along the $y$-axis in the unit cell. Using Me-plastic as an
example, we computed its energy map as the function of the rotation angles
($R_{1}$ and $R_{2}$) as shown in Fig. 8. In this map, it consists of two main
energy basins. The basin around (0, 0) represent the global minima (GM)
configuration around the equilibrium state. In the Me crystal, the molecules
are aligned with nonzero inclination angles. The arrow from (0, 0) to (-30,
30) represents a bending of two molecules along the opposite directions
(namely, clockwise and anti-clockwise directions). On the other hand, there
exists another local minimum (LM) of energy at (30, 30), which represents a
flattened configuration with both molecules being aligned horizontally. Such a
state can be achieved through barrier crossing by adding a large tensile or
compressive strains as shown in the arrow from (0, 0) to (30, 30).

Consequently, we applied this approach to compute the rotational-dependent
energy maps for all systems and attempted to find a predictive model to link
the possible deformation mechanism with our atomistic simulations. The results
are summarized as follows,

  * •

Pr has a very stiff GM (see Fig. 9a). This indicates that even a slight
rotation can lead to a high energy penalty. The energy basin of GM is aligned
diagonally. In this energy basin, the total energy increases over 1000
kcal/mol if two molecules bend symmetrically from (0, 0) to ($\pm 10,\mp 10$).
Such a high energy penalty would eventually lead to the formation of crack. In
addition, there is a LM centered around (20, 20). But this state is nearly
inaccessible from the GM due to a high energy barrier. Hence, Pr has a limited
rotational freedom, which is consistent with its brittle nature.

  * •

Et has more spreads around the GM (Fig. 9b). As shown in Fig. 5, two peaks are
symmetrically distributed at $\pm$20 degrees when the system reaches the
elastic limit. The rotation from (0, 0) to ($\pm 20,\mp 20$) would lead to a
penalty energy of 500 kJ/mol. Therefore, the Et molecules can rotate more than
Pr before the crack event starts. Similarly, Et has a LM around (30, 30) with
a high energy barrier.

  * •

Me-elastic (see Fig. 9c) has a shape similar to Et (Fig. 9b), except that it
has a wider bending region as denoted by the dotted line. Similarly, it has a
high energy barrier that prevents the phase transition to the adjacent LM
through the flattening motion. Therefore, Me molecules can bend more easily
than Et, but they cannot reach the LM state due to a high barrier.

  * •

Me-plastic has the flattest GM basin (Fig. 9d). The energy barrier of
symmetric bending from (0, 0) to ($\pm 30,\mp 30$) is only about 500 kcal/mol.
Hence the Me molecules can bend more easily than Et. More interestingly, there
is a low energy pathway that connects the LM at (30, 30) to the GM basin.
Under the bending deformation, the molecules in a large non-periodic supercell
may access other states due to the thermal fluctuation. The required barrier
crossing from GM to LM can be further reduced due to the surface molecules,
strain and other defects. Hence, it is possible to trigger the nucleation of a
secondary domain with re-aligned molecules in the LM state. According to the
nature of bending, such phase transition is more likely occur in either the
upper or lower surface due to extra tensile/compressive strains. And the
reoriented molecules (near the LM state) result in a stronger peak around
$\beta$ = 30∘ as compared to that around -30∘ in Fig. 5.

From the above analysis, it is clear that each type of deformation has its own
characteristics in their rotation-dependent energy maps despite that the model
is restricted to single unit cell assumption. First, a brittle deformation
should correspond to a stiff GM with high energy penalty to bend. When the GM
becomes less stiffer, the system tends to have more elastic region and becomes
more compliant. Finally, the key to achieve a plastic deformation is to have a
low energy barrier between the GM and its adjacent LM states. Clearly, such a
simplified energy map is instructive to understand the trend of bending
deformation when molecular rotation is the major factor. In addition, we
checked the rotation-dependent energy maps with two other energy models
(OpenFF and DFTB-TS). Encouragingly, the results are overall very similar (see
extended analysis and discussion in the Appendix A). Hence, we may be able to
employ this model to predict the deformation behaviors for new organic
crystals without performing expensive large scale MD simulations in the future
work.

##  IV Concluding Remarks

In this work, we have performed extensive molecular dynamics simulations to
directly model the mechanical bending of organic crystals. Using three
recently reported naphthalene diimide derivatives as the examples, our
simulation successfully produced different deformation behaviors from brittle
fracture to elastic/plastic deformation upon mechanical bending. By analyzing
the atomistic trajectories from our simulations, we discovered that molecular
rotational freedom is the key factor determining a material’s bendability,
which arises from the delicate interplay between geometric packing and
intermolecular interactions. Furthermore, we found the rotation-dependent
potential energy surface can be used to clarify the origin of different
mechanical deformation for organic materials. Although the role of molecular
rotation in driving the plastic bending have been recognized in several
previous experiments based on the observation of twin formation and phase
transitions Takamizawa and Miyamoto (2014); Karothu  _et al._ (2016);
Takamizawa  _et al._ (2018), our work extends this mechanism to more general
cases in which the rotated molecules do not necessarily form a large domain to
allow the plastic deformability.

While we focused on only three naphthalene diimide derivatives in this study,
the proposed three-point bending setup is entirely general to handle different
organic systems as long as the crystal structures and orientations are known.
In future, we will continue to test this approach on other crystallographic
directions Devarapalli  _et al._ (2019) and other systems Reddy  _et al._
(2005b); Raju  _et al._ (2018); Zhang  _et al._ (2021); Panda  _et al._
(2015); Reddy  _et al._ (2005a); Saha  _et al._ (2018). Additionally, the
impacts of model size, strain rate and indenter shape on other systems need to
be studied to ensure the simulation pipeline is transferable to other systems.

In parallel to this work, we recently proposed a crystal packing similarity
model Zhu  _et al._ (2022) that can rapidly identify the organic crystals with
similar packing and intermolecular interaction. Combining it with the present
atomistic modelling approach, we hope to develop a full simulation pipeline to
screen new mechanically flexible organic crystals from the available
databaseTaylor and Wood (2019) for future device applications.

###### Acknowledgements.

This research is sponsored by the NSF (DMR-2142570) and Sony Group
Corporation. The computing resources are provided by ACCESS (TG-DMR180040).
The authors also thank Changquan Sun, Reddy Malla, Pance Naumov, Liang Li and
Sinisa Mesarovic for helpful discussions.

## Code availability

The codes used to calculate the results of this study are available in
https://github.com/MaterSim/OST.

##  Appendix A Validation with Other Energy Models

As discussed in the main text, the GAFF model, as well as other energy models,
can describe the equilibrium lattice constants reasonably well. Since our
simulations also involve samples with large deformation, it is necessary to
validate the feasibility of GAFF in describing the configurations with large
deformation. Hence, we performed additional validations from the following
aspects.

First, we repeated the same bending simulations with the OpenFF model and the
results are qualitatively similar to the simulations based on GAFF. Namely, we
observed the same behaviors of brittle fracture, elastic and plastic
deformations for the three systems. Given that the GAFF and OpenFF models are
parameterized from completely different protocols, the agreement from two
independent FF parameters indicates the observed phenomena should be general
and invariant with the choice of force field models.

Figure 10:  The comparison of rotation-dependent energy maps from different
models.

Second, we proposed the idea of rotation-dependent energy maps to understand
the atomistic mechanism of bending in the main text. Fig. 10 displays the
comparison of rotation-dependent energy maps from different models, including
GAFF (upper panels), OpenFF (middle panels), and DFTB-TS (lower panels). It
can be clearly seen that the GAFF results are remarkably similar to the OpenFF
results for all three systems. Due to the convergence issue in self-consistent
field method of DFTB, we omitted the configurations with rotations smaller
than -15∘. The DFTB-TS approach, as a more first-principle model, also yields
consistent GM and LM shapes for each system, despite that the overall energy
surfaces are much smoother. If we compare the GM-LM transitions, the elastic
system generally requires higher energy barrier as compared to the plastic
systems for all three energy models, thus confirming our interpretation that
the nucleation of LM in the plastic system is easier due to a smaller energy
barrier. There is only one notable difference that DFTB-TS predict that the
GM-LM transition in the brittle system requires comparable barrier than that
in the plastic system. However, such transition should be prevented by the
interlocking molecular packing. Hence, it does not impact our main conclusion
that the transition between LM and GM in the plastic system is most favorable.

## References

  * Naumov  _et al._ (2015) Panče Naumov, Stanislav Chizhik, Manas K Panda, Naba K Nath, and Elena Boldyreva, “Mechanically responsive molecular crystals,” Chem. Rev. 115, 12440–12490 (2015).
  * Saha  _et al._ (2018) Subhankar Saha, Manish Kumar Mishra, C Malla Reddy, and Gautam R Desiraju, “From molecules to interactions to crystal engineering: mechanical properties of organic solids,” Acc. Chem. Res. 51, 2957–2967 (2018).
  * Reddy  _et al._ (2006) C Malla Reddy, K Anantha Padmanabhan, and Gautam R Desiraju, “Structure-property correlations in bending and brittle organic crystals,” Cryst. Growth Des. 6, 2720–2731 (2006).
  * Takamizawa and Miyamoto (2014) Satoshi Takamizawa and Yasuhiro Miyamoto, “Superelastic organic crystals,” Angew. Chem. Int. Ed. 53, 6970–6973 (2014).
  * Ghosh and Reddy (2012) Soumyajit Ghosh and C Malla Reddy, “Elastic and bendable caffeine cocrystals: implications for the design of flexible organic materials,” Angew. Chem. Int. Ed. 51, 10319–10323 (2012).
  * Panda  _et al._ (2015) Manas K Panda, Soumyajit Ghosh, Nobuhiro Yasuda, Taro Moriwaki, Goutam Dev Mukherjee, C Malla Reddy, and Panče Naumov, “Spatially resolved analysis of short-range structure perturbations in a plastically bent molecular crystal,” Nat. Chem. 7, 65–72 (2015).
  * Krishna  _et al._ (2016) Gamidi Rama Krishna, Ramesh Devarapalli, Garima Lal, and C Malla Reddy, “Mechanically flexible organic crystals achieved by introducing weak interactions in structure: supramolecular shape synthons,” J. Am. Chem. Soc. 138, 13561–13567 (2016).
  * Yadav  _et al._ (2019) Jay Prakash Yadav, Ram Naresh Yadav, Piyush Uniyal, Hongbo Chen, Chenguang Wang, Changquan Calvin Sun, Navin Kumar, Arvind K Bansal, and Sanyog Jain, “Molecular interpretation of mechanical behavior in four basic crystal packing of isoniazid with homologous cocrystal formers,” Cryst. Growth Des. 20, 832–844 (2019).
  * Raju  _et al._ (2018) K Bal Raju, Subham Ranjan, VS Vishnu, Manjima Bhattacharya, Biswajit Bhattacharya, Anoop K Mukhopadhyay, and C Malla Reddy, “Rationalizing distinct mechanical properties of three polymorphs of a drug adduct by nanoindentation and energy frameworks analysis: role of slip layer topology and weak interactions,” Cryst. Growth Des. 18, 3927–3937 (2018).
  * Mishra and Sun (2020) Manish Kumar Mishra and Changquan Calvin Sun, “Conformation directed interaction anisotropy leading to distinct bending behaviors of two roy polymorphs,” Cryst. Growth Des. 20, 4764–4769 (2020).
  * Zhang  _et al._ (2021) Keke Zhang, Changquan Calvin Sun, Yu Liu, Chenguang Wang, Peng Shi, Jun Xu, Songgu Wu, and Junbo Gong, “Structural origins of elastic and 2d plastic flexibility of molecular crystals investigated with two polymorphs of conformationally rigid coumarin,” Chem. Mater. 33, 1053–1060 (2021).
  * Devarapalli  _et al._ (2019) Ramesh Devarapalli, Sourabh Bhagwan Kadambi, Chun-Teh Chen, Gamidi Rama Krishna, Bal Raju Kammari, Markus J Buehler, Upadrasta Ramamurty, and C Malla Reddy, “Remarkably distinct mechanical flexibility in three structurally similar semiconducting organic crystals studied by nanoindentation and molecular dynamics,” Chem. Mater. 31, 1391–1402 (2019).
  * Root  _et al._ (2017) Samuel E Root, Suchol Savagatrup, Adam D Printz, Daniel Rodriquez, and Darren J Lipomi, “Mechanical properties of organic semiconductors for stretchable, highly flexible, and mechanically robust electronics,” Chem. Rev. 117, 6467–6499 (2017).
  * Li  _et al._ (2019) Liang Li, Patrick Commins, Marieh B Al-Handawi, Durga Prasad Karothu, Jad Mahmoud Halabi, Stefan Schramm, James Weston, Rachid Rezgui, and Panče Naumov, “Martensitic organic crystals as soft actuators,” Chem. Sci. 10, 7327–7332 (2019).
  * Mutai  _et al._ (2021) Toshiki Mutai, Toshiyuki Sasaki, and Satoshi Takamizawa, “Organic soft crystals exhibiting spontaneously reversible mechano-responsive luminescence,” J. Photochem. Photobiol. C Photochem. Rev. , 100479 (2021).
  * Di  _et al._ (2022) Qi Di, Liang Li, Xiaodan Miao, Linfeng Lan, Xu Yu, Bin Liu, Yuanping Yi, Panče Naumov, and Hongyu Zhang, “Fluorescence-based thermal sensing with elastic organic crystals,” Nat. Commun. 13, 5280 (2022).
  * Hayashi and Koizumi (2016) Shotaro Hayashi and Toshio Koizumi, “Elastic organic crystals of a fluorescent $\pi$-conjugated molecule,” Angew. Chem. Int. Ed. 128, 2751–2754 (2016).
  * Wang  _et al._ (2019) Yu Wang, Lingjie Sun, Cong Wang, Fangxu Yang, Xiaochen Ren, Xiaotao Zhang, Huanli Dong, and Wenping Hu, “Organic crystalline materials in flexible electronics,” Chem. Soc. Rev. 48, 1492–1530 (2019).
  * Samanta  _et al._ (2023) Ranita Samanta, Susobhan Das, Saikat Mondal, Tamador Alkhidir, Sharmarke Mohamed, Satyaprasad P Senanayak, and C Malla Reddy, “Elastic organic semiconducting single crystals for durable all-flexible field-effect transistors: insights into the bending mechanism,” Chem. Sci. 14, 1363–1371 (2023).
  * Liu  _et al._ (2018) Huapeng Liu, Zhuoqun Lu, Zuolun Zhang, Yue Wang, and Hongyu Zhang, “Highly elastic organic crystals for flexible optical waveguides,” Angew. Chem. Int. Ed. 57, 8448–8452 (2018).
  * Sun (2017) Changquan Calvin Sun, “Microstructure of tablet—pharmaceutical significance, assessment, and engineering,” Pharm. Res. 34, 918–928 (2017).
  * Wang and Sun (2018) Chenguang Wang and Changquan Calvin Sun, “Identifying slip planes in organic polymorphs by combined energy framework calculations and topology analysis,” Cryst. Growth Des. 18, 1909–1916 (2018).
  * Wang and Sun (2019) Chenguang Wang and Changquan Calvin Sun, “Computational techniques for predicting mechanical properties of organic crystals: A systematic evaluation,” Mol. Pharmaceutics 16, 1732–1741 (2019).
  * Wang and Sun (2020) Chenguang Wang and Changquan Calvin Sun, “The landscape of mechanical properties of molecular crystals,” CrystEngComm 22, 1149–1153 (2020).
  * Ootani and Kubo (2022) Yusuke Ootani and Momoji Kubo, “Density-functional tight-binding molecular dynamics simulation of the bending mechanism of molecular crystals,” J. Phys. Chem. C 126, 10554–10565 (2022).
  * Matveychuk  _et al._ (2022) Yury V Matveychuk, Artem S Yurchenko, Artëm E Masunov, and Ekaterina V Bartashevich, “Quantitative comparison of structural and mechanical properties of 6-chloro-2, 4-dinitroaniline polymorphs,” Cryst. Growth Des. 22, 6472–6482 (2022).
  * Bryant  _et al._ (2018) MJ Bryant, AGP Maloney, and RA Sykes, “Predicting mechanical properties of crystalline materials through topological analysis,” CrystEngComm 20, 2698–2704 (2018).
  * Zhu  _et al._ (2012) Wenpeng Zhu, Hongtao Wang, and Wei Yang, “Orientation- and microstructure-dependent deformation in metal nanowires under bending,” Acta Mater. 60, 7112–7122 (2012).
  * Zhang  _et al._ (2014) J J Zhang, Y D Yan, X Liu, T Sun, and Y C Liang, “Influence of coherent twin boundaries on three-point bending of gold nanowires,” J. Phys. D: Appl. Phys. 47, 195301 (2014).
  * Nöhring  _et al._ (2016) Wolfram G. Nöhring, Johannes J. Möller, Zhuocheng Xie, and Erik Bitzek, “Wedge-shaped twins and pseudoelasticity in fcc metallic nanowires under bending,” Extreme Mech. Lett. 8, 140–150 (2016).
  * Zhuo and Beom (2018) X.R. Zhuo and H.G. Beom, “Atomistic study of the bending properties of silicon nanowires,” Comput. Mater. Sci. 152, 331–336 (2018).
  * Katakam and Yedla (2020) Krishna Chaitanya Katakam and Natraj Yedla, “Crack velocities and microstructural investigations in nickel nanowires with crack, crack-defect under mode-i and bending using large-scale molecular dynamics simulations,” Superlattice. Microst. 146, 106674 (2020).
  * He and Ma (2022) Yuxiao James He and Bin Ma, “Molecular dynamics analysis on bending mechanical behavior of alumina nanowires at different loading rates,” Trans. Nonferrous Met. Soc. China 32, 3687–3698 (2022).
  * Reddy  _et al._ (2010) C Malla Reddy, G Rama Krishna, and Soumyajit Ghosh, “Mechanical properties of molecular crystals—applications to crystal engineering,” CrystEngComm 12, 2296–2314 (2010).
  * Reddy  _et al._ (2005a) C Malla Reddy, Ravi C Gundakaram, Srinivas Basavoju, Michael T Kirchner, K Anantha Padmanabhan, and Gautam R Desiraju, “Structural basis for bending of organic crystals,” Chem. Commun. , 3945–3947 (2005a).
  * Mathew  _et al._ (2013) Nithin Mathew, Catalin R Picu, and Peter W Chung, “Peierls stress of dislocations in molecular crystal cyclotrimethylene trinitramine,” J. Phys. Chem. A 117, 5326–5334 (2013).
  * Mathew and Picu (2013) N Mathew and RC Picu, “Slip asymmetry in the molecular crystal cyclotrimethylenetrinitramine,” Chem. Phys. Lett. 582, 78–81 (2013).
  * Olson  _et al._ (2018) Isabel A Olson, Alexander G Shtukenberg, Bart Kahr, and Michael D Ward, “Dislocations in molecular crystals,” Rep. Prog. Phys. 81, 096501 (2018).
  * Plimpton (1995) Steve Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” J. Comput. Phys. 117, 1–19 (1995).
  * Schneider and Stoll (1978) T. Schneider and E. Stoll, “Molecular-dynamics study of a three-dimensional one-component model for distortive phase transitions,” Phys. Rev. B 17, 1302–1322 (1978).
  * Case  _et al._ (2020) D. A. Case, K. Belfon, I. Y. Ben-Shalom, S. R. Brozell, D. S. Cerutti, T. E. Cheatham, III, V. W. D. Cruzeiro, T. A. Darden, R. E. Duke, G. Giambasu, M. K. Gilson, H. Gohlke, A. W. Goetz, R. Harris, S. Izadi, S. A. Izmailov, K. Kasavajhala, A. Kovalenko, R. Krasny, T. Kurtzman, T. S. Lee, S. LeGrand, P. Li, C. Lin, J. Liu, T. Luchko, R. Luo, V. Man, K. M. Merz, Y. Miao, O. Mikhailovskii, G. Monard, H. Nguyen, A. Onufriev, F. Pan, S. Pantano, R. Qi, D. R. Roe, A. Roitberg, C. Sagui, S. Schott-Verdugo, J. Shen, C. L. Simmerling, N. R. Skrynnikov, J. Smith, J. Swails, R. C. Walker, J. Wang, L. Wilson, R. M. Wolf, X. Wu, Y. Xiong, Y. Xue, D. M. York, and P. A. Kollman,  _AMBER 2020_ (2020).
  * Wang  _et al._ (2004) Junmei Wang, Romain M Wolf, James W Caldwell, Peter A Kollman, and David A Case, “Development and testing of a general amber force field,” J. Comput. Chem. 25, 1157–1174 (2004).
  * Jakalian  _et al._ (2000) Araz Jakalian, Bruce L Bush, David B Jack, and Christopher I Bayly, “Fast, efficient generation of high-quality atomic charges. am1-bcc model: I. method,” J. Comput. Chem. 21, 132–146 (2000).
  * Wagner  _et al._ (2023) Jeff Wagner, Matt Thompson, David L. Mobley, John Chodera, Caitlin Bannan, Andrea Rizzi, trevorgokey, David Dotson, Jaime Rodríguez-Guerra, Josh A. Mitchell, Camila, Pavan Behara, Christopher Bayly, JoshHorton, Iván Pulido, Victoria Lim, Sukanya Sasmal, Lily Wang, SimonBoothroyd, Andrew Dalke, Daniel Smith, Josh Horton, Lee-Ping Wang, Ziyuan Zhao, and Yutong Zhao, “openforcefield/openff-toolkit: 0.12.0 Minor feature and bugfix release,”  (2023).
  * Boothroyd  _et al._ (2023) Simon Boothroyd, Pavan Kumar Behara, Owen Madin, David Hahn, Hyesu Jang, Vytautas Gapsys, Jeffrey Wagner, Joshua Horton, David Dotson, Matthew Thompson, Jessica Maat, Trevor Gokey, Lee-Ping Wang, Daniel Cole, Michael Gilson, John Chodera, Christopher Bayly, Michael Shirts, and David Mobley, “Development and benchmarking of open force field 2.0.0 — the sage small molecule force field,”  (2023).
  * Shirts  _et al._ (2017) Michael R Shirts, Christoph Klein, Jason M Swails, Jian Yin, Michael K Gilson, David L Mobley, David A Case, and Ellen D Zhong, “Lessons learned from comparing molecular dynamics engines on the SAMPL5 dataset,” J. Comput. Aided Mol. Des. 31, 147–161 (2017).
  * Hourahine  _et al._ (2020) Ben Hourahine, Bálint Aradi, Volker Blum, F Bonafé, A Buccheri, Cristopher Camacho, Caterina Cevallos, MY Deshaye, T Dumitrică, A Dominguez, _et al._ , “Dftb+, a software package for efficient approximate density functional theory based atomistic simulations,” J. Chem. Phys. 152, 124101 (2020).
  * Tkatchenko and Scheffler (2009) Alexandre Tkatchenko and Matthias Scheffler, “Accurate molecular van der waals interactions from ground-state electron density and free-atom reference data,” Phys. Rev. Lett. 102, 073005 (2009).
  * Li  _et al._ (2018) Qing-Jie Li, Bin Xu, Shotaro Hara, Ju Li, and Evan Ma, “Sample-size-dependent surface dislocation nucleation in nanoscale crystals,” Acta Mater. 145, 19–29 (2018).
  * Takamizawa  _et al._ (2018) Satoshi Takamizawa, Yuichi Takasaki, Toshiyuki Sasaki, and Noriaki Ozaki, “Superplasticity in an organic crystal,” Nat. Commun. 9, 3984 (2018).
  * Karothu  _et al._ (2016) Durga Prasad Karothu, James Weston, Israel Tilahun Desta, and Panče Naumov, “Shape-memory and self-healing effects in mechanosalient molecular crystals,” J. Am. Chem. Soc. 138, 13298–13306 (2016).
  * Reddy  _et al._ (2005b) C Malla Reddy, Srinivas Basavoju, and Gautam R Desiraju, “Sorting of polymorphs based on mechanical properties. trimorphs of 6-chloro-2, 4-dinitroaniline,” Chem. Commun. , 2439–2441 (2005b).
  * Zhu  _et al._ (2022) Qiang Zhu, Weilun Tang, and Shinnosuke Hattori, “Quantification of crystal packing similarity from spherical harmonic transform,” Cryst. Growth Des. 22, 7308–7316 (2022).
  * Taylor and Wood (2019) Robin Taylor and Peter A. Wood, “A million crystal structures: The whole is greater than the sum of its parts,” Chem. Rev. 119, 9427–9477 (2019).


* [41] T. L. van Noorden.  Crystal precipitation and dissolution in a porous medium: effective equations and numerical experiments.  Multiscale Modeling & Simulation, 7(3):1220–1236, 2009.
  * [42] T. L. van Noorden.  Crystal precipitation and dissolution in a thin strip.  European Journal of Applied Mathematics, 20(1):69–91, 2009.
  * [43] Lars von Wolff and Iuliu Sorin Pop.  Upscaling of a Cahn–Hilliard Navier–Stokes model with precipitation and dissolution in a thin strip.  Journal of Fluid Mechanics, 941:A49, 2022.
  * [44] Tianfu Xu, Javier Samper, Carlos Ayora, Marisol Manzano, and Emilio Custodio.  Modeling of non-isothermal multi-component reactive transport in field scale porous media flow systems.  Journal of Hydrology, 214(1):144–164, 1999.
  * [45] Zhijie Xu, Hai Huang, Xiaoyi Li, and Paul Meakin.  Phase field and level set methods for modeling solute precipitation and/or dissolution.  Computer Physics Communications, 183(1):15–19, 2012.


11institutetext: Escola de Ciências e Tecnologia
Universidade de Trás-os-Montes e Alto Douro, Vila Real 5000–811, Portugal
Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
#  Subtractive mountain clustering algorithm applied to a chatbot to assist
elderly people in medication intake

Neuza Claro and Paulo A. Salgado and T-P Azevedo Perdicoúlis

###### Abstract

Errors in medication intake among elderly people are very common. One of the
main causes for this is their loss of ability to retain information. The high
amount of medicine intake required by the advanced age is another limiting
factor. Thence, the design of an interactive aid system, preferably using
natural language, to help the older population with medication is in demand. A
chatbot based on a subtractive cluster algorithm, included in unsupervised
learned models, is the chosen solution since the processing of natural
languages is a necessary step in view to construct a chatbot able to answer
questions that older people may pose upon themselves concerning a particular
drug. In this work, the subtractive mountain clustering algorithm has been
adapted to the problem of natural languages processing. This algorithm version
allows for the association of a set of words into clusters. After finding the
centre of every cluster — the most relevant word, all the others are
aggregated according to a defined metric adapted to the language processing
realm. All the relevant stored information is processed, as well as the
questions, by the algorithm. The correct processing of the text enables the
chatbot to produce answers that relate to the posed queries. To validate the
method, we use the package insert of a drug as the available information and
formulate associated questions.

###### Keywords:

chatbot, medicine intake aid system, natural language processing, subtractive
mountain clustering.

##  1 Introduction

Physicians often prefer the treatment of illness using medication. One main
reason for this is the non-invasiveness of this method of cure and another is
the advance in science, namely in pharmacological engineering, that has given
rise to new drugs and techniques with great effectiveness.

Older adults are more susceptible to the use of medication because they are
also more prone to having chronic disorders such as high blood pressure,
cardiac arrhythmia and diabetes. The use of multiple drugs to treat cumulative
diseases — multiple pharmacy — and also the use of numerous medications to
treat a single condition — polypharmacy — are very common situations among
this age group. A study carried out in the period 2010–2011 in the United
States [1] showed that almost 90% of elderly adults regularly take at least
one prescription drug, approximately 80% take at least two, and 36% take at
least five prescribed medications.

With the upsurge in medication, errors associated with their use have also
raised, especially in older people [2]. The error associated with taking drugs
is especially problematic for these generations because, besides being the
ones with the most complex clinical conditions, they also have cognitive
problems related to memory and assimilation of information, making it
difficult to take the right medication at the right time and in the manner
prescribed by their doctor. It is estimated that close to half of the older
adults do not take their medication according to the doctor’s prescription
[3]. Incorrect drug administration can happen by (i) taking medication at a
different time of the day than it was prescribed; (ii) taking a different dose
from the prescribed one; and (iii) changing the medication in course on a
particular occasion [4].

These events may be caused by the similarity of the medicine box, the shape
and colour of the pill, the similarity between the names of the medicines and
the complexity and length of the medical prescription [4]. One study confirmed
that 25% of the medication errors are associated with confusion with the name
and 33% with confusion with the medication box and package insert [5].
Medication intake errors can lead to loss of treatment effectiveness and
increased risk of new complications may induce a new disease state, a new
hospitalisation or even death.

Technology information in the health area has been growing [6]. The use of
systems to avoid medication errors at the hospital and home medication at home
is vast. For example, Ahmed and coworkers [7] developed an automatic drug
dispensing system. Mobile phone applications are also numerous [8] [9].

This work aims to develop an automatic conversational agent, a chatbot,
capable of accompanying older people with taking medication. The person should
be able to interact with this aid system in natural language.

The increase in the processing capacity of computers and smartphones,
associated with the great development in artificial intelligence, has allowed
the development of largescale software to facilitate the daily lives of
humans. A bot, diminutive of a robot, is an automated hardware or software
machine with the capacity to simulate human behaviour [10]. Bots are powered
by advances in Artificial Intelligence (AI) technologies. Inside a bot, an
algorithm can produce a certain answer according to the input data. An example
of a bot is the chatbot, a program capable of having an online conversation
with a human being.

The evolution in machine learning algorithms, such as deep learning and deep
reinforcement learning, has improved natural language processing (NLP)
performance. These advances have made intelligent conversational systems gain
more and more popularity. Since then, chatbots have been applied in the most
diverse areas: from commercial use [11] to medicine [12] [13].

Various NLP tasks are carried out to resolve the ambiguity in speech and
language processing. Before machine learning techniques, all NLP tasks are
carried out using various rules-based approaches. In rule-based systems, rules
were constructed manually by linguistic experts or grammarians for particular
tasks. Machine learning and statistical techniques are everywhere in today’s
NLP. In literature, the implementation of various machine learning techniques
for various NLP tasks has been investigated extensively. The machine learning
systems start analysing the training data to build their knowledge and produce
their own rules and classifiers.

Machine learning techniques can be divided into three categories [14] [15]:

  * •

supervised learning whose models are trained using the labelled data set,
where the model learns about each type of data, that includes models like
hidden Markov model (HMM) [16] and support vector machines (SVM) [17].

  * •

semi-supervised learning that involves a small degree of supervision, and one
example is bootstrapping [18].

  * •

unsupervised learning whose model is not trained. The most common approach of
the unsupervised category is clustering [19] [20].

Deep learning is a subfield of machine learning based on artificial neural
networks which try to learn from the layered model of inputs. In the deep
learning approach concept, learning of a current layer is dependent on the
previous layer input. Deep learning algorithms can fall into both supervised
and unsupervised categories [21]. The main applications of deep learning
include pattern recognition and statistical classification. For example, to
combat the increasing amount and reduce the threat of malicious programs,
novel deep learning was developed [22].

In this work, we develop an aid system for medication intake — a chatbot —
that gives answers to queries related to the medicine intake. Previously, the
system had to been informed about the user’s prescription, that is, the daily
medicine intake routine and detailed information about every drug. We use the
subtractive mountain clustering (SMC) algorithm to perform NLP tasks.

The main objective of this work is to apply the SMC algorithm within the
chatbot to facilitate communication. Namely to process natural language using
the SMC algorithm.

The topics covered in this paper relate mainly to chatbots and NLP. To
conclude this section, we start by recalling some important facts in the
history of chatbots in Subsection 1.1 and in Subsection 1.2 we review some
existing chatbot solutions in healthcare. In Section 2, the chatbot
configuration is described. Then, in Section 3, an overview of NLP and its
implementation is discussed. In Subsection 3.1, we present the required steps
for text pre-processing and, in Subsection 3.2, the SMC algorithm, necessary
to build the response, is described. Finally, in Section 4, we present and
discuss the results obtained for the case study. We conclude with Section 5
where the main outcomes of the work are outlined and some guidelines for
future work are stated.

###  1.1 History of chatbots

The first chatbot was created in 1966 by Joseph Weizenbaum. It was called
ELIZA [23], and the objective was to pretend to be a psychologist. To do this,
it used simple rules of conversation and rephrased most of what the users said
to simulate a Rogerian therapist — person-centreed therapy. In 1991, the
Loebner Prize, an annual competition in artificial intelligence, was launched.
The contest awards the computer programs considered to be the most human-like.
The competition takes the format of a standard Turing test, i.e., in each
round, a human judge simultaneously holds textual conversations with a
computer program and a human being via computer and, based upon the responses,
the judge must decide which one is which [24]. In 2014, a chatbot named Eugene
Goostman managed to fool 33% of the judges, thereby beating the test.

Another example of a chatbot is using Natural Language Interface to Database
(NLIDB) to access information in the databases instead of Structured Query
Language (SQL). An NLIDB system is proposed as a solution to the problem of
accessing data in a simple way: any user can access the information contained
in the database and get the answer in natural language [25]. Nowadays,
multiple virtual assistants already exist, being the most complex and more
widely used. Siri (from Apple), Google Assistant (from Google) and Alexa (from
Amazon). At the moment, they are mainly used to call people, ask for
directions or search the Internet for information [26].

###  1.2 Chatbots in healthcare

These days, chatbots also have increased use in healthcare to treat disorders
such as cancer [27] or induce behavioural changes such as quitting smoking
[28] and weight control [29]. Chatbots are increasingly being adopted to
facilitate access to information from the patient side and reduce the load on
the clinician side. In the field of medicine, there are already chatbots for
the most varied purposes. Next, some illustrative examples are reviewed:

  1. 1.

OneRemission is a healthcare chatbot to help cancer survivors, fighters, and
supporters to learn about cancer and post-cancer health care [30].

  2. 2.

Wysa is an emotionally intelligent chatbot. Its purpose is to help the user to
build mental resilience and promote mental well-being with a test-based
interface [31].

  3. 3.

Florence is a chatbot related to medication intake that can remember taking
medication, monitor certain biomedical parameters, and find information about
diseases [32].

After reviewing the chatbots already developed in healthcare, we did not find
one that focuses on taking medication by older adults, like the chatbot
presented here.

##  2 Chatbot configuration

The chatbot herein presented was designed to help older adults with their
medication. The chatbot can provide information about the physician’s
prescription, the medicine’s package insert, and also extra information, such
as the colour of the box and the colour of the pill, to avoid confusion in
taking medication. Regarding the medical prescription, the chatbot can inform
about the dose (how many pills per day) and when to take it (at which part of
the day). Another set of information relates to the medicine’s package insert,
for example, indications, side effects, and what to do in case of forgetting
to take it. Finally, the chatbot can also provide an image of the medicine box
and the pill so that the patient can easily identify the medicine to be taken.
Once the medicine in question is identified, then the related information can
be retrieved.

The set of questions and answers must be as similar as possible to a
conversation between human beings. Chatbots are developed to connect with the
users and feel that they are communicating with a human and not a bot. In our
chatbot, the possible answers are predefined and designed to emulate daily
human communication.

Since a main step in the construction of the chatbot is NLP, this topic will
be discussed in Section 3.

##  3 Natural language processing

NLP is an area of computer science, more specifically, in artificial
intelligence (AI), concerned with giving computers the human ability to
understand text and spoken words. This processing generally involves
translating natural language into data (numbers) that a computer can use and
generate a certain answer. NLP is applied in several areas, such as, machine
translation [33], text summarisation [34] and spam detection [35]. One of the
NLP’s uses is chatbots.

A chatbot system analyses a query posed by a person and generates an answer
from an organised collection of data stored and accessed electronically from a
computer system. Usually, the answer is retrieved based on the basic keyword
matching, and a selected response is then given as the output. When we talk
about natural language, there are many ways to say the same information. So,
when the chatbot is faced with several alternative sentences requesting the
same information, it is necessary to use an algorithm that can import what is
truly relevant. After selecting this information, the chatbot will be able to
insert the phrase into a context to produce the proper answer. Furthermore,
the data needs to be pre-processed so that the chatbot can easily understand
it.

In Figure 1, one can see a flowchart that shows the process necessary for the
chatbot to generate an answer fitting the user input text. In the following
subsection, we will detail the pre-processing of text.

We apply word processing to both (i) the drug’s package insert, which is used
to define the algorithm that leads to the answers, and (ii) the user’s
queries.

Figure 1: Flowchart of the process needed from input to response in the
chatbot.

###  3.1 Text pre-processing

Natural Language Understanding (NLU) converts natural language utterances into
a structure that the computer can deal with. As a sentence cannot be processed
directly to the model, it needs to perform some NLP to further operations.

To achieve this, we implement an algorithm in MATLAB that uses a set of
functions included in the Text Analytics Toolbox. Hence:

  1. 1.

Read the file with raw text.

  2. 2.

Segment the text into paragraphs — documents.

  3. 3.

Remove empty documents.

  4. 4.

Segment the text into tokens — tokenisation. This process consists in
splitting the text into tokens, which are the basic units.

  5. 5.

Remove unwanted words (stop words) or irrelevant punctuation from those
tokens. Stop words are the words present in a sentence that does not make much
difference if removed. For example, the words and, off, for.

  6. 6.

Use normalising techniques, which consist in finding the root/stem of a word.
For example, the words ”ends” and ”ending” are represented by the same stem -
”end”.

The various steps above are not always an easy task to carry out. Non-standard
words are often ambiguous. Some words can have different meanings depending on
the context. Moreover, there are also some acronyms and abbreviations that can
be misleading. For example, should an acronym be read as a word (IKEA) or
using each letter in a sequence (IBM)? The abbreviation ”Dr.” has a full stop
that can lead to the wrong separation of sentences.

The following step is word embedding, which consists in converting text to
numbers. Converting words into numbers will make the algorithm easier to apply
[36] [37] [38].

###  3.2 Building the response

Creating an algorithm that could relate the user’s queries to the answers that
contain the right information is a big deal for chatbots.

Before selecting the correct answer is necessary to identify to which
class/topic the query belongs. Here, the SMC algorithm is used to define word
sets (clusters) [39].

Based on the medication package’s insert is defined a proximity degree between
the words related to medication intake. The calculus of the proximity degree
is based on the distance between any two words of the package insert. Only
relevant words are taken into account to calculate the distance. The word
frequency and size are taken into account to evaluate relevance. One possible
criterion is, for instance, a word appearing more than twice and having more
than two letters being classified as relevant.

The calculus of the distance between two relevant words is based on their
relative position in the text, where $r$ and $s$ are the codes of the words
that occupy positions $i$ and $j$, respectively. Hence, we distinguish three
different cases:

  * •

Two words are in the same sentence: the distance reflects the number of words
that separate them. The end of sentence recognition is done by punctuation
marks. Then, $D(r,s)=\left|j-i\right|$.

  * •

Two words are in the same document but in different sentences: the distance
reflects the number of sentences that separate them and the number of words.
Let $S_{n}$ be the number of separation sentences, then their pair distance is
given by $D(r,s)=\left|j-i\right|\textperiodcentered S_{n}\textperiodcentered
a$, where $a$ is an adjustable parameter according to the average number of
words in each sentence.

  * •

Two words are in different documents: the distance only reflects the number of
documents that separate them. Let $P_{n}$ be the number of separation
paragraphs, then their pair distance is given by
$D(r,s)=P_{n}\textperiodcentered b$, where $b$ is an adjustable parameter
according to the number of words in sentences and sentences in the paragraph.

This is indeed a distance, since is non-negative, symmetric and verify the
triangular inequality.

As the words may appear many times, the minimum distance between them is
considered.

In addition, we also need a factor $B(r,s)$ for each word pair that reflects
the number of times each pair appears together in the same sentence and/or the
same paragraph.

In the SMC algorithm, we use the distance and $B(r,s)$ to calculate the
potential. The algorithm is presented next.

####  3.2.1 Subtractive mountain cluster algorithm

Sometimes, it is necessary to reduce the size of the data set to a set of
representative points. For example, fuzzy logic algorithms are very complex
and are not applicable in a large data set [40].

The SMC approach, developed by Chiu [39], assumes that enormous data sets are
partitioned into subsets called clusters, and each cluster is represented by
one representative element called the cluster centre. Initially, all data set
points are potential cluster centres and each point potential is calculated by
equation (1). Hence

$P(I)=\sum_{j=1}^{N}e^{-\alpha{\left|x_{i}-x_{j}\right|^{2}}},$ (1)

for i = 1, …, N, ($N$ is the number of points) and $\alpha=4/(r_{a})^{2}$ for
constant $r_{a}>0$. Equation (1) shows that a data point with many
neighbouring data points will have a high potential value. Parameter $r_{a}$
defines the data points influence. Data points outside this radius have little
influence on potential.

The SCM algorithm allows for the association of a set of words into clusters.
After finding the centre of every group — the most relevant term, all the
others are aggregated according to a defined metric adapted to the language
processing realm.

We adapted equation (1) to NLP and define equation (2). So the greater the
number of times a word pair appears, the greater the potential of that word
pair. Hence

$P(r,s)=B(r,s)\sum_{s=1}^{N}e^{-\alpha{{D(r,s)}^{2}}},$ (2)

where $D$ is the symmetric distance matrix, $N$ the number of words. We assume
the potential of equals words pairs is approximately zero.

After the potential of every distance has been computed, as the matrix $D$ is
symmetric, we sum up all the potentials for each and every word, and we
obtained a vector with $N$ columns. Thus, we were able to obtain the potential
for every word and choose the one with the greatest potential as the centre of
the first cluster. Let $I_{1}^{*}$ be the first word cluster centre index and
$P_{1}^{*}$ its potential value.

Equation (3) is an adaptation of the potential equation for new cluster centre
words after determining the first cluster centre word. Hence

$P^{*}(r,s)=P(r,s)-P_{1}^{*}B(r,s)e^{-\beta{{D({I_{1}^{*}},s)}^{2}}},$ (3)

where $\beta=4/(r_{b})^{2}$ for constant $r_{b}>0$. With Equation (3), we
subtract a portion of potential at each word pair according to the distance
from each word to the word chosen as the centre of the first cluster. Words
close to the first cluster centre will have significantly reduced potential
and are unlikely to be selected as the next cluster centre. Parameter $r_{b}$
defines the radius affected to the potential reduction.

Now, the word with more potential is selected as the second cluster centre. We
then further reduce the potential of each word according to their distance to
the second cluster centre. In general, each time we select the next centre of
the next cluster, we revise the value of the potential in the following
manner:

$P_{k}^{*}(r,s)=P^{*}(r,s)-P_{k-1}^{*}B(r,s)e^{-\beta{{{D({I_{k}^{*}},s)}^{2}}}}.$
(4)

This iterative process ends when the word with the most potential $P_{k}^{*}$
is less than $\epsilon P_{1}^{*}$ where $\epsilon$ is a small fraction.

The minimum distance $d_{min}$ between every two cluster centres also needs to
be defined. If the following inequality is verified, the point is accepted as
a possible centre of a cluster:

$\frac{d_{min}}{r_{a}}+\frac{P_{k}^{*}}{P_{1}^{*}}\geq 1.$

Next, the belonging degree of each word to the cluster is calculated using
equation (5):

$U(I)=\frac{1}{\sum_{k=1}^{c}{\left(\frac{\parallel
x_{i}-c_{j}\parallel}{\parallel
x_{i}-c_{k}\parallel}\right)}^{\frac{2}{m-1}}},$ (5)

where $m$ is the hyper-parameter that controls how fuzzy the cluster is.

After applying the SMC algorithm, we can get groups of words (clusters)
according to the distance they appear in the medicine package insert. It is
necessary to pre-process the text of the questions, extracting the relevant
information. Therefore, we can attribute them to a cluster and generate an
appropriate answer.

##  4 Results and discussion

The parameters used are available in Table 1. We used the Xarelto drug package
insert to validate the algorithm, which is used to prevent blood clots.

Table 1: Parameters Parameters | Value
---|---
$a$ | 10
$b$ | 20
$r_{a}$ | 12
$r_{b}$ | 14
$\epsilon$ | 0.1
$m$ | 2

We obtained 297 documents (paragraphs) and 838 tokens (include words and two
important punctuation marks).

Figure 2 shows the number of times the first 40 most frequent words appear.
The most frequent word is the name of the drug — ”Xarelto” and then ”patient”.
After the first ten words, the frequency of each expression remains
approximately the same.

Figure 2: Number of times the first 40 most frequent words appear.

After applying the word relevance metrics, we get 468 relevant terms
representing 55.85% of the total words.

The distance between relevant words is represented in Figure 3. On each axis,
we have the indices of the words. The closer to white the colour, the smaller
the distance between words.

Figure 3: Representative diagram of word distance. The white colour represents
the closing words.

We get 20 clusters and the first five are represented in Figure 4. In the
graphs, the degree of belonging is also represented. We can see that in each
cluster there are at least 14 words with a degree of belonging equal to 1.

Figure 4: First five word clusters in function of their normalising belonging
degree.

In addition, we calculate the degree of belonging of each word to the cluster
and in Figure 5 is the graph that shows the number of words belonging to each
cluster with a degree of belonging greater than 0.5.

Figure 5: The number of words belonging to each cluster with a degree of
belonging greater than 0.5

We can get the most significant words, with question pre-processing, and
associate them to a cluster or set of clusters. The answer is in accordance
with the belonging relationship of words to clusters.

For example, one question with the stemming words ”risk”, ”foetal” and ”bleed”
has multiple possible answers. We can sort answers by the relative relevance.
Answers correspond to the paragraphs (documents) whose profile of belonging to
the words in the question and the document are closest. The answer with 1.00
as relative relevance is ”Xarelto increases the risk of bleeding and can cause
serious or fatal bleeding. In deciding whether to prescribe Xarelto to
patients at increased risk of bleeding the risk of thrombotic events should be
weighed against the risk of bleeding.”

##  5 Conclusion and future work

Chatbots are present in different areas, including healthcare. In this field,
there are chatbots for many purposes, from psychologist chatbots to medication
reminders.

Chatbots have the potential to help elderly people in taking medication since
they can emulate human conversation. However, after a thorough literature
review, we conclude that they are still not widely used to solve and avoid
medication errors.

A fundamental step to enable communication between the users and chatbots is
the processing of natural languages. The NLP area involves a large scientific
community and contains a variety of associated algorithms.

We adapt the SMC algorithm to NLP. To the best of our knowledge, this
algorithm has not been used before to resolve this type of problem. This
algorithm does not define a priori the number of clusters, leaving rather more
freedom to the association. Instead, through word clusters creation based on
distance, we were able to define a degree of belonging among words and thus
build adequate answers to the posed queries.

In this work, an aid system to help elderly people with medicine intake is
proposed. Namely, a chatbot that is able to interact with the user in natural
language. To process natural languages, a fundamental step is to achieve a
communication system- -user. To do this, we apply the SCM algorithm to define
the relationship between words. Furthermore, we implement the proposed
solution and applied it to a simple problem, obtaining results that we
consider promising.

To complete this work, building a user interface might be the next step. The
refinement of the algorithm is also in place. Moreover, the study of the
accuracy of the SCM model should be done as well as to investigate whether
other algorithms can solve the same problem. In addition, a study of the
adherence of people of more advanced age to chatbots might be also an
interesting point of investigation.

## References

  * [1] E. D. Kantor, C. D. Rehm, J. S. Hass, A. T. Chan, and E. L. Giovannucci, “Trends in Prescription Drug Use Among Adults in the United States From 1999-2012,” _JAMA_ , vol. 314, no. 17, pp. 1818–1830, 2015\.
  * [2] S. Reddy. (2017, Jul.) Patients Make More Medication Mistakes. [Online]. Available: https://www.wsj.com/articles/patients-make-more-medication-mistakes-1500913414
  * [3] J. M. Ruscin and S. A. Linnebur. (2018, Dec.) Aging and Drugs. [Online]. Available: https://is.gd/ulrQx1
  * [4] J. J. Mira, S. Lorenzo, M. Guilabert, I. Navarro, and V. Pérez-Jover, “A systematic review of patient medication error on self-administering medication at home,” _Expert Opinion on Drug Safety_ , vol. 14, no. 6, pp. 815–838, 2015.
  * [5] A. Berman, “Reducing Medication Errors Through Naming, Labeling, and Packaging.” _Journal of medical systems_ , no. 28, pp. 9–29, 2004\.
  * [6] M. E. Reed, “The Health Information Technology Special Issue: New Real-World Evidence and Practical Lessons,” _The American Journal of Managed Care_ , vol. 25, no. 1, p. 12, 1 2019.
  * [7] A. Ahmed, M. R. R. A., and Barua, “Home Medication for Elderly Sick People,” in _2019 5th International Conference on Advances in Electrical Engineering (ICAEE)_ , 2019, pp. 369–374.
  * [8] M.-Y. Wang, P. Tsai, J. Liu, and J. K. Zao, “Wedjat: A Mobile Phone Based Medicine In-take Reminder and Monitor,” in _2009 Ninth IEEE International Conference on Bioinformatics and BioEngineering_ , 2009, pp. 423–430.
  * [9] J. Zao, M. Wang, P. Tsai, and J. Liu., “Smart phone based medicine in-take scheduler, reminder and monitor,” in _The 12th IEEE International Conference on e-Health Networking, Applications and Services_ , 2010, pp. 162–168.
  * [10] A. S. Gillis. (2020, Jan.) Bot (robot). [Online]. Available: https://whatis.techtarget.com/definition/bot-robot
  * [11] M.-H. Wen, “A conversational user interface for supporting individual and group decision-making in stock investment activities,” in _2018 IEEE International Conference on Applied System Invention (ICASI)_. IEEE, April 2018, pp. 216–219.
  * [12] T. Schmidlen, M. Schwartz, K. DiLoreto, H. L. Kirchner, and A. C. Sturm, “Patient assessment of chatbots for the scalable delivery of genetic counseling,” _Journal of Genetic Counseling_ , vol. 26, no. 6, pp. 1166–1177, 2019.
  * [13] E. Sandalova, J. G. Ledford, M. Baskaran, and S. Dijkstra, “Translational Medicine in the Era of Social media: A Survey of Scientific and Clinical Communities,” _Frontiers in Medicine_ , vol. 6, p. 152, 2019\.
  * [14] S. S. Dash, S. K. Nayak, and D. Mishra, “A Review on Machine Learning Algorithms,” in _Intelligent and Cloud Computing_ , vol. 153, 2021, pp. 495–507.
  * [15] W. Khan, A. D. J. Nasir, and T. Amjad, “A survey on machine learning models for Natural Language Processing (NLP),” _Kuwait Journal of Science_ , vol. 43, pp. 95–113, 10 2016.
  * [16] N. Ebrahimi, A. Trabelsi, S. Islam, A. Hamou-Lhadj, and K. Khanmohammadi, “An HMM-based approach for automatic detection and classification of duplicate bug reports,” _Information and Software Technology_ , vol. 113, pp. 98–109, 2019.
  * [17] M. S. Desai, S. P. Gururaj, and T. L. Prakash, “Analysis of Health Care Data Using Natural Language Processing,” in _2020 Second International Conference on Inventive Research in Computing Applications (ICIRCA)_ , 2020, pp. 242–246.
  * [18] S. Deepika and T. Geetha, “Pattern-based bootstrapping framework for biomedical relation extraction,” _Engineering Applications of Artificial Intelligence_ , vol. 99, p. 104130, 2021.
  * [19] A. Chokor, H. Naganathan, W. K. Chong, and M. Asmar, “Analyzing Arizona OSHA Injury Reports Using Unsupervised Machine Learning,” _Procedia Engineering_ , vol. 145, pp. 1588–1593, 2016.
  * [20] G. Liu, M. Boyd, M. Yu, S. Z. Halim, and N. Quddus, “Identifying causality and contributory factors of pipeline incidents by employing natural language processing and text mining techniques,” _Process Safety and Environmental Protection_ , vol. 152, pp. 37–46, 2021.
  * [21] S. Weidman, _Deep Learning From Scratch: Building With Python From First Principles_ , 1st ed. O’Reilly Media, Incorporated, 2019.
  * [22] J. Zhang, “CLEMENT: Machine Learning Methods for Malware Recognition Based on Semantic Behaviours,” in _2020 International Conference on Computer Information and Big Data Applications (CIBDA)_ , 2020, pp. 233–236.
  * [23] J. Weizenbaum, “ELIZA — a computer program for the study of natural language communication between man and machine,” _Communications of the ACM_ , vol. 9, no. 1, pp. 36–45, 1966.
  * [24] S. Janarthanam, _Hand-On Chatbots and Conversational UI Development_ , 1st ed. Packt Publishing Ltd., 2017.
  * [25] A. Kumar and K. Vaisla, “Natural Language Interface to Databases: Development Techniques,” _Elixir Comp. Sci. and Engg._ , vol. 58, pp. 14 724–14 727, 05 2013.
  * [26] L. Burbach, P. Halbach, N. Plettenberg, J. Nakayama, M. Ziefle, and A. C. Valdez, “”Hey, Siri”, ”Ok, Google”, ”Alexa”. Acceptance-Relevant Factors of V]irtual Voice-Assistants,” in _2019 IEEE International Professional Communication Conference (ProComm)_ , 2019, pp. 101–111.
  * [27] R. V. Belfin, A. J. Shobana, M. Manilal, A. A. Mathew, and B. Babu, “A Graph Based Chatbot for Cancer Patients,” in _2019 5th International Conference on Advanced Computing Communication Systems (ICACCS)_ , 2019, pp. 717–721.
  * [28] D. Calvaresi, J. Calbimonte, F. Dubosson, A. Najjar, and M. Schumacher, “Social Network Chatbots for Smoking Cessation Agent and Multi Agent Frameworks,” in _IEEE/WIC/ACM International Conference on Web Intelligence Web Intelligence_ , Thessaloniki, Grece, 10 2019, p. 7.
  * [29] C.Huang, M. Yang, C.Huang, Y. Chen, M. Wu, and K. Chen, “A Chatbot supported Smart Wireless Interactive Healthcare System for Weight Control and Health Promotion,” in _2018 IEEE International Conference on Industrial Engineering and Engineering Management (IEEM)_ , 12 2018, pp. 1791–1795.
  * [30] Keenethics. (2021, Jul.) OneRemission. [Online]. Available: https://keenethics.com/project-one-remission
  * [31] B. Inkster, S. Sarda, and V. Subramanian, “An Empathy-Driven, Conversational Artificial Intelligence Agent (Wysa) for Digital Mental Well-Being: Real-World Data Evaluation Mixed-Methods Study,” _JMIR Mhealth Uhealth_ , vol. 6, no. 11, p. e12106, Nov 2018.
  * [32] M. Bates, “Health Care Chatbots Are Here to Help,” _IEEE Pulse_ , vol. 10, pp. 12–14, 05 2019.
  * [33] Z. Zong and C. Hong, “On Application of Natural Language Processing in Machine Translation,” in _2018 3rd International Conference on Mechanical, Control and Computer Engineering (ICMCCE)_ , 2018, pp. 506–510.
  * [34] S. Jugran, A. Kumar, B. S. Tyagi, and V. Anand, “Extractive Automatic Text Summarization using SpaCy in Python NLP,” in _2021 International Conference on Advance Computing and Innovative Technologies in Engineering (ICACITE)_ , 2021, pp. 582–585.
  * [35] P. Garg and N. Girdhar, “A Systematic Review on Spam Filtering Techniques based on Natural Language Processing Framework,” in _2021 11th International Conference on Cloud Computing, Data Science Engineering (Confluence)_ , 2021, pp. 30–35.
  * [36] M. Virkar, V. Honmane, and S. U. Rao, “Humanizing the Chatbot with Semantics based Natural Language Generation,” in _2019 International Conference on Intelligent Computing and Control Systems (ICCS)_ , 2019, pp. 891–894.
  * [37] P. C. Gaigole, L. H. Patil, and P. M. Chaudhari, “Preprocessing Techniques in Text Categorization,” _IJCA Proceedings on National Conference on Innovative Paradigms in Engineering & Technology 2013_, vol. NCIPET 2013, no. 3, pp. 1–3, December 2013.
  * [38] D. Jurafsky and J. H. Martin, _Speech and Language Processing: An Introduction to Natural Language Processing, Computational Linguistics, and Speech Recognition_ , 1st ed. Prentice Hall PTR, 2000.
  * [39] S. Chiu, “Fuzzy Model Identification Based on Cluster Estimation,” _Journal of the Intelligent and Fuzzy Systems_ , vol. 2, pp. 267–278, 01 1994\.
  * [40] A. Chmielowiec, “Implementation of the mountain clustering method and comments on its practical use for determining clusters centres,” in _Materials, Technologies, Constructions - Constructions and design_. Oficyna Wydawnicza Politechniki Rzeszowskiej, 05 2019, ch. 3, pp. 45–56.

## Authors

Neuza Claro is a Biomedicine 3rd student at University of Trás-os-Montes e
Alto Douro. She enjoys programming and being in nature.

Paulo A. Salgado received the B.S. degree and the M.S. degree in electrical
and computer science from the Faculty of Engineering, University of Porto,
Portugal, in 1989 and 1993, respectively, and the Ph.D. degree in electrical
and computer science from UTAD University, Portugal, in 1999. He is currently
an Associate Professor with the Department of Engineering, School of Science
and Engineering of UTAD, Portugal. His research interests include artificial
intelligence, control, and robotics science.

T-P Azevedo Perdicoúlis graduated in Mathematics (Computer Science) at the
University of Coimbra in 1991, and pursued graduate studies at the University
of Salford, Manchester (UK), where she obtained an MSc degree in Electronic
Control Engineering in 1995 and a Ph.D. degree in Mathematics and Computer
Science in 2000. She works for the Universidade of Trás-os-Montes e Alto
Douro, Portugal, since 1991 where she is currently an Associate Professor. She
is a founding member of the research institute ISR, University of Coimbra, and
also of APCA, IFAC. Her main research interests are in the fields of
differential games, optimisation and simulation of gas networks as well as
identification theory. Besides research she has got an intense pedagogic
activity, having taught many different courses in applied mathematics to
engineers.