# Simple Periodic Boundary Conditions for Molecular Simulation of Uniaxial Flow

We present rotating periodic boundary conditions (PBCs) for the simulation of nonequilibrium molecular dynamics (NEMD) under uniaxial stretching flow (USF) or biaxial stretching flow (BSF). Such nonequilibrium flows need specialized PBCs since the simulation box deforms with the background flow. The technique builds on previous models using one or lattice remappings, and is simpler than the PBCs developed for the general three dimensional flow. For general three dimensional flows, Dobson <cit.> and Hunt <cit.> proposed schemes which are not time-periodic since they use more than one automorphism remapping. This paper presents a single automorphism remapping PBCs for USF and BSF which is time periodic up to a rotation matrix and has a better minimum lattice spacing properties.

## Authors

• 3 publications
• 1 publication
10/15/2019

### Numerical multiscale methods and effective boundary conditions

We develop numerical multiscale methods for viscous boundary layer flow....
12/26/2021

### Stability analysis and error estimates of local discontinuous Galerkin method for convection-diffusion equations on overlapping mesh with non-periodic boundary conditions

A new local discontinuous Galerkin (LDG) method for convection-diffusion...
12/08/2017

### Periortree: An Extention of R-Tree for Periodic Boundary Conditions

Searching spatial data is an important operation for scientific simulati...
08/31/2020

### An Expedient Approach to FDTD-based Modeling of Finite Periodic Structures

This paper proposes an efficient FDTD technique for determining electrom...
11/02/2021

### FC-based shock-dynamics solver with neural-network localized artificial-viscosity assignment

This paper presents a spectral scheme for the numerical solution of nonl...
09/24/2021

### Clipping over dissipation in turbulence models

Clipping refers to adding 1 line of code A=minA,B to force the variable ...
10/18/2017

### Review of Data Structures for Computationally Efficient Nearest-Neighbour Entropy Estimators for Large Systems with Periodic Boundary Conditions

Information theoretic quantities are extremely useful in discovering rel...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

Nonequilibrium molecular dynamics (NEMD) [5, 25] techniques are one tool used to study molecular fluids under steady flow, and for instance, some of recent applications can be found in [9, 20, 19, 16, 21, 17, 22, 15, 6, 13]. However, there are special challenges in formulating the periodic boundary conditions (PBCs) in the nonequilibrium setting [3, 26, 2, 23, 24].

We consider a molecular simulation where the particles have an average flow consistent with a homogeneous background flow matrix This flow is used to simulate the micro-scale motion of a fluid with local strain rate

We denote the coordinates of the simulation box via three linearly independent vectors coming from the origin, and we write the vectors as the columns of the matrix

 Lt=[v1tv2tv3t]∈R3×3,t∈[0,∞).

To be consistent with the background flow, a particle with a coordinates , where q is the position and v is the velocity, has images with coordinates , where are triples of integers. Since the image velocity is the time derivative of its position we have

 ddt(q+Ltn)= v+ALtn,

which implies that the simulation box deforms with the flow

 ddtLt=ALt.

If the initial lattice is not chosen carefully, the resulting lattice deformation

 Lt=eAtL0

can become degenerate and lead to a particle and some of its images becoming arbitrarily close. We want to ensure that the minimum distance between a particle and its images is nonzero for all time,

 d=infn∈Z3∖0t∈R≥0||Ltn||2>0. (1)

This is necessary in order to have long-time stable periodic boundary conditions for NEMD flows.

We consider a class of PBCs based on remapping the simulation box at various times during the simulation by choosing a new set of basis vectors for the lattice that describes the simulation box. This remapping is called a lattice automorphism and can be represented as a integer matrix with determinant one. This was first used for the case of shear flow by Lees-Edwards [11] and was then extended to the case of planar elongational flow by Kraynik and Reinelt (KR) [8]. Those algorithms result in remappings which are periodic in time, though KR showed that a time-periodic remapping to the original simulation box using such matrices is impossible for general three dimensional flows. Dobson [4] and Hunt [7]

developed PBCs for general three dimensional diagonalizable flow using similar remapping technique to the KR scheme. Those schemes use more than one automorphism matrix and result in a remapping that is not time periodic. In this paper we present a rotating box algorithm applicable to uniaxial stretching flow (USF) and biaxial stretching flow (BSF) which features advantageous properties. Namely, we will show that using the class of automorphism matrices that has a pair of complex conjugate eigenvalues, we can construct a single remapping matrix algorithm which is time periodic up to a rotation matrix and whose minimum distance (

1) is larger than those of the GenKR and Hunt algorithms.

The outline of this paper is as follows. Section 2 gives the background for PBCs especially shear flow, planar elongational flow, and general three dimensional flows. Section 3 presents the rotating box algorithm, and Section 4 gives the prove that the deformed lattice obtained is not time periodic. Section 5 compares the rotating box algorithm with the existing three dimensional flow PBCs.

## 2. Background

In this section, we give a description of the existing remapping PBCs, starting with the two dimensional flows, especially, shear flow and planar elongational flow. In the case of three dimensional flows, the generalized KR (GenKR) algorithm developed by Dobson and Hunt are presented. All the algorithms follow the same procedure: given a background flow , for each time we find the appropriate integer power of the chosen automorphism matrix (or matrices) to remap the lattice .

### 2.1. Shear Flow

We first consider the shear flow case where the background matrix is given by

 A=⎡⎢⎣0ϵ0000000⎤⎥⎦.

At time , the lattice is given by

 Lt=⎡⎢⎣1tϵ0010001⎤⎥⎦L0 where L0=⎡⎢⎣100010001⎤⎥⎦.

A highly sheared box makes the computation of interparticle interactions more difficult, however this problem can be overcome by looking at the geometry of shears that are integer multiples of the box length. The Lees-Edwards (LE) boundary conditions [11] is used to prevent the simulation box from becoming too deformed. Whenever the simulation time is an integer multiple of the inverse shear rate, the simulation box is sheared by box lengths. We remap the simulation box with the matrix

 Mn=⎡⎢⎣1−10010001⎤⎥⎦n=⎡⎢⎣1−n0010001⎤⎥⎦,n∈Z

such that at a time , the simulation box lattice is

 LtMn=⎡⎢⎣1tϵ−n0010001⎤⎥⎦,n∈Z.

Since is an integer matrix with determinant equal to one, that is, the matrices and generate the same lattice. Throughout the simulation, we choose so that the stretch is at most half of the simulation box, and that this remapping process is time-periodic with period , where denote rounded to nearest integer.

### 2.2. Planar Elongational Flow

Here, the background flow matrix is

 A=⎡⎢⎣ϵ000−ϵ0000⎤⎥⎦,

meaning that the simulation box elongates in the direction and shrinks in the direction of the standard coordinate plane. To treat this case, KR proposed the use of a diagonalizable automorphism matrix that has the form

 MV=VΛ,Λ=⎡⎢⎣λ000λ−10001⎤⎥⎦,λ>0,λ≠1,

to remap the simulation box. We consider the initial lattice so that at time when we apply to the lattice basis vectors

 LtMn=etAL0Mn=etϵDΛnV−1=eAtV−1, where D=⎡⎢⎣1000−10000⎤⎥⎦,

and . Letting , the stretch of the flow remains bounded during the simulation, and in addition, it is time periodic with period For instance,

 M=⎡⎢⎣2−10−110001⎤⎥⎦

is an example of matrix which gives a good minimum distance between a particle and its images.

### 2.3. General three-dimensional (3D) flow PBCs

For a general 3D flow

 A=⎡⎢⎣ϵ1000ϵ2000−ϵ1−ϵ2⎤⎥⎦,

Dobson and Hunt proposed equivalent algorithms to control the deformation.

#### 2.3.1. Dobson’s Approach

In [4], the author develops PBCs which use two commutative automorphism matrices which have positive eigenvalues for the remapping of the simulation box. Since the matrices are commutative, they are simultaneously diagonalizable, An example of the pair of the automorphism matrices are

 M1=⎡⎢⎣111122123⎤⎥⎦and M2=⎡⎢⎣2−21−23−11−11⎤⎥⎦.

The algorithm requires that the diagonal of the logarithm of the eigenvalue matrices must be linearly independent, thus there exists solving . Now by considering the initial lattice and picking , we remark that the remapping of the simulation box with results in the remapped lattice

 ~Lt=LtMn11Mn22 =eAtL0Mn11Mn22=etAΛn11Λn22V−1=eAtV−1,

where the remaining stretch matrix

is clearly bounded for every time . Thus the minimum distance of the remapped lattice is bounded away from zero during the entire simulation.

#### 2.3.2. Hunt’s Approach

Hunt’s approach is similar to Dobson’s, using the Lenstra-Lenstra-Lovász ( [12] in place of a second automorphism matrix. As convention in this paper, we will describe Hunt’s algorithm using column vectors instead of the row vectors used in the original paper. In fact, Hunt’s PBCs consists of remapping the simulation box with the automorphism

 M=⎡⎢⎣00110−5016⎤⎥⎦, where MV=VΛ

and choosing the initial lattice basis . After applying , the remapped lattice becomes

 ~Lt=etAL0Mn1=etAΛn1V−1=eAtV−1,

where . This singe matrix is not enough to control the deformation. The reduction algorithm [12] is used to reduce the remapped lattice by finding a matrix using a high precision reduction,

 ^Lt=LLL(~Lt)=eAtV−1M2.

In comparison to the GenKR’s approach, such is automatically found on the earlier stage of the method while considering the communicative matrices. On this point, we can improve the Hunt PBCs by finding the commutative matrix manually and apply the GenKR to produce remapped lattice which minimum distance is bounded before the reduction step. The combination of this algorithm is presented in Algorithms 2 will be presented later in the paper.

## 3. Rotating Box Algorithm

In this section, we will develop PBCs for USF and BSF that are time periodic up to a rotation. We write the background flow as

 A=ϵD,where D=⎡⎢⎣10001000−2⎤⎥⎦.

Here, rather than choosing a pair of matrices with real spectrum, we will use a single matrix which has a pair of complex conjugate eigenvalues and use it to remap the simulation box.

Let us consider , a matrix that has a pair of complex conjugate eigenvalues and write its real Jordan form

 MV=VΛ where Λ=⎡⎢ ⎢⎣~η−~β0~β~η000(~η2+~β2)−1⎤⎥ ⎥⎦,

where , and in order to avoid a full rotation. Taking the logarithm of , we have

which can be decomposed as:

 log(Λ)=βB+ηD, where B=⎡⎢⎣0−10100000⎤⎥⎦.

For all time , by choosing the initial lattice , we can keep the lattice bounded by remapping the simulation box with

 ~Lt =eAtL0Mn=eϵDtΛnV−1=enβBeAtV−1,

where for . We have already seen in the planar elongational flow case (Section 2.2) that the remapped lattice is bounded and time-periodic of period . In this case, the remapped lattice is time-periodic up to the rotation matrix . For a forward simulation in time, the algorithm reads

Since the rotation matrix is bounded, we observe that the remapped lattice is also bounded during all the simulation. In the next section, we show that the rotating algorithm is not time periodic using the fact that the rotation matrix is never equal to the identity matrix for any automorphism chosen.

## 4. Non time periodicity of the lattice in the rotating box algorithm

As mentioned above, for the class of automorphism matrices with real eigenvalues, it has been shown in [8] that it is impossible to construct KR PBCs with a time periodic lattice for USF or BSF. In this section, we will extend this demonstration to the class of automorphism matrices which have complex conjugate eigenvalues. Namely, we show in the following corollary that although the rotation algorithm applied to USF or BSF is time-periodic up to a rotation matrix, there is no choice of where the period of the remapping aligns with that of the rotation. In other word, we show that rotation matrix is not equal to the identity matrix for , or is not a equal of times a rational number

 β=tan−1(~β~η)≠2πmn,n,m∈Z,n≠0,

for any consider in Section 3, i.e with complex eigenvalues one of the eigenvalue of is not equal to .

We start by reminding that are the roots of the characteristic polynomial of , and write these roots in the polar coordinate as , . Let us first show the following lemma:

###### Lemma 1.

A matrix with complex eigenvalues as define above has if it has at least an eigenvalue equal to .

The proof of Lemma 1 requires the use of the following results. Let us consider , the Euler totient function where is the number of positive integers that are relatively prime to . A scalar is said to be algebraic over a field if there exists elements of , not equal to zero, such that

 α0+αa1+⋯+αiai=0,

and is the degree of the irreducible characteristic polynomial. We refer the reader to [10, Chapter 4] or any introduction to Algebra book for the background about the definitions used in this section. Then we have:

###### Theorem 1.

[18, Theorem 3.11] For and ,

 deg{tan2mπn}=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩φ(n) for gcd(n,8)<4,φ(n)2 for gcd(n,8)=4,φ(n)4 for gcd(n,8)>4.
###### Theorem 2.

[1, Theorem 16.8.5] For the splitting field of an irreducible cubic polynomial over a field and the discriminant of ,

• If is a square in , the degree of the extension field over is three

• If is not a square in , the degree of the extension field over is six.

We determine the degree of the algebraic integer in the following lemma:

###### Lemma 2.

is an algebraic integer of degree at most six.

###### Proof.

Since and are elements of the splitting field of the irreducible polynomial , we have that is also an element of . By Theorem 2, has a degree as most six in and so does . ∎

Let us prove Lemma 1 by finding all coefficients of the characteristic polynomial of for which

###### Proof.

Let us assume that and find the possible by using Theorem 1 and the Theorem 2 which guarantee that is an algebraic integer of degree at most six. Thus using [14], we find all that satisfy the following

 φ(n)≤6 for gcd(n,8)<4, φ(n)2≤6 for gcd(n,8)=4, φ(n)4≤6 for gcd(n,8)>4,

and report all and in Table 1.

Then after few computing we find that

 h =1+2r3cos2πmnr2,k=r3+2cos2πmnr,

and plugging in from the Table 1 and such that , we remark that if . In result, has at least one eigenvalue equal to , since , or for the latter values of . ∎

In sum, we derive the main corollary of this section:

###### Corollary 1.

The rotating box algorithm cannot give a time-periodic simulation box for any choice of integer commutative complex conjugate matrix.

###### Proof.

Using Lemma 1, we know that only matrices with an eigenvalue equal to one 1 have a rotational part that is a root of unity. However, those matrices are themselves pure rotations and have no use for the PBCs since they cannot control the stretching caused by the underlying background flow. ∎

## 5. Comparison of the three dimensional algorithms

In this section, we compute the minimum distance of the particle images for our algorithm and compare it with the GenKR using, Hunt’s and Dobson’s automorphism matrices.

To compute the minimum distance between a particle and its images when the rotating box PBCs is applied, we propose the matrix

 M =⎡⎢⎣0−21110010⎤⎥⎦,

which has a pair of complex conjugate eigenvalues with positive real part. Then, the initial lattice is given by

 L0=⎡⎢⎣0.77260−0.2083−0.260860.434420.48424−0.35555−0.141060.84978⎤⎥⎦,

which implies that, given the standard lattice with the coordinate , the plane is rotated counterclockwise by approximately degrees, and by degrees.

Moreover for the GenKR algorithm, we keep the automorphism matrices and the initial lattice given in the original paper. The commutative matrices and associated with the initial orthonormal lattice basis which determinant is equal to one, are given by

 M1=⎡⎢⎣111122123⎤⎥⎦,M2=⎡⎢⎣2−21−23−11−11⎤⎥⎦and, L0=⎡⎢⎣0.59101−0.736980.327990.736980.32799−0.591010.327990.591010.73698⎤⎥⎦.

For Hunt’s formulation, we find a second automorphism matrix which has positive eigenvalues and is commutative with the matrix given in the original paper. The commutative matrices and the normalized initial lattice respectively read

 M1=⎡⎢⎣00110−5016⎤⎥⎦,M2=⎡⎢⎣311−5−2−4114⎤⎥⎦and, L0=⎡⎢⎣0.522762.639413.32590.522760.336190.21620.522760.1610.049584⎤⎥⎦.

Then we graph the minimum distance for the three dimensional algorithms in figure 1, when the stretch is . We can observe in the graph that the minimum distance curve in rotating box PBCs case presents a clear pattern of periodicity. In addition, the minimum distance for all the simulation for the rotating box algorithm is approximately 1.0271 compare to 0.9054 in GenKR’s case.

## 6. Conclusion

Kraynik-Reinelt proved that it is impossible to find time periodic PBCs for general three dimensional flow, using matrices with real eigenvalues. In this paper, we show that by using an matrix with complex eigenvalues, we can create an algorithm that is time-periodic up to a rotation matrix for USF and BSF. Although we show that the rotations never align, the regularity of the remapping make this algorithm more straightforward than the existing ones. These PBCs also offer a better minimum distance between a particle and its images.

## References

• [1] M. Artin. Algebra. Pearson Prentice Hall, 2011.
• [2] András Baranyai and Peter T. Cummings. Steady state simulation of planar elongation flow by nonequilibrium molecular dynamics. The Journal of Chemical Physics, 110(1):42–45, 1999.
• [3] P.J. Daivis, M.L. Matin, and B.D. Todd. Nonlinear shear and elongational rheology of model polymer melts by non-equilibrium molecular dynamics. Journal of Non-Newtonian Fluid Mechanics, 111(1):1 – 18, 2003.
• [4] Matthew Dobson. Periodic boundary conditions for long-time nonequilibrium molecular dynamics simulations of incompressible flows. The Journal of Chemical Physics, 141(18):184103, 2014.
• [5] Denis J. Evans and Gary P. Morriss. Statistical mechanics of nonequilibrium liquids. ANU E Press, Canberra, 2007.
• [6] James Ewen, D. Heyes, and Daniele Dini. Advances in nonequilibrium molecular dynamics simulations of lubricants and additives. Friction, 6, 02 2018.
• [7] Thomas A. Hunt. Periodic boundary conditions for the simulation of uniaxial extensional flow of arbitrary duration. Molecular Simulation, 42(5):347–352, 2016.
• [8] A.M. Kraynik and D.A. Reinelt. Extensional motions of spatially periodic lattices. Int. J. Multiphase Flow, 18(6):1045 – 1059, 1992.
• [9] Philipp S. Lang, Benedikt Obermayer, and Erwin Frey. Dynamics of a semiflexible polymer or polymer ring in shear flow. Phys. Rev. E, 89:022606, Feb 2014.
• [10] Serge Lang. Algebra. Springer, New York, NY, 2002.
• [11] A W Lees and S F Edwards. The computer study of transport processes under extreme conditions. J. Phys. C Solid State, 5(15):1921, 1972.
• [12] Arjen Lenstra, H. Lenstra, and László Lovász. Factoring polynomials with rational coefficients. Mathematische Annalen, 261, 12 1982.
• [13] Zhen Li, Shiyun Xiong, Charles Sievers, Yue Hu, Zheyong Fan, Ning Wei, Hua Bao, Shunda Chen, Davide Donadio, and Tapio Ala-Nissila. Influence of thermostatting on nonequilibrium molecular dynamics simulations of heat conduction in solids. The Journal of Chemical Physics, 151(23):234105, 2019.
• [14] N. S. Mendelsohn. The equation . Mathematics Magazine, 49(1):37–39, 1976.
• [15] A. G. Menzel, P. J. Daivis, and B. D. Todd. Equilibrium and nonequilibrium molecular dynamics methods to compute the first normal stress coefficient of a model polymer solution. Phys. Rev. Fluids, 5:084201, Aug 2020.
• [16] David A. Nicholson and Gregory C. Rutledge. Molecular simulation of flow-enhanced nucleation in n-eicosane melts under steady shear and uniaxial extension. The Journal of Chemical Physics, 145(24):244903, 2016.
• [17] Akihiro Nishioka, Tatsuhiro Takahashi, Yuichi Masubuchi, Jun-ichi Takimoto, and Kiyohito Koyama. Description of uniaxial, biaxial, and planar elongational viscosities of polystyrene melt by the k-bkz model. Journal of Non-Newtonian Fluid Mechanics, 89:287–301, 03 2000.
• [18] Ivan Niven. Irrational Numbers, volume 11. Mathematical Association of America, 1 edition, 1985.
• [19] Thomas C. O’Connor, Nicolas J. Alvarez, and Mark O. Robbins. Relating chain conformations to extensional stress in entangled polymer melts. Phys. Rev. Lett., 121:047801, Jul 2018.
• [20] Thomas C. O’Connor, Ting Ge, Michael Rubinstein, and Gary S. Grest. Topological linking drives anomalous thickening of ring polymers in weak extensional flows. Phys. Rev. Lett., 124:027801, Jan 2020.
• [21] AnaSofia Oliveira, Giovanni Ciccotti, Shozeb Haider, and Adrian Mulholland. Dynamical nonequilibrium molecular dynamics reveals the structural basis for allostery and signal propagation in biomolecular systems. The European Physical Journal B, 94:144, 07 2021.
• [22] Clark Templeton, R. Elber, Mauro Ferrario, and Giovanni Ciccotti. A new boundary driven nemd scheme for heat and particle diffusion in binary mixtures. Molecular Physics, page e1892849, 03 2021.
• [23] B. D. Todd and Peter J. Daivis. Nonequilibrium molecular dynamics simulations of planar elongational flow with spatially and temporally periodic boundary conditions. Phys. Rev. Lett., 81:1118–1121, Aug 1998.
• [24] B. D. Todd and Peter J. Daivis. The stability of nonequilibrium molecular dynamics simulations of elongational flows. The Journal of Chemical Physics, 112(1):40–46, 2000.
• [25] Billy D. Todd and Peter J. Daivis. Nonequilibrium Molecular Dynamics: Theory, Algorithms and Applications. Cambridge University Press, 2017.
• [26] Wen-Sheng Xu, Jan-Michael Y. Carrillo, Christopher N. Lam, Bobby G. Sumpter, and Yangyang Wang. Molecular dynamics investigation of the relaxation mechanism of entangled polymers after a large step deformation. ACS Macro Letters, 7(2):190–195, 2018.