Local Structure of Atactic Polystyrene Investigated by Molecular

Dynamics Method

ANDREI V. KOMOLKIN, SERGEY G. POLUSHIN, VYACHESLAV B. ROGOZHIN,

ALEXANDRA A. LEZOVA, GALINA E. POLUSHINA, IRINA A. SILANTEVA

Department of Physics,

Saint Petersburg State University,

Saint Petersburg, 199034,

RUSSIAN FEDERATION

Abstract: Molecular dynamics computer simulation of three substances ethylbenzene (EB), pentastyrene (PS-5),

and polystyrene-25 (PS-25) was performed to investigate the local order of the phenyl rings in monomers and

side-chain polymers. Monomer molecules (EB) tend to be in T-configuration, which corresponds to isotropic

local structure. Phenyl rings in chained molecules PS-5 and PS-25 partly cooperate in both parallel-displaced and

“sandwich” configuration with π–πstacking. These configurations are locally anisotropic and lead to the increas-

ing of Kerr constant K. Analysis of the local structure was performed by calculating the cylindrical distribution

function.

Key-Words: molecular dynamics simulation, comb polymer, local structure, pi-pi conjugation, phenyl rings

Received: June 9, 2023. Revised: November 19, 2023. Accepted: December 22, 2023. Published: February 2, 2024.

1 Introduction

Isotropic liquids have local (short-range) orienta-

tional order. Local orientational order exists due to

the mutual orientation of densely packed anisometric

molecules. The correlation length characterizing this

order is small; it corresponds to two to three molec-

ular sizes. The situation is different in the isotropic

phase of liquid crystalline substances (LC). The cor-

relation length in them increases many times when ap-

proaching the temperature Tcof the phase transition to

the mesophase. This short-range orientational order

in an isotropic melt of mesogens appears in macro-

scopic effects, which are called pre-transition. De-

spite this name, the effects are detected by the electri-

cal birefringence method (Kerr effect) over a temper-

ature range tens of degrees wide.

The structure of polymers is different from the

structural organization of low molecular weight LC.

Polymer melts are formed by long flexible macro-

molecules in a coil conformation. For this reason, ob-

taining the LC state in polymers previously seemed

unlikely. The problem was solved when comb-

shaped polymers were chosen as the basis for creating

LC polymers. The roles of side groups were played

by liquid crystal molecules. At the initial stage of

their work, chemists attached mesogenic molecules

directly to the main polymer chain. Mesophase did

not arise in such polymers. Then the idea was born

to attach mesogenic molecules to the chain through

flexible spacers such as aliphatic chains. In this case,

the side groups acquire relative independence from

the chain and the ability to self-organize. This is

how liquid crystalline polymers appeared,[1].In-

side the medium of disordered polymer coils, there

is a subsystem formed by mutually ordered combs –

providers of orientational order.

Studies of the isotropic phase of comb-shaped LC

polymers have shown that their equilibrium electro-

optical properties are well described by the Landau–

de Gennes theory [2], [3]. In this, they are similar

to low molecular weight LC [4], [5]. This suggests

that the presence of polymer chains does not manifest

itself noticeably in the short-range order in the meso-

genic side group subsystem.

The results of these works drew our attention to

the need to study polymers with the direct addition

of anisotropic side groups to the chain. It is obvious

that in their melts the short-range order and the macro-

scopic effects associated with it will be different.

Atactic polystyrene, in which phenyl rings are

bonded to the main polymer chain, was chosen as a

model for such a comb polymer, despite this is not

LC. The phenyl rings make the main contribution to

short-range order and local anisotropy. The electro-

optical properties of melts of four fractions of styrene

and polystyrene oligomers were studied earlier [6]. It

turned out that for fractions with the highest degree

of polymerization, the Kerr constant Kdoes not de-

crease with increasing temperature, but, on the con-

trary, increases at temperatures above 120 ◦C. The ob-

served effect may be associated with the liquid–liquid

transition (ll-transition). It was observed in the equi-

librium and dynamic properties of polymers, includ-

ing polystyrene [7], [8]. The transition is explained

by the fact that upon heating and reaching the ll tran-

sition temperature, the cooperativity of the movement

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2024.23.9

Andrei V. Komolkin, Sergey G. Polushin,

Vyacheslav B. Rogozhin, Alexandra A. Lezova,

Galina E. Polushina, Irina A. Silanteva

E-ISSN: 2224-2678

Volume 23, 2024

of the monomer units of the chain decreases; instead

of segmental movement, the monomer unit becomes

the main kinetic unit.

The increase in Kmay be due to the mutual or-

dering of the optically and dielectrically anisotropic

phenyl rings included in the polystyrene monomer

unit. The molecular dynamics method was used to

simulate liquid benzene as an analog of the phenyl

rings of polystyrene [6]. The configurations of the ad-

jacent rings have been determined. Calculations have

shown that at low temperatures, neighboring rings are

predominantly in a mutual T-configuration, which has

the lowest local anisotropy. Consequently, the in-

crease in local anisotropy should occur due to an in-

crease in the contributions of the sandwich configura-

tion and the planar configuration.

2 Objects and Methods of

Investigation

2.1 Objects and Simulation Details

In the present work, the molecular dynamics method

[9] was used to simulate the liquid state of three

substances: ethylbenzene (EB) as a monomer, pen-

tastyrene (polystyrene-5, PS-5, strictly speaking, it

is an oligomer), and polystyrene-25 (PS-25) as poly-

mers consisted of 5 or 25 monomer units which are

similar to EB. The purpose of the simulation was to

determine the possibility of the ordering of phenyl

rings in the substances in the isotropic phase.

Each molecule of polymer PS-5 or PS-25 was gen-

erated as an atactic isomer from monomer units which

are conventionally called left (L) and right (D) with

the use of a random number generator. In general, for

polystyrene the meaning of L and D enantiomers is

conditional, but two adjusting units form two types of

dyads which are distinguishable from each other: the

pair of the same units (either LL or DD) form meso

(m) dyad, the pair of the different units (either LD or

DL) form racemic (r) dyad. The number of L and D

units and their sequences in PS-5 and PS-25 is in con-

sistence with Bernoulli statistic [10], [11], [12] with

probabilities P(m) = P(r)=0.5(Table 1). This

is a reasonable model of industrially synthesized at-

actic polystyrene [13], [14], [15]. Examples of three

molecules of PS-25 are shown in Fig.1.

All-atoms model was chosen. Parameters of in-

termolecular atom-atom interactions for carbon and

hydrogen atoms were taken from the OPLS-AA [16]

force field. All molecules are considered as flexi-

ble. Intramolecular interactions included bond an-

gle bending, torsional (dihedral angle) rotations, and

intramolecular Lennard-Jones and Coulombic pair-

wise interaction [17]. To speed-up simulation, the

bond lengths were fixed by the SHAKE algorithm

[18]. Verlet leap-frog algorithm [19] was used with

Table 1.Chain statistics for generated atactic PS-5

and PS-25 molecules. Statistical weights of indistin-

guishable sequences of mand rdyads are summed.

Data for the sequences up to pentads are presented.

Statistical weight of dyads (%)

Dyads PS-5 PS-25 Bernoulli

sequence low

m50.3 49.5 50

r49.7 50.5 50

mm 25.6 24.2 25

mr,rm 48.8 50.4 50

rr 25.6 25.4 25

mmm 13.2 11.7 12.5

mmr,rmm 25.2 24.9 25

rrm,mrr 24.7 25.0 25

rmr 11.2 12.7 12.5

rrr 13.6 12.8 12.5

mrm 12.2 12.9 12.5

mmmm 6.6 5.4 6.25

mmmr,rmmm 13.1 12.4 12.5

mmrr,rrmm 12.7 11.8 12.5

mmrm,mrmm 11.8 13.5 12.5

mrrm 6.8 6.2 6.25

mrrr,rrrm 13.4 12.5 12.5

mrmr,rmrm 12.5 12.4 12.5

rrmr,rmrr 9.9 12.9 12.5

rrrr 6.9 6.5 6.25

rmmr 6.5 6.3 6.25

a time step of 1fs. Periodic boundary conditions

with the Ewald method [20], [21] for the accounting

of long-range electrostatic interactions were applied.

The computer program AKMD was used.

All the systems contained about an equal number

of monomer units (3550–3555 units) combined in dif-

ferent numbers of molecules, which leads to similar

cell sizes during simulation (Table 2).

Each molecule of PS-5 and PS-25 was generated

in all-trans conformation as shown in Fig.1. In the

first stage of simulation (annealing) every molecule

was equilibrated in vacuo at a temperature of 175 ◦C

during 0.5ns. During this time, the energy of the

molecules decreased, and they left all-trans conforma-

tion. After the equilibration, molecules were placed

in the simulation cell. The orientation of molecules

was random; centers of mass were placed in a reg-

ular network. The initial size of the simulation box

was chosen big enough in order the molecules do not

touch each other. The simulation box was slowly

compressed in the NVT ensemble at t= 175 ◦C until

the total energy of the system (potential plus kinetic)

reached its minimum.

After the equilibration, the systems were simu-

lated in an NpT ensemble [22] at pressure of 1atm

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2024.23.9

Andrei V. Komolkin, Sergey G. Polushin,

Vyacheslav B. Rogozhin, Alexandra A. Lezova,

Galina E. Polushina, Irina A. Silanteva

E-ISSN: 2224-2678

Volume 23, 2024

Fig.1: Examples of three molecules of atactic PS-25 generated in all-trans conformation. Meso (m) and

racemic (r) dyads are superscribed for the upper molecule

Table 2.Simulated systems and their parameters. Fluctuations of density and cell size in the NpT ensemble

during simulation are shown. Molecular mass is calculated by accounting natural abundance of isotopes as it is

recommended in the OPLS-AA force field.

Species Number of Molecular Temperature, Density, Cell size,

molecules mass, ◦C g cm−3Å

in a cell g mol−1

EB 3554 106.2 25 0.864 ±0.002 89.84 ±0.05

PS-5 711 536.8 25 1.005 ±0.001 85.75 ±0.04

175 0.926 ±0.003 88.11 ±0.08

PS-25 142 2619.8 25 1.002 ±0.002 85.06 ±0.04

175 0.969 ±0.002 86.04 ±0.06

Table 3.Calculated diffusion coefficients.

Species Temperature, diffusion coefficients,

◦C m2s−1

EB 25 1.4·10−9

PS-5 25 1.1·10−12

175 7.8·10−11

PS-25 25 9.5·10−13

175 1.6·10−12

and temperature of 175 ◦C during 25.5ns; than trajec-

tories were stored for 140.0ns; than temperature was

slowly changed to 25 ◦C and the systems were equili-

brated for 30.0ns; finally, trajectories were stored for

200.0ns.

EB was simulated at 25 ◦C. Equilibration time was

2.4ns; trajectories were stored for 4.7ns. Diffusion

coefficient of 1.4·10−9m2/s (Table 3) and mean dis-

placement of EB molecules from their initial positions

during the analyzed period (54.7 Å) proved that the

simulation time was enough to obtain averaged struc-

tural and dynamical parameters of EB molecules.

2.2 Cylindrical Distribution Function

The local structure of the substances around phenyl

rings was analyzed in terms of cylindrical distribu-

tion function (CDF). The CDF is a kind of particle-

particle distribution function gAB

2(Z, R). It shows the

distribution of particle B in the local frame of parti-

cle A. Both particles may be either a “real” particle

(carbon or hydrogen atom) or “virtual” particle as,

for example, the center of mass of carbon atoms of

phenyl rings. Widely used radial distribution func-

tion (RDF) gAB

2(R)is another kind of these func-

tions. In contrast to one-dimensional isotropic RDF,

the CDF is a two-dimensional function. This distri-

bution shows the positions of particle B in the local

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2024.23.9

Andrei V. Komolkin, Sergey G. Polushin,

Vyacheslav B. Rogozhin, Alexandra A. Lezova,

Galina E. Polushina, Irina A. Silanteva

E-ISSN: 2224-2678

Volume 23, 2024

Fig. 2: On the definition of cylindrical distribution

function (CDF) gAB

2(Z, R). Atoms of a phenyl ring

are shown in the ball-and-sticks model, part of the

main chain of the polymer is shown with sticks in the

left part of the figure. The center of mass of the car-

bon atoms of the ring is the origin O,Zaxis coincides

with the C6axis of symmetry of the ring. Explanation

of the used values of height ∆zand thickness ∆rof

cylindrical segments seen in the text

cylindrical frame of particle A. Originally, CDF was

developed for the investigation of the local structure

of nematic liquid crystals which molecules have got

rotation axis of symmetry and recently were adapted

to toluene (methylbenzene) [23]. In this work, the

function was applied to the investigation of the local

structure of phenyl rings of side-chain polymers and

ethylbenzene (as monomer unit of the polymers).

In the present work, the “virtual particles A”

are the phenyl rings of polystyrene or ethylbenzene

molecules. The local cylindrical frame (Z, R)con-

nected to the carbon atoms of the rings: six-fold sym-

metry axis (C6) of a phenyl ring is taken as Z-axis, the

origin Ois in the center of mass of carbon atoms. This

definition is shown in Fig.2. Coordinates Zand Rof

particle B are calculated in this frame. To calculate

the distribution function, the space is separated into

cylindrical segments of height ∆zand thickness ∆r.

In this work, the values of height ∆zand thickness

∆rwere chosen as 0.1Å. Both coordinates, Zand

R, are stepped with this value of 0.1Å. Number of

particles B nBcalculated inside each cylindrical seg-

ment. This corresponds to the calculation of a two-

dimensional histogram. The distribution function is a

normalized histogram according to the formula:

gAB

2(Z, R) = nB(Z, R, ∆z, ∆r)

Vsegment (Z, R, ∆z, ∆r)/NB

where nBis the average number of particles B in the

cylindrical segment Vsegment(Z, R, ∆z, ∆r). The nor-

malization coefficient is NB

V, where NBis the total

number of particles B in the simulation cell, and Vis

the volume of the cell. During a simulation, the num-

ber nBis averaged for all the phenyl rings (particles

A) of all the molecules (ensemble averaging) on each

step of the trajectory (time averaging). The value of

gAB

2(Z, R) = 1 corresponds to the mean density of the

particles B in the segment (assuming they are equally

distributed in the space). The value gAB

2(Z, R)=0

near particle A shows so-called “excluded volume”

of particle A. Excess of unity (gAB

2(Z, R)>1) shows

the position of solvation shell of the particle A by the

particles B.

To picture the CDF, a two-dimensional plot of

square pixels of size 0.1×0.1Å2is performed as

the gray-scale map. The medium gray pixel shows

the value gAB

2(Z, R)=1, the white pixel shows

gAB

2(Z, R) = 0, black pixel corresponds to the value

of 2 or more. Axis Zis plotted vertically, axis Ris

duplicated from right to left. As a result, the symmet-

ric map is been shown.

3 Results and Discussion

CDFs of centers of mass of phenyl rings (as “vir-

tual particles B”) were calculated for three systems

at 25 ◦C. They are presented in Fig.3. Excess of the

density of particles B under and below center of the

ring A is at the distance of 4.4–5.1Å from the plane of

ring A. Such distance between centers of rings allows

to ring B turn perpendicular to ring A as schemat-

ically shown in Fig.3a. This corresponds to the T-

configuration of adjacent rings [24].This result is sim-

ilar to that obtained for benzene molecules in liquid

phase [6]. Liquid benzene as well as ethylbenzene do

not show local anisotropy.

Compared to the CDF of EB, in the CDF of poly-

meric PS-5, new regions of high probability of the

neighboring rings appear diagonally left and right of

the central ring (one of them is marked with a cir-

cle in Fig.3b). Such position of the center of ring

B corresponds to the parallel-displaced configuration

(possibly, tilted) of two adjacent rings. The scheme

is included in the figure. Such a configuration is

anisotropic.

Further, in CDF of PS-25 new high-density re-

gions directly under and below the central ring at a

distance of 3.5–4.1 Å exist. These positions of the

centers of mass of ring B correspond to the “sand-

wich” configuration of two rings. Also, this config-

uration assumes π–πconjugation of aromatic rings.

The parallel orientation of phenyl rings shows their

local anisotropic ordering. The scheme is included in

the Fig.3c.

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2024.23.9

Andrei V. Komolkin, Sergey G. Polushin,

Vyacheslav B. Rogozhin, Alexandra A. Lezova,

Galina E. Polushina, Irina A. Silanteva

E-ISSN: 2224-2678

Volume 23, 2024

(a) (b) (c)

Fig.3: Cylindrical distribution functions (CDFs) of centers of mass of phenyl rings as virtual “particles” B:

(a) EB, (b) PS-5, (c) PS-25 at 25 ◦C. Inside the excluded volume of the central ring (virtual “particle” A), the

distribution of hydrogen atoms of the central ring is shown as reference points (black). Green lines schematically

show the most probable positions and orientations of two phenyl rings. The length of the green line corresponds to

the distance H–H in the para position. In the system of EB, the only possible orientations of the rings are T-shaped

(a), in the PS-5 system new parallel-displaced positions are added to the T-shaped (b), and in the PS-25 system

other new positions of sandwich-type added with respect to PS-5 (c)

At higher temperatures of 175 ◦C additional high-

density regions become less populated (Fig.4), but

still exist.

These results are following the experimental data

obtained earlier in the work [6]. The Kerr constant

K= 3.3·10−12 (cm/300 V)2of PS-5 has the lowest

average value among the measured polystyrene sam-

ples (PS-5, PS-25, PS-45, and PS-88, the two latest

are not simulated in this work) at temperature 20 ◦C,

and it is close to the constant Kof benzene (K=

2.1·10−12 (cm/300 V)2) and toluene (K= 3.9·10−12

(cm/300 V)2) [6]. So, the local orientational order of

phenyl rings in PS-5 and benzene molecules are com-

parable to each other. The fraction of phenyl rings in

the T-configuration for both substances is sufficient to

maintain the value of Kat a low level because of this

configuration is isotropic.

The formation of a coil in PS-25 leads to the re-

striction of the freedom of orientational mobility of

the phenyl rings due to steric effects. As a result, the

number of rings nBin T-configuration decreases in

favor of π–πconjugated configurations: “sandwich”

and parallel-displaced. The decreasing of nBleads

to a decrease of both the area of T-configuration and

the value of gAB

2(Z, R), which is shown in Fig.3 and

Fig.4. Therefore, the local anisotropy in polymer PS-

25 increases, and the average value Kbecomes 1.5

times greater than Kof the PS-5. Temperature de-

pendence of Kfor both PS-5 and PS-25 has the clas-

sic form K∝(1/T)[6].

4 Conclusion

In this work, the cylindrical distribution function

(CDF) gAB

2(Z, R)was adapted to analyze the mu-

tual positions and orientation of phenyl rings in the

polystyrene PS-5 and PS-25. Molecules were gener-

ated and simulated as natural atactic polymers. All-

atom model of flexible molecules was used. This

model better describes local structure of polymeric

molecules than united atoms and coarse grain mod-

els.

The modeling results and previously obtained ex-

perimental data allow us to make assumptions about

the molecular mechanisms causing changes in short-

range order (local structure) in polystyrene. Among

them:

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2024.23.9

Andrei V. Komolkin, Sergey G. Polushin,

Vyacheslav B. Rogozhin, Alexandra A. Lezova,

Galina E. Polushina, Irina A. Silanteva

E-ISSN: 2224-2678

Volume 23, 2024

(a) (b)

Fig.4: Cylindrical distribution functions (CDFs) of centers of mass of phenyl rings as virtual “particles” B:

(a) PS-5 and (b) PS-25 at 175 ◦C

1. strong influence of the length of the main chain

on the orientation of phenyl rings;

2. an increase in the freedom of rotation of

monomer units around single bonds of the main

chain and an increase in the micro-Brownian mo-

tion of the units (an increase in kinetic flexibility)

of the chain as a result of the liquid-liquid transi-

tion with increasing temperature;

3. the T-configuration is predominant in the ethyl-

benzene, PS-5, and PS-25;

4. the increase of the length of the polymeric chain

leads to the decreasing probability of isotropic

T-configuration and the increase of anisotropic

“sandwich” configuration and parallel-displaced

configuration of phenyl rings.

The future investigation will touch on unusual

temperature behavior Kerr constant of PS-45 and PS-

88 in the polystyrene melt.

Acknowledgment:

We sincerely thank A. E. Fedorov for his great

technical assistance in preparing the published

material.

References:

[1] V. Shibaev, Liquid Crystalline Polymers, in

Reference Module in Materials Science and

Materials Engineering, Elsevier, 2016, pp.

1-46.

[2] M. Eich, K. Ullrich, J. H. Wendorff, and H.

Ringsdorf, Pretransitional phenomena in the

isotropic melt of a mesogenic side chain

polymer, Polymer, Vol.25, No.9, 1984, pp.

1271-1276.

[3] Th. Fuhrmann, M. Hosse, I. Lieker, J. Rubner,

A. Stracke, and J. H. Wendorff, Frustrated

liquid crystalline side group polymers for

optical storage, Liquid Cryst., Vol.26, No.5,

1999, pp. 779-786.

[4] J. Prost and J. R. Lalanne, Laser-Induced

Optical Kerr Effect and the Dynamics of

Orientational Order in the Isotropic Phase of a

Nematogen, Phys. Rev. A, Vol.8, No.4, 1973,

pp. 2090-2093.

[5] R. Yamamoto, S. Ishihara, S. Hayakawa, and

K. Morimoto, The Kerr constants and

relaxation times in the isotropic phase of

nematic homologous series, Phys. Lett. A,

Vol.69, No.4, 1978, pp. 276-278.

[6] S. G. Polushin, V. B. Rogozhin, G. E.

Polushina, and A. V. Komolkin,

Nanostructuring Polystyrene in a Melt,

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2024.23.9

Andrei V. Komolkin, Sergey G. Polushin,

Vyacheslav B. Rogozhin, Alexandra A. Lezova,

Galina E. Polushina, Irina A. Silanteva

E-ISSN: 2224-2678

Volume 23, 2024

Nanobiotechnology Reports, Vol.17, No.1,

2022, pp. 93-97.

[7] C. A. Glandt, H. K. Toh, J. K. Gillham, and R.

F. Boyer, Effect of dispersity on the T|| (>Tg)

transition in polystyrene, J. Appl. Polym.

Science, Vol.20, No.5, 1976, pp. 1277-1288.

[8] R. F. Boyer, Pressure dependence of secondary

transitions in amorphous polymers. 1. T|| for

polystyrene, poly(vinyl acetate), and

polyisobutylene, Macromolecules, Vol.14,

No.2, 1981, pp. 376-385.

[9] M. P. Allen and D. J. Tildesley, Computer

Simulation of Liquids, Oxford University Press,

2017.

[10] G. Moad and D. Solomon, The Chemistry of

Radical Polymerization, Elsevier, 2005.

[11] M. Kobayashi, Section 2 - Structure of gels,

characterization techniques, in Gels Handbook,

Academic Press, 2001, pp. 172-412.

[12] F. Bovey, High-Resolution NMR of

Macromolecules, Elsevier, 2012.

[13] C. Ayyagari, D. Bedrov, and G. D. Smith,

Structure of Atactic Polystyrene: A Molecular

Dynamics Simulation Study, Macromolecules,

Vol.33, No.16, 2000, pp. 6194-6199.

[14] K. Matsuzaki, T. Uryu, K. Osada, and T.

Kawamura, Stereoregularity of

polystyrene-β,β-d2,Journal of Polymer

Science: Polymer Chemistry Edition, Vol.12,

No.12, 1974, pp. 2873-2879.

[15] K. Matsuzaki, T. Uryu, T. Seki, K. Osada, and

T. Kawamura, Stereoregularity of polystyrene

and mechanism of polymerization, Die

Makromolekulare Chemie, Vol.176, No.10,

1975, pp. 3051-3064.

[16] W. L. Jorgensen, D. S. Maxwell, and J.

Tirado-Rives, Development and Testing of the

OPLS All-Atom Force Field on

Conformational Energetics and Properties of

Organic Liquids, J. Am. Chem. Soc., Vol.118,

No.45, 1996, pp. 11225-11236.

[17] M. Chalaris, A. Koufou, K. Kravari, Dipole

Moment of A-agents series via Molecular

Dynamics Simulations, Molecular Sciences

and Applications, Vol.3, 2023, pp. 1-4.

[18] J.-P. Ryckaert, G. Ciccotti, and H. J. C.

Berendsen, Numerical integration of the

cartesian equations of motion of a system with

constraints: molecular dynamics of n-alkanes,

J. Comput. Phys., Vol.23, No.3, 1977, pp.

327-341.

[19] L. Verlet, Computer ”Experiments” on

Classical Fluids. I. Thermodynamical

Properties of Lennard-Jones Molecules, Phys.

Rev., Vol.159, No.1, 1967, pp. 98-103.

[20] P. P. Ewald, Die Berechnung optischer und

elektrostatischer Gitterpotentiale, Ann.

Phys.-Berlin, Vol.369, No.3, 1921, pp. 253-287.

[21] S. W. de Leeuw, J. W. Perram, and E. R. Smith,

Simulation of electrostatic systems in periodic

boundary conditions. I. Lattice sums and

dielectric constants, Proc. R. Soc. A – Math.

Phys., Vol.373, No.1752, 1980, pp. 27-56.

[22] S. Melchionna, G. Ciccotti, and B. L. Holian,

Hoover NPT dynamics for systems varying in

shape and size, Molecular Physics, Vol.78,

No.3, 1993, pp. 533-544.

[23] V. Y. Bazaikin, A. V. Komolkin, D. A.

Markelov, Journal of Molecular Liquids 383,

122188 (2023).

[24] S. E. Wheeler and J. W. G. Bloom, Toward a

More Complete Understanding of Noncovalent

Interactions Involving Aromatic Rings, J. Phys.

Chem. A, Vol.118, No.32, 2014, pp. 6133-6147.

Contribution of Individual Authors to the

Creation of a Scientific Article (Ghostwriting

Policy)

The authors equally contributed in the present re-

search, at all stages from the formulation of the prob-

lem to the final findings and solution.

Sources of Funding for Research Presented in a

Scientific Article or Scientific Article Itself

No funding was received for conducting this study.

Conflicts of Interest

The authors have no conflicts of interest to

declare that are relevant to the content of this

article.

Creative Commons Attribution License 4.0

(Attribution 4.0 International , CC BY 4.0)

This article is published under the terms of the

Creative Commons Attribution License 4.0

https://creativecommons.org/licenses/by/4.0/deed.en

_US

WSEAS TRANSACTIONS on SYSTEMS

DOI: 10.37394/23202.2024.23.9

Andrei V. Komolkin, Sergey G. Polushin,

Vyacheslav B. Rogozhin, Alexandra A. Lezova,

Galina E. Polushina, Irina A. Silanteva

E-ISSN: 2224-2678

Volume 23, 2024