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 microscale 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
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
which implies that the simulation box deforms with the flow
If the initial lattice is not chosen carefully, the resulting lattice deformation
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,
(1) 
This is necessary in order to have longtime 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 LeesEdwards [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 timeperiodic 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
At time , the lattice is given by
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 LeesEdwards (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
such that at a time , the simulation box lattice is
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 timeperiodic with period , where denote rounded to nearest integer.
2.2. Planar Elongational Flow
Here, the background flow matrix is
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
to remap the simulation box. We consider the initial lattice so that at time when we apply to the lattice basis vectors
and . Letting , the stretch of the flow remains bounded during the simulation, and in addition, it is time periodic with period For instance,
is an example of matrix which gives a good minimum distance between a particle and its images.
2.3. General threedimensional (3D) flow PBCs
For a general 3D flow
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
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
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 LenstraLenstraLová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
and choosing the initial lattice basis . After applying , the remapped lattice becomes
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,
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
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
where , and in order to avoid a full rotation. Taking the logarithm of , we have
which can be decomposed as:
For all time , by choosing the initial lattice , we can keep the lattice bounded by remapping the simulation box with
where for . We have already seen in the planar elongational flow case (Section 2.2) that the remapped lattice is bounded and timeperiodic of period . In this case, the remapped lattice is timeperiodic 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 timeperiodic 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
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
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 ,
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
and report all and in Table 1.
1  1, 2 

2  3, 6, 12, 16, 24 
4  5, 10, 20, 32, 40, 48 
6  7, 9, 14, 18, 28, 36, 56, 72 
Then after few computing we find that
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 timeperiodic 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
which has a pair of complex conjugate eigenvalues with positive real part. Then, the initial lattice is given by
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
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
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
KraynikReinelt 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 timeperiodic 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 nonequilibrium molecular dynamics. Journal of NonNewtonian Fluid Mechanics, 111(1):1 – 18, 2003.
 [4] Matthew Dobson. Periodic boundary conditions for longtime 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 AlaNissila. 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 flowenhanced nucleation in neicosane melts under steady shear and uniaxial extension. The Journal of Chemical Physics, 145(24):244903, 2016.
 [17] Akihiro Nishioka, Tatsuhiro Takahashi, Yuichi Masubuchi, Junichi Takimoto, and Kiyohito Koyama. Description of uniaxial, biaxial, and planar elongational viscosities of polystyrene melt by the kbkz model. Journal of NonNewtonian 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] WenSheng Xu, JanMichael 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.
Comments
There are no comments yet.