PULSEDYN—A dynamical simulation tool for studying strongly nonlinear chains

Computer Physics Communications 239 (2019) 134–149

Contents lists available at ScienceDirect

Computer Physics Communications journal homepage: www.elsevier.com/locate/cpc

PULSEDYN—A dynamical simulation tool for studying strongly nonlinear chains✩ ∗

Rahul Kashyap a , , Surajit Sen a,b a b

Department of Physics, State University of New York at Buffalo, Amherst, NY 14260-1500, USA Department of Physics, Brock University, St. Catharines, Ontario, Canada L2S 3A1

article

info

Article history: Received 26 July 2018 Received in revised form 24 January 2019 Accepted 26 January 2019 Available online 15 February 2019 Keywords: Fermi–Pasta–Ulam system Nonlinear dynamics Relaxation Toda lattice

a b s t r a c t We introduce PULSEDYN, a particle dynamics program written in C++, to solve the equations of motion of many-body nonlinear systems in one dimension. PULSEDYN contains a suite of commonly used potentials (Fermi–Pasta–Ulam–Tsingou, Toda, Morse and Lennard-Jones) and solvers (velocity-Verlet and Gear 5th order predictor–corrector) for particle dynamics. PULSEDYN is designed to perform scientifically accurate simulations using these calculations. For accessing the built-in features of PULSEDYN, no knowledge of programming is expected, apart from the creation of a parameter file using predefined commands for running the executable. Therefore, simulations for research projects can be performed with minimal code writing using PULSEDYN. PULSEDYN is distributed under the GNU GPL v3 license and the source code and the executable are available on Github. The writing style emphasizes organization and legibility. Therefore, we anticipate that it would serve as a good template for users who may wish to adapt the code to their specific needs. In this manuscript, we first discuss PULSEDYN and its features in detail. Then, we show results reproduced from literature using PULSEDYN for the following cases (i) Soliton propagation and collisions in the (integrable) Toda lattice (ii) Recurrence phenomena, decay of localized excitations and solitary wave collision in the Fermi–Pasta–Ulam–Tsingou lattice and (iii) Solitary wave propagation in the Morse and Lennard-Jones lattices. Finally, we present a new result on a problem of fundamental historical importance using PULSEDYN. We show by means of explicit specific heat calculations, that in the limit of strong nonlinearity, the quartic Fermi–Pasta–Ulam–Tsingou system approaches equipartition at late times. Program summary Program Title: PULSEDYN Program Files doi: http://dx.doi.org/10.17632/fc8483t7jz.1 Licensing provisions: GNU General Public License 3 (GPL) Programming Language: C++ Supplementary Material: user manual, documentation Nature of problem: We solve a Newtonian particle dynamics problem in one dimension for a variety of commonly used nonlinear manybody systems. The goal is that the user should be able to run scientifically accurate simulations with minimal effort. Benchmarking should also be easy to perform by different users running simulations across different platforms. Solution method: The code uses a parameter file interface with a set of commands for a range of prebuilt functionalities. By providing the code and the needed parameter file to the user we hope to make scientifically accurate simulations easier to perform. We provide the Gear 5th order predictor– corrector algorithm as well as the velocity-Verlet algorithm and the following potentials — the Toda potential, Morse Potential, Lennard-Jones Potential and a polynomial potential up to 4th order. To make benchmarking consistent we release the software as open source and free to distribute. © 2019 Elsevier B.V. All rights reserved.

1. Introduction ✩ This paper and its associated computer program are available via the Computer Physics Communication homepage on ScienceDirect (http://www.sciencedirect. com/science/journal/00104655). ∗ Corresponding author. E-mail address: [email protected] (R. Kashyap). https://doi.org/10.1016/j.cpc.2019.01.025 0010-4655/© 2019 Elsevier B.V. All rights reserved.

The importance of computational power and techniques to analyze many-body nonlinear systems (MBNS) cannot be overstated. Major discoveries in many-body physics such as the famous Fermi– Pasta–Ulam–Tsingou (FPUT) paradox [1–3] have been explored via numerical means. Sometimes, numerical techniques have been

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

used in conjunction with analytical descriptions of nonlinear phenomena [4–11]. A great number of times, numerical techniques are our only tool to explore the physics of a system when analytical techniques become unwieldy [12–22]. Some books on nonlinear systems such as Ref. [23,24] give a more detailed discussion. Given the importance of computational methods to investigate MBNS, it is important to develop a particle dynamics (PD) code to accurately solve the coupled dynamical equations in MBNS over long time scales while retaining strict control over error propagation. PULSEDYN can be used to carry out accurate dynamical simulations of MBNSs (see for example Ref. [25,26]). The differential equation solvers currently available in the code are the Gear 5th order predictor–corrector [27] and the velocity-Verlet algorithms [28]. The latter is a commonly used modification of the Störmer–Verlet algorithm [29]. Interested readers can find comparisons of these algorithms in standard texts such as [30]. PULSEDYN is currently available for Windows (testing performed on Windows 10) and Linux (testing performed on Ubuntu 14.04 LTS). In addition to containing the potentials and solvers listed above, it also contains features such as external forcing and dissipation [17]. It would allow the user to vary a wide variety of parameters and set up simulations with publication level precision. The code is open source under the GNU GPL v3 license. The executable and code details can be found in Tables 3 and 4 in the Appendix. Hence, users with experience in C++ will be able to modify PULSEDYN for specific purposes across different platforms. Further, we believe it could serve as a tool to benchmark other codes designed for the same purpose. In Section 2 we discuss the class of problems that the program is designed to solve. We then describe the programming details and details of the parameter file interface used to set up simulations in Sections 3 and 4. We show results obtained from the code in Section 6 for test cases of the potentials built into the system. We test the four potentials provided in the program and show that the numerical results match the results in the literature. First, we show results for the Toda chain system, which is integrable, and compare the well known Toda soliton solution from simulations to theoretical predictions [3,4,31]. Here, by soliton we mean a non-dispersive energy bundle that does not interact with another identical energy bundle except for a trivial phase lag. The lag is induced by a slowdown during the collision, which is a property of integrable systems. Next, we consider non-integrable systems starting with the FPUT system and show the recurrence phenomenon, localized excitations and solitary wave (SW) collisions in the lattice [1,2,32,33]. We then show SW propagation through the Morse and Lennard-Jones lattices and demonstrate that the simulations done using PULSEDYN are in agreement with established results from the literature [25,26]. We distinguish here between solitons and SWs. Typically, SWs show significant scattering when they collide [33–35]. Solitons are also typically obtained as exact solutions of integrable systems such as the Toda chain [3,4,31]. In Section 7, we turn our attention to a problem of historical importance — namely the road to equipartitioning of energy in the FPUT system [1,2]. We show that the results we have obtained using PULSEDYN provide strong evidence that the purely nonlinear β -FPUT system goes to equilibrium at late times and energy is equipartitioned in the system [36–40] 2. Physical systems considered We describe MBNSs by defining a Hamiltonian and then writing down the equations of motion. This gives rise to a set of coupled differential equations which are solved numerically using a well chosen integration algorithm depending on the problem at hand [30,41]. We have included the Toda, FPUT, Morse and LJ potentials in PULSEDYN. The systems are described in the code

135

by an inter-particle potential of the form Vi,i+1 (k1 , k2 , k3 ) where i and i + 1 are adjacent particles in the chain. The form of Vi,i+1 is described below for each of the potentials with nomenclature that is different than that in the original works for reasons pertaining to the code architecture, which is addressed in Section 3. (a) Toda Potential [3,4], Vi,i+1 (k1 , k2 ) =

k1 k2

e−k2 (xi+1 −xi ) + k1 (xi+1 − xi ) −

k1 k2

, k3 = 0.

(1)

(b) FPUT Potential (α + β model) [1,2], where α , β were used in the original and many subsequent works to denote the prefactors of the cubic and quartic terms in the potential. We use the constants k1 , k2 and k3 for the quadratic, cubic and quartic terms, respectively, as shown below, Vi,i+1 (k1 , k2 , k3 ) = k1 (xi+1 − xi )2 + k2 (xi+1 − xi )3 + k3 (xi+1 − xi )4 . (2) (c) Morse Potential [26,42], Vi,i+1 (k1 , k2 , k3 ) = k1 (e−k2 (xi+1 −xi ) − 1)2 , k3 = 0,

(3)

(d) and the LJ Potential [26,43],

[( Vi,i+1 (k1 , k2 , k3 ) = k2

k1

)12

k1 + xi+1 − xi

( − 2

k1

)6

k1 + xi+1 − xi

] + 1 , k3 = 0.

(4)

In the potentials listed above, k1 , k2 and k3 set the parameters. Observe that k3 ̸ = 0 only in Eq. (2). For the LJ potential we have used the form in [26]. The code solves the following force equation mi d2 xi dt 2

= −V ′ (xi − xi−1 ) + V ′ (xi+1 − xi ),

(5)

where, mi and xi are the mass and displacement of the ith particle. V ′ (xi+1 − xi ) is the derivative dV (r)/dri where ri ≡ xi+1 − xi . PULSEDYN is capable of handling driven dissipative systems described by mi d2 xi dt 2

= −V ′ (xi − xi−1 ) + V ′ (xi+1 − xi ) + Fdiss + Fext .

(6)

In Eq. (6), Fdiss and Fext are the velocity dependent dissipative and sinusoidal driving forces, respectively. These additional terms may be added to individual particles as well as to the entire system. We briefly discuss the issue of numerical instability below. If two particles come sufficiently close to each other, the magnitude of the restoring force can become very large in a manner that is dictated by the equations of motion. This rapid increase can potentially cause numerical instability due to floating point calculational errors. For instance, if the energies are set too high, the time scales of interest in the dynamical studies may become too small to be resolved by the time step set being used in the parameter file, which can cause the simulation to output unreliable results. Care must be taken while setting the values of the parameters since they can easily lead to instabilities. For instance, in the case of the LJ potential, k1 sets the bond length in the chain explicitly. In the case of the FPUT potential, the cubic term can potentially cause instabilities. In our studies presented here, we have set only positive k1 , k2 and k3 . However, negative values of the parameters are not necessarily unphysical and the code does not prohibit the user from setting negative parameter values. Hence, the responsibility of setting physically meaningful parameters rests on the user. The user must also distinguish between Hamiltonian and nonHamiltonian systems when running simulations. The velocityVerlet algorithm is unstable for non-conservative systems [41].

136

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

Table 1 Parameter file commands used to run the code.

model: method: systemsize: timestep: recsteps: printint: init: boundary: mass: force: dissipation:

Specify potential type and parameters Choose integration algorithm Set number of particles in the chain Set the value of dt used by the integration algorithm Set total number of recorded steps Set the time iterations after which data is recorded Set initial conditions Set boundary conditions Set masses of individual particles Set external force on the system Set a velocity dependent dissipation

Therefore, adding external forces to the system would cause instability in the simulations if the velocity-Verlet algorithm is implemented naïvely. In PULSEDYN, if the user wishes to add external driving or dissipation to the system, the Gear 5th order corrector– predictor method must be specified to integrate the equations of motion correctly. If the velocity-Verlet algorithm is called in the parameter file, the code ignores Fdiss and Fext and solves Eq. (5) using the initial conditions and parameters set in the parameter file.

where ki , i = 1, 2, 3 are the parameters of the potentials listed in Section 2. The options for modelType are (a) (b) (c) (d)

toda for the toda potential, fput for the FPUT potential, morse for the Morse potential, And lennardjones for the LJ potential.

Only the FPUT potential uses k3 during calculations. However, for the rest of the potentials, the value of k3 must be set to zero explicitly to avoid any unexpected behavior. The total number of particles in the system may be set using the systemsize command using the following format.

systemsize: N As an example, to set up a 100 particle system enter the command as shown below.

systemsize: 100 The integration algorithm can be set by the command method using the format below.

method: methodName 3. Software functionalities PULSEDYN is designed to set up simulations with minimal effort. A simulation can be started by setting up the required simulation details using a few pre-defined commands in the parameter file, which must be named parameters.txt, and double clicking on the executable or by running the executable in a terminal window. The options which can be set in the parameter file are (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k)

Potential type, Integration algorithm, Number of particles N, Integration algorithm, Integration time step dt, Data sampling interval Tsamp , Number of recorded data snapshots in time Trec , Initial conditions, Boundary types, Particle masses, Driving forces and dissipation.

A detailed explanation of how to specify these commands is given in the user manual provided with the code on the Github repository. The list of commands is shown in Table 1. The standard format of the parameter file commands is as follows.

command: A B C D . . . In the above template, the command name is always followed by colon. A, B, C, D etc. are values to be passed to PULSEDYN, are specific to the commands and they are detailed below. Spaces must be given between the semicolon and value of A and between the different values themselves. Care must be taken to enter all the values required by each command. If some values need not be set, placeholder values must be given. Specifying less or more number of values than the default number of values following a command can cause unexpected behavior and crashes. It is also recommended that the particular order of commands in Table 1 be followed when setting up the parameters file as unexpected behavior may arise if commands are issued in a conflicting order. Any exceptions to this recommendation will be pointed out as they arise. The command model is used to set up the potential type along with the parameters in the following format.

model: modelNamek1 k2 k3

For the velocity-Verlet algorithm set methodName to velocityverlet and for the Gear 5th order predictor–corrector method, set methodName to gear5. As an example, we may set the integration algorithm as the velocity-Verlet algorithm in the following manner.

method: velocityverlet The integration time step dt can be set using the timestep command. Setting an integration time-step for a simulation can be a tricky matter for nonlinear systems due to the multiple possible timescales depending on the system and the phenomena under investigation. To maintain strict control over error, the integration time-step may have to be reduced quite drastically. This causes a corresponding increase in computational cost. Further, there also arises the issue of large data accumulation due to the high time resolution. While smaller time-steps may be needed for error controlling purposes, data at every single iteration may not be needed to perform further analyses. Therefore, it is necessary to have the ability to sample data at well-chosen time intervals. For this purpose, the commands printint and recsteps are provided. Together, they can be used to set the time sampling interval in time, Tsamp and total number of recorded data snapshots in time Trec so that the size of the data files remains manageable. The commands can be set using the format given below.

timestep:dt printint:Tsamp recsteps:Trec These three quantities set the total time in system time units Tsys , up to which the simulation is to be performed. With an integration time-step of dt and a data sampling interval of Tsamp , data is recorded at intervals of Tsamp dt in system time. If a total of Trec such data points are recorded, then Tsys = Trec Tsamp dt amount of time has passed in system units. Consider the example below.

timestep: 0.001 printint: 10 recsteps: 10000 In this example, dt = 0.001, Tsamp = 10 and Trec = 10000. Each iteration advances system time by amount dt. Data is recorded every Tsamp number of iterations. Therefore, each recorded time snapshot in the data files correspond to an advancement of system

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

time by Tsamp dt = 0.01. Since a total of Trec = 10000 such snap shots are recorded, the final simulation time in system time units is then Tsys = 100. By default, Tsamp is set to be 1/dt. This records data every single unit of time in system units. Initial conditions can be specified with the command init using one of three formats. The first format is for setting initial conditions at each of the particles. This instance uses the format

Table 2 Format for initial conditions to be read from file. x1 x2 x3 . . . . . xN

init: particleNumber perturbationType strength To give particles initial velocities perturbationType is set to vel. For initial displacements, perturbationType is set to pos. The following two examples show examples of an initial displacement of particle 40 by 0.1 and an initial velocity of −0.1 assigned to particle 60.

init: 40 pos 0.1 init: 60 vel -0.1 The second format also assigns initial conditions to individual particles with values drawn from a uniform distribution of numbers over an interval. The format of the command in this case is given below.

init: particleNumber type random lowVal highVal In the above format, particleNumber is the particle for which initial conditions need to be specified. The keyword type refers to either initial displacement or velocity assigned by pos and vel respectively. The keyword random specifies that a random number must be assigned. The random number to be used as the initial condition value is chosen between the two numbers lowVal and highVal. An example is shown below for assigning a random velocity drawn from a uniform distribution of numbers between 0.1 and 0.5 to particle 20.

init: 20 vel random 0.1 0.5 Finally, initialization can also be done by reading in values from a file containing initial conditions in the format shown in Table 2. In this case, the data for the initial conditions must be present as three columns. The first column corresponds to position, the second to velocity and the third to acceleration. The number of rows corresponds to the number of particles in the chain. The format for inputting initial conditions from file is shown below.

137

v1 v2 v3 . . . . .

vN

a1 a2 a3 . . . . . aN

Boundary conditions are set using the boundary command in the following format.

boundary: side type The keyword side refers to either the left end of the chain i.e. boundary condition for particle 1, or the right end of the chain i.e. boundary condition for particle N in a chain of size N. For the right boundary, side is set to right and for the left boundary, side is set to left. Three types of boundaries are available — fixed, open and periodic. They can set by assigning type to be fixed, open and periodic respectively. The left and the right boundaries can be set independently to be either fixed or open. However, choosing periodic boundaries for either one of the left or right boundaries will force the second boundary to be periodic as well. An example is shown below with the left and right boundaries set to be open and fixed respectively.

boundary: left open boundary: right fixed The masses of individual particles are set to 1.0 by default. However, mass impurities can be introduced by changing masses of individual particles using the mass command. The command format is shown below.

mass: particleNumber massValue In the above format, massValue assigns a new mass to the particle specified by particleNumber. An example is shown below for assigning the mass of 2.0 to particle 10.

init: file fileName.extension

mass: 10 2.0

Notice that the file name must be provided with its full extension. As an example, if the initial conditions are stored in the file init.txt, the command would be given in the format shown below.

Finally, external driving (Fext ) and dissipation (Fdiss ) can be added to the system. The format for specifying external forces on the system is given below.

init: file init.txt

force: particleNum forceType A t1 t2 t3 ramp

In this case, it is recommended to not specify the number of particles explicitly using the systemsize command. If the number of particles are set correctly, then no unexpected behavior would be seen. If the number of particles are not equal to the number of rows in the initial conditions file, the systemsize command is given precedence. If fewer particles are specified than the number of rows in the initial conditions file, the values from the file are read in till the number of specified particles are reached and the rest of the values are ignored. If the number of particles are specified to be more than the number of rows in the initial conditions file, the initial conditions file is read to the last row and the remaining particles are assigned values of 0 for initial displacement and velocity. It must be mentioned that no provisions for including comments in the initial conditions file is provided. Blank lines will be skipped, but the usage is strongly discouraged. The format for the initial conditions file is shown in Table 2.

In the above format, particleNum is the index of the particle to which the force has to be applied. One can specify an individual number for a particle, or set particleNum to all to apply the driving force on all particles. The type of driving force function can be set by specifying forceType. Two driving functions are provided - sine and cosine. The amplitude of the driving force can be set by specifying a value for A. The quantities t1 , t2 and t3 set the start time, time period and the stop time of the driving force. These times must be specified in system time units. The value of t2 is used to calculate f0 = 1/t3 where f is the frequency of the driving force. The quantity ramp is used to specify a time changing frequency of the form f (t) = f0 + ramp × t, where t is system time. When ramp is set to 0, we recover the usual sinusoidal driving force. The case ramp ̸ = 0 is also referred to as ‘chirping’ and is a more general form of the sinusoidal driving force. It has been used to precipitate intrinsic localized modes in nonlinear chains [16,17].

138

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

Fig. 1. An example of the parameter file used for the seeded LNE in the FPUT chain is shown here. The results are discussed in Section 6.2.

An example is shown below for the driving force of the form, Fext = A sin(2π ft) on particle 10.

force: 10 sine 0.2 0.0 0.5 100.0 0.10 In the above example, A = 0.2, t1 = 0.0, t2 = 0.5, t3 = 100.0 and ramp = 0.10. Therefore, f0 = 1/t2 = 2.0 and ramp = 0.10 and f = f0 + ramp × t = 2.0 + 0.10t. The force then becomes Fext = 0.2 sin(2π ft) and is applied from t = 0 to t = 100.0 in system time units. Dissipation can be set as a velocity dependent damping force given by Fdiss = −γ v where γ is the damping coefficient and v is the instantaneous velocity. Dissipation can be set up in the following manner.

dissipation: particleNum γ Dissipation can be set to individual particles by specifying the particle index as the particleNum value and the corresponding γ value. The keyword all can be used to set dissipation for every particle in the chain. An example is shown for dissipation applied to the entire chain with γ = 0.1.

dissipation: all 0.1 Data is written to file every printint number of iteration for recsteps number of time snap shots. The following data files are written (a) (b) (c) (d) (e) (f)

position.dat for displacement data, velocity.dat for velocity data, acceleration.dat for acceleration data, ke.dat for kinetic energy data, pe.dat for potential energy data, restart.dat records displacement, velocity and accelera-

tion at the last time iteration, (g) mass.dat for mass values of each particle and (h) totalEnergy.dat for total system energy at each recorded iteration. A sample parameter file is also shown in Fig. 1 for setting the parameters for a seeded localized nonlinear excitation (LNE) in the FPUT chain. The system is discussed in Section 6.2. The files are written for position, velocity, acceleration and kinetic energy as arrays of size Trec × N, where N is the size of the chain. The masses of the particles are recorded as an (N × 1) array of numbers. The potential energies of the springs are recorded as a Trec × (N + 1) array. As detailed previously, in the case of fixed boundaries at both ends, there are (N + 1) springs. If one of the boundaries is open or if both are open, then there are N and N − 1 springs respectively. The code simply prints zeros on the side of the open boundary in the file. In the case of the periodic boundaries, there are N springs. The spring connecting particle 1 and particle N is considered the leftmost spring and the (N + 1)th column in the file is given a value of 0 to avoid over counting.

It must be noted that imposing open boundaries or periodic boundaries can make the chain susceptible to drifting. For instance, a 2 particle simple harmonic oscillator with open ends would show a center-of-mass drifting due to a drift velocity as well as periodic oscillations if one of the particles is given an initial velocity perturbation in one direction. The code also prints out total energy of the system at each recorded step as an array of size Trec × 1. The last file recorded is a restart file. It records the position, velocity and the acceleration as an (N × 3) array. The first column corresponds to position, the second to velocity and the third to acceleration. This is identical to the format required by the code to read initial conditions from file which was detailed earlier. Since the restart file follows the same format as that for the initial conditions file, the feature to import data from file may be also leveraged to continue runs from the last known point in time of a previous simulations. This can be achieved by copying the restart file as the initial conditions file for the continuation of a simulation. 4. Code architecture 4.1. Organization This section discusses the implementation details of PULSEDYN and may be skipped by a reader who wishes to use PULSEDYN as a black-box. Users who wish to adapt PULSEDYN to their specific needs can gather the implementation details of PULSEDYN in this section along with the documentation provided on the Github repository to make suitable changes. The code is organized as follows. First, the parameter file is parsed to extract values of the variables required for running the simulation. Any parameter not specified is simply set to its default value. If no parameter file is found, the code defaults to preset parameters given in the user manual. The user manual is provided with the code on the Github repository (see Appendix). Once the variables have been set, the objects for the classes in the software are initialized. After initialization, the model chosen is fed into the starter function for the simulation. Data is written to file periodically as explained in Section 3, until the simulation is completed. After the simulation is over the program exits. A flowchart illustrating the code organization is shown in Fig. 2. The code is written in an object oriented style to allow for good organization and legibility. The code contains 11 source files and 11 header files and the list is provided in the documentation file (in pdf) accompanying the software on the webpage provided in Table 3. The code is distributed as a Code::Blocks project fully compiled and accompanied by a stand alone executable. Any user with a basic knowledge of C ++ can modify and add more features to these source files and customize PULSEDYN. A Makefile is also provided to compile the source files if the user does not wish to rely on Code::Blocks to perform compilation. The Makefile may be used to generate executables using the make all command as well as clean up object files using the make clean command. The command make clean removes only the executable and the object files in the \bin and the \obj folders. There is also an option to use the make clean-dir command which deletes the \bin and \obj folders along with all their contents. It has been written keeping in mind the organization of the subdirectories of the source files of PULSEDYN. The difference between Code::Blocks and the Makefile compilation is that the executable is stored in the \bin sub-folder when the Makefile is used while Code::Blocks allows the user to generate the executables in either the \bin\Release or \bin\Debug sub-folders for the Debug and Release versions respectively. The Makefile is written in the Linux style. On Windows, we recommend using the Cygwin environment to issue the make command. Keep in mind, that when using Cygwin,

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

139

4.3. Potential classes There are four Potential classes corresponding to each of the four spring potentials provided in the code. They are todaPotential, fpuPotential, morsePotential and lennard-

JonesPotential. Note that the combinations of different kinds of springs cannot be made in the same chain. The information about the parameters in the potentials i.e., k1 , k2 and k3 , are stored by an object of the chosen potential’s class. The functions for force and potential energy corresponding to each of these models are written as member functions of the sub-classes corresponding to the potentials. An object of any of the Potential classes stores the following members, (a) parameters of the potential and (b) potential energy. Each of the potential classes contains the following functions, (a) getter and setter functions for each member, (b) acceleration function and (c) potential energy function.

Fig. 2. A high level flowchart for the software is shown in this figure.

the executables created are not standalone and the necessary dll files must be placed next to the executable if it is to be run by simply double clicking on it. These dll files are provided on the Github repository and may be downloaded for use. On Linux, the make utility comes standard and the make command can be issued at the terminal. PULSEDYN as it is currently written, requires a version of the GNU GCC compilers greater than 5.0. We request any user to give attribution to this work when reporting results using PULSEDYN or any edited version of the same. While PULSEDYN does not prefer any particular compiler or environment, we recommend using the GNU GCC compiler (Linux) or the MinGW compiler (Windows) to re-compile PULSEDYN if changes are made to the source files. The implementation also emphasizes encapsulation. For instance, a user need only modify the boundary condition files and update the input file interface in order to impose custom boundaries on the system, while leaving the rest of the code unchanged. The organization of the code is described next. 4.2. Particle class The Particle class, as the name suggests, defines a particle in a chain. A Particle object carries the following information about each particle in the chain, (a) (b) (c) (d) (e) (f)

position, velocity, acceleration, mass, kinetic energy and boundary information (useful only for boundary particles).

The following functions are included in the Particle class (a) getter and setter functions for each member and (b) kinetic energy calculator. For a chain of particles, a C + + vector of such objects is created and the dynamical quantities are updated in each iteration of the algorithm. Further details are given in Section 4.5. The organization of the class is shown in Fig. 3.

The acceleration function calculates the right hand side of Eq. (5) i.e., it excludes the driving and dissipation forces. When required, the function is able to call the boundary conditions functions (discussed later) to calculate the accelerations for the edge particles. The potential energy function calculates potential energy for each spring in the system. Similar to the acceleration function, the potential energy of the boundary springs is calculated by calling the external boundary conditions functions during the calculation. The information about which boundary conditions are to be used is stored in the Particle objects detailed before. The layout of the Potential classes is shown in Fig. 4. 4.4. External force While the Potential classes implement the potential energy and the resulting conservative force on the right hand side of Eq. (5), the Force class provides Fdiss and Fext as the velocity dependent damping force and a sinusoidal driving force on the right hand side of Eq. (6). The Force class members are (a) (b) (c) (d) (e) (f) (g) (h)

type of force (sine or cosine), start time, period of the periodic force, end time of the force, frequency (calculated from period above), ramp (for chirping), amplitude of driving force and velocity dependent damping γ .

The member functions of the Force class are (a) (b) (c) (d)

getter and setter functions for each member, force select function, sine driving force function and cosine driving force function.

The organization of the Force class is shown in Fig. 5. 4.5. Simulation class The Simulation Class contains the solvers used to integrate the equations of motion. The members of the class are as follows, (a) time step for the integration algorithms, (b) sampling interval,

140

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

Fig. 3. Shown here is the Particle class layout in PULSEDYN.

Fig. 4. Shown here is the layout for each of the Potential classes in PULSEDYN.

(c) total number of time snapshots of system evolution recorded, (d) size of the chain and (e) integration algorithm name. The functions implemented in this class are as follows, (a) (b) (c) (d)

getter and setter functions for each member, starter function for simulations, Gear 5th order predictor–corrector algorithm and Velocity-Verlet algorithm.

The solver functions in the class are implemented as function templates. Function templates in C++ are functions whose operations can be defined on objects of generic type. The template function performs the operations on the object passed to it, provided the class to which the object belongs contains the appropriate member functions with names identical to those the function template expects. Therefore, when an object belonging to any of the Potential sub-classes is passed through the solver function, the sub-class must have member functions defined with the same name: in this case the acceleration and potential energy functions. In PULSEDYN, the potential sub-classes already follow a common naming scheme and are compatible with the solvers and each other. Furthermore, the definitions of the member functions are added in their headers to allow the compiler to optimize code effectively. We get three advantages from this approach. First, it enables code reusability. Second, the runtime efficiency of the code is better

compared to the other options we can avail given the structure of this code. Third, it provides modularity to the code. However, if the user adds any new potential classes to the code, they must follow the same naming convention for the different member functions so that compatibility is maintained. The organization of the Simulation class is shown in Fig. 6. The starter function in the simulation class checks which integration algorithm is chosen in the parameter file and directs the code to the correct solver function. In each of the integration algorithms, the following routine is employed. First, the initial conditions are written to file. Then the algorithm variables needed to execute the update steps are initialized. After initialization, a time loop is started and the dynamical variables are updated in each iteration by the solver. During the update step, calculations of force and acceleration are performed to obtain the values of the position and velocity at the next iteration. The total force is calculated as a sum of the conservative and non-conservative parts of the right hand side (RHS) of Eq. (6). The conservative part is calculated by calling the acceleration functions inside Potential class and the non-conservative parts are calculated by calling the functions in the Force class. Inside the time loop, a counter is implemented which checks if the sampling interval is reached. This determines the writing condition. When the writing condition is satisfied in the counter, the energies are calculated and the data is written to file. Once file writing is complete, the counter for file writing is reset. The details of file output functions are given in Section 4.7. 4.6. Boundary conditions Boundary conditions are implemented as stand-alone functions in the software. The functions for boundary conditions return the value of the spring extension or compression at the right end of the Nth spring and the left end of the 1st spring. Three kinds of boundaries have been provided — open, fixed and periodic. The functions for boundary conditions are, (a) (b) (c) (d)

left boundary function for acceleration, right boundary function for acceleration, left boundary function for potential energy and right boundary function for potential energy.

The functions first evaluate the type of boundary chosen for both ends of the chain. Once the boundary type is evaluated, they return the spring extension or compression corresponding to the edge springs. In the case of the fixed boundary, the edge spring is fixed to a wall on one end and the boundary particle on the other. The functions then return the value of the spring extension

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

141

Fig. 5. Shown here is the layout for the Force class in PULSEDYN.

Fig. 6. Shown here is the layout for the Simulation class in PULSEDYN.

or compression by assuming that the wall displacement is zero at all times. If the boundary is open, the edge spring is connected only to the boundary particle i.e., the compression or extension in the boundary spring is always zero. In the case of periodic boundaries, there is only one boundary spring and it connects particles 1 and N. Therefore, the magnitude of change in the length of the boundary spring is |x1 − xN |. To avoid counting the single periodic boundary spring twice, the boundary spring energy is calculated only at the left chain edge while the energy of the periodic boundary spring is calculated to be zero at the right end. It is also worth mentioning that in PULSEDYN, while the open and fixed boundaries can be set independently on the ends of the chain, periodic boundaries cannot. If one of the boundaries is set to be periodic the other boundary is also forced to be periodic.

4.7. Output class

To write data, we have implemented a static class Output that can be called from anywhere in the code to write data to file. The following data is written to file, (a) (b) (c) (d) (e) (f)

position.dat for displacement data, velocity.dat for velocity data, acceleration.dat for acceleration data, ke.dat for kinetic energy data, pe.dat for potential energy data, restart.dat records displacement, velocity and accelera-

tion at the last time iteration, (g) mass.dat for mass values of each particle and (h) totalEnergy.dat for total energy at each recorded iteration.

142

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

Fig. 7. Shown here is the layout for the Output class in PULSEDYN.

The class contains the following members, (a) file names of data files and (b) file streams to data files. The functions in the class are the write functions to the files. Each of the functions take as arguments the value to be written to file and an end of line check. Once called, the functions open the associated files and write the values. These functions are called from the update functions in the Simulation class when the file writing counter determines that data writing condition is satisfied. The details of the format of data contained in the files is given in Section 3. The layout of the Output class is shown in Fig. 7. 5. Benchmarks PULSEDYN emphasizes legibility, organization and modularity over raw computational speed. That being said, accurate dynamical simulations of MBNSs are notoriously expensive computationally. So, it is necessary to get an estimate of the efficiency of the code. Therefore, a short discussion of PULSEDYN’s performance is provided here. The memory requirements of the code are not significant given the capacity of modern machines. As an example, a 5000 particle system requires around 2 MB of RAM. The typical scaling of running time of the program is shown in Fig. 8. The running time has been obtained for a system of 101 particles with the FPUT potential. The timer starts as soon as the program is launched and stops counting after the restart file is written. The code tested here was compiled with the GNU GCC compiler (version 7.3.0) on Linux. The testing was performed on a machine with an Intel Xeon E5-2630 CPU, 32 GB of RAM and a 7200 rpm hard drive running Ubuntu 14.04. The parameters are k1 = k2 = 0 and k3 = 1. The simulations end at the same point in system time for various values of dt and data is collected every 1/dt iterations. Therefore, as dt is reduced, file writing becomes less frequent. The system time at the end of each of the simulations is 2000 (arbitrary

Fig. 8. Shown here is the running time (in seconds) of PULSEDYN as a function of the time step dt for the velocity-Verlet and Gear 5th order predictor–corrector algorithms. These running times are obtained for the FPUT chain with 101 particles and single initial velocity perturbation.

units). The initial conditions are a single velocity perturbation to a particle at the center of the chain. There are two extremes to be noted here. Both the Gear algorithm and the velocity-Verlet algorithm show similar running times when dt 10−1 − 10−3 . This is due to the fact that the main bottleneck for dt in this range of values is file writing. As dt is decreased to less than 10−4 , velocity-Verlet is observed to be faster as expected since there are fewer update steps in the velocityVerlet algorithm. However, the difference is not significant unless late time simulations are performed. This is due to the fact that the primary computational costs of the algorithms are in the calculation of the inter-particle force terms. We do not perform late time simulations using the Gear 5th order predictor–corrector algorithm due to its inherent issues with energy drifting. Addition of external driving and dissipative forces also does not add significant load on the performance of the Gear 5th order predictor–corrector algorithm.

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

143

using various combinations of the boundary conditions, solvers and initial conditions. 6.1. Toda lattice In the limit N → ∞, the Toda lattice from Eq. (1) admits a soliton solution of the following form [31], xn =

Fig. 9. The logarithm of the root mean squared fluctuations in the total energy as a function of logarithm of time-step dt is shown in this figure for the velocity-Verlet and Gear 5th order predictor–corrector algorithms. We have used the FPUT chain with k1 = k2 = 0, k3 = 1 and a velocity perturbation given to a single particle in the chain.

The code used for testing was compiled with the -O3 flag which allows the compiler to perform SIMD vectorization, which streamlines the calculations. Vectorization is successfully performed on loops for the position and velocity update steps (excluding the force calculation) which would normally be carried out serially on a processor. The compiler flag -fopt-info-vec-optimized gives information about which parts of the code were successfully vectorized. However, we have found the speed up to be negligible for our systems compared to the -O2 level optimization which does not perform the vectorization (data not shown here). There are two possible reasons for this. First, as mentioned above, calculation of the interaction terms is the dominant computational cost assuming file writing is infrequent enough. Therefore, even after vectorization is successfully performed on the update steps, the speed up is not noticeable. Second, for SIMD vectorization to be a significant enough improvement, some amount of finetuning of the code for the processor being used may be necessary. Since, performing the fine tuning on our test machine may reduce the portability of the code, we leave it up to the user to implement such a step based upon the calculational needs. Higher optimization levels such as -Ofast are not used here since they typically implement the -ffast-math directive which is not strictly standard compliant. However, depending upon the user’s needs, further optimization can be done by making the necessary changes to the code and then adding the appropriate flags to the Makefile provided with the code. Alternatively, the modified files may also be saved as part of the Code::Blocks project file and the project file may be rebuilt in Code::Blocks with the appropriate flags added to the compiler options. We use the accuracy associated with conservation of energy as a way to test the accuracy of the algorithms. Fig. 9 shows how the rms value of the total energy fluctuations in the system depends on the time step used for the integration. We expect that the velocityVerlet algorithm roughly scales as O(dt 2 ) while the Gear 5th order predictor–corrector scales as O(dt 5 ) [30]. To verify this behavior, we have overlaid the fits of the form log(dERMS ) = a log(dt) + b

(7)

to the data obtained from the velocity-Verlet and the Gear 5th order predictor–corrector algorithm. From our fits, we find a = 2.01 ± 0.01 and a = 5 ± 1.1 for velocity-Verlet and the Gear 5th order predictor–corrector algorithms, respectively, which match our expectations. 6. Results: test cases In this section, we discuss results that we have reproduced from literature for each of the four systems in the PULSEDYN code

1 k2

[ ln

1 + exp(2(κ n − κ + β t)) 1 + exp(2(κ n + β t))

] ,

(8)

where xn refers to the displacement of the nth mass. The velocity of the nth particle can be found by taking a derivative with time of xn and is

vn = −

2β exp(2(κ n + 2β t))(exp(2κ ) − 1) k2 (1 + exp(2(κ n + β t)))(exp(2κ ) + exp(2(κ n + β t)))

. (9)

In the above solution, κ is a free parameter which √ controls the width of the soliton, width ∝ 1/κ and β = ±( k1 k2 /m sinh κ ). Here, κ > 1. Solitons of smaller widths propagate faster and the sign of β controls the direction of soliton propagation. For our simulations, we used k1 = 1.0, k2 = 10.0 and varied κ to seed soliton profiles into the Toda lattice. We use Eqs. (8) and (9) to seed the soliton profile into the Toda chain. The Toda lattice soliton solution holds true for an infinite system. Due to the fact that the soliton is kink shaped i.e., the right end is in a different configuration than the left end, the finite sized chain recoils as soon as the system is released. This leads to some radiative emissions from the soliton profile seeded into the chain (too small to be seen in Fig. 10, but visible in Fig. 11). In our simulations, x1 = 0 and xN = −0.06. Due to emission of this radiation, the soliton slows down slightly as compared to the theoretical solution. The effects of boundaries on the system dynamics is an important issue and has been the subject of previous study where the slowing down of the numerical solutions as compared to the theoretical solution has been reported [44]. In our simulations, we have used open boundaries and seeded the kink closer to the right boundary to minimize the recoil. The time-step used is 0.01 and the data is recorded at every 100 iterations. The results are shown in Fig. 10. The energy is conserved to one part in 106 using a Gear 5th order predictor–corrector algorithm. The Toda lattice also allows for multiple soliton solutions. We show in Fig. 11 two solitons seeded in the chain at opposite ends and traveling towards each other. As soon as the simulation starts, the solitons radiate some energy as stated above. After the radiation separates from the solitons, they move with constant speed towards each other and meet at the center of the chain. As expected, the two solitons show no scattering when they collide with each other. Further, the radiation moves slightly slower than the solitons themselves. This is due to the fact that the soliton speed in this system is always greater than the speed of sound which is the fastest speed at which the radiation can propagate [3,4,31]. For the soliton collision problem we have used κ = 0.6, k1 = 1.0 and k2 = 10.0 for both solitons. 6.2. FPUT Lattice The FPUT system is a nonlinear spring–mass model with a polynomial potential given by Eq. (2). The model was originally proposed and studied to show that in contrast to the finite sized harmonic oscillator chain, a weakly nonlinear spring–mass chain would thermalize i.e., the energy of the system would be shared between different modes and this would lead to the eventual equipartitioning of energy [1,2]. What they found however, was

144

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

Fig. 10. We show the Toda soliton propagating on a lattice (moving to the left) at times t = 0, 10, 40 in and 100 in (a) through (d). The numerical solution is shown as the solid line while the exact solution is shown as a dashed line.

Fig. 11. The collision of two Toda solitons as they propagate on the lattice is shown in the panel. The figure shows the logarithm of kinetic energy plotted against time on the x-axis and particle index on the y-axis. The shade of gray sets the scale of the energy plotted. Darker shades correspond to higher energy as shown in the color bar.

something they never expected. Their simulations showed that the energy sharing between modes is at best very weak. The system showed no observable equipartitioning of energy in the time window across which it was probed. Using PULSEDYN, we recreated the recurrence phenomenon observed by Fermi and his collaborators [1,2]. For our simulations, we used k1 = 0.5, k2 = 0.0833 and k3 = 0 in a 64 particle system as described by Eq. (2) with mi = 1. We initialized the system with all the energy fed into the lowest mode, l = 1 using the following equation, xj =

N ∑

x(l) sin

l=1

[ jlπ ] . N +1

(10)

The mode energies are then given by El =

1 2

a˙ 2l + ωl2 a2l .

(11)

Here,

ωl2 =

4l m

sin2

[

lπ 2(N + 1)

]

,

(12)

Fig. 12. The energy stored in each of the first 4 modes of the α -FPUT chain as a function of time is shown here. The recurrence of energy in the first mode is seen at t ≈ 60000. Here, k1 = 0.5, k2 = 0.0833, k3 = 0 and the energy is initialized in the first mode of the system i.e. l = 1.

and al =

N 2 ∑

N

j=1

xj sin

[ jlπ ] . N +1

(13)

While the energy is initially shared with higher modes l = 2, 3, . . ., it returns back to the original mode after some time as shown in Fig. 12. The total energy here is conserved to 1 part in 108 . We have used a time step of 0.001 with the velocity-Verlet algorithm [28]. We have used fixed boundaries for the system. We now turn our attention to phenomena in the FPUT chains which have been studied more recently. It has been shown that a bond squeezing or stretching in the β model of the FPUT chain i.e. k2 = 0 leads to a localized nonlinear excitation (LNE) [21,32]. A squeeze in a bond is effected by setting the following initial condition, xi = −xi+1 = A. The LNE is not stable and it delocalizes after some time. The time scale depends on the parameters of the system as well as A [32]. In the absence of phonons, i.e. k1 = 0, it has been also shown that the destabilization takes place via compression pulses, which we call SWs, dilation pulses, which we call anti-SWs (ASWs) and other metastable nonlinear excitations [32].

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

145

Fig. 13. A contour plot of log10 KE in the quartic FPUT system with k1 = 0 is shown here. The LNE is in bond 51 of a 100 particle chain with fixed boundaries. The color sets the energy scale as shown in the color bar.

Following the destabilization of the LNE, the system then enters a state called quasi-equilibrium (QEQ) where there is no memory of initial conditions, there is a Gaussian distribution of velocities but no equipartitioning of energy [20,22,45,46]. The potential is given by k1 = k2 = 0 and k3 > 0. A simulation of a seeded LNE in the FPUT system is shown in Fig. 13. The LNE initially emits weak SWs and ASWs followed by other metastable excitations. A detailed discussion of the LNE decay and the evolution of the system to the energy equipartitioned state is presented in Section 7. The time step used here is 10−6 . A small time step is required for a problem such as this one owing to the high frequency quasi-periodic oscillations of the LNE before delocalization. The energy is conserved to better than one part in 1011 using the Gear 5th order corrector–predictor algorithm [32]. Another test for the quartic FPUT system is the dynamics of SW collisions in the system. A SW can be seeded in the chain by giving a velocity perturbation to a particle at t = 0 and letting it evolve. We seed two SWs in a quartic system (with k1 = k2 = 0). The SWs are identical and propagate towards each other from opposite sides of the chain. It has been shown that two such interacting waves undergo scattering when they collide and they leave behind residual energy [33]. The results are shown below for a 2000 particle chain in Fig. 14. Two SWs were created at particle 400 and 1600 by giving them velocity perturbations vo = 0.3 and −0.3, respectively. The boundaries are set to be free boundaries and the simulation was performed with the velocity-Verlet algorithm with a time step of 10−3 . The total energy of the system is conserved to 1 part in 107 . At particle 1000 where the two SWs meet they collide and leave behind weak oscillations as expected. Due to energy conservation, the heights of the SWs post collision are slightly reduced. If this process happens repeatedly, the SWs break up into many small SWs which we have called secondary SWs (SSWs) and the system enters into the QEQ phase. 6.3. Morse and Lennard-Jones potentials Exact solutions to the equations of motion for Morse and the Lennard-Jones (LJ) chains do not exist to our knowledge. However, numerical studies have shown evidence of SW propagation in these lattices (for instance, see Ref. [47,48]). Therefore, to test propagation of SWs through chains with the Morse and LJ potentials, we follow the numerical scheme employed by Flytzanis et al. [25,26]. In their 1989 study, the authors reduce the equations of each of the potentials under consideration into their corresponding α + β FPUT forms as truncated polynomial expansions with V (0) = 0. They show that in the continuum limit, these potentials can be approximated by the generalized Boussinesq equation which has well known soliton solution. We use the solution as a generator to

Fig. 14. In (a) we show two SWs in the quartic FPUT system traveling towards each other at t = 1150. In (b) we show the two SWs moving away from each other towards the ends of the chain at t = 1560 after colliding with each other. The collision takes place at particle 1000. Notice that the collision is inelastic and the SWs lose energy. The formation of two SSWs is also seen in (b) at particles 970 and 1030. More low energy SWs would be formed at the region centered around particle 1000.

seed a SW into the chains. Please note that Flytzanis et al. refer to these excitations as solitons in their study. However, we refer to the excitations created as SWs keeping in mind our original definition. Our calculations using PULSEDYN can be described as follows. Consider the Morse potential given by Eq. (3). The potential can be expanded and written as a truncated polynomial potential of the form Vi+1,i =

1 2

G(xi+1 − xi )2 +

1 3

A(xi+1 − xi )3 +

1 4

B(xi+1 − xi )4 .

(14)

In the continuum limit, the RHS of Eq. (14) can be written as V ′ (ri+1 ) − 2V ′ (ri ) + V ′ (ri−1 ) ≈ D2

∂2 V

∂ y2 +

1 12

D4

∂ 4V . ∂ y4

(15)

Here, ri = xi+1 − xi and y = iD where D is the bond length of each bond and i is the site index. The corresponding Boussinesq equation takes the form utt − co2 uyy − p(u2 )yy − q(u3 )yy − huyyyy = 0,

(16)

with = GD /m, p = AD /m, h = GD /(12m) and q = BD /m, where G > 0, B > 0, and A is real and positive, negative or 0. The soliton solution of this equation is given by co2

2

y(x, t) = ±2sgn(h)

3

( 2h )(1/2) q

4

arctan

(1 w

4

( x − v t ))

tanh

L

,

(17)

where

( w=

) [4p2 + 18(v 2 − co2 )q]1/2 ± 2p , [4p2 + 18(v 2 − co2 )q]1/2 ∓ 2p

(18)

146

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

and

[

L=2

]1/2

h

v 2 − co2

.

(19)

Since this solution is not for the discrete system, it decomposes into a radiative part and a kink shaped SW. In the study [26], the authors found that the propagating kink SW does not fit well to polynomial and other simple functions. In fact, the stable part of the propagating kink is fitted best to a generalized Toda soliton form given by the equation xn =

1 2

Aw m log

1 + exp[±2(n − no − 0.5)/w] 1 + exp[±2(n − no + 0.5)/w]

+ C,

(20)

where, Am is the height of the kink, no is the position of the center of the kink and 2w is the width of the kink soliton. Flytzanis et al. [26] were able to obtain fits of the non-radiative part of the propagating kink in the LJ and the Morse chains using the above forms for chains with N ≫ 1. The similarity of the numerical solitons to the Toda soliton was also observed and reported elsewhere, prior to the study by Flytzanis et al. (for instance, see Ref. [49,50]). We now use the same approach to seed SWs in the Morse potential. We then verify that the propagating constant velocity front of the kink can indeed be approximated well by the generalized Toda soliton solution in Eq. (20). However, we note that while the approximate solution works well, the continuum approximation becomes inaccurate when the width of the soliton is comparable to the lattice spacing. For the Morse potential, we used the parameters k1 = 0.01 and k2 = 7. This leads to the values of G = 1, A = −10.5 and B = 57. Then we use the values of p, q, h and co2 to calculate the generator from Eq. (17) and seed this shape into the chain. We track the propagating kink and fit it to the generalized Toda soliton solution from Eq. (20). The time step used is 0.1 and the energy conservation holds to 1 part in 105 using the velocity-Verlet algorithm. The plot is shown in Fig. 15(a). As with Flytzanis’ study, the fit obtained confirms that PULSEDYN is able to reproduce Flytzanis et al.’s well known results for the Morse chain. We calculate the parameters of the fit as Am = 0.117 ± 0.005, w = 1.00 ± 0.02 and no = 368. The goodness of the fit R2 = 0.999. It must be noted that the fit was obtained after the soliton suffers multiple collisions with the boundaries of the system and has undergone considerable scattering. However, for the LJ potential, we have taken a different approach. We consider open boundaries and give a velocity perturbation to the first particle in the chain. This approach is similar to the one used for generating SWs in the FPUT system [25]. The rationale is that the velocity perturbation also generates a kink shaped SW. Therefore, we anticipate that the generalized Toda solution provided by Eq. (17) still provides a good approximation to the SW generated through this means. We use k2 = 0.01389 and k1 = 1 and set D = 1 for the simulations. A velocity perturbation of v0 = 0.2 is given to particle 1 and the boundaries are kept open. Our simulations confirm that the generalized Toda solution provides a good description for the propagating SW front in the LJ chain. For the case of the LJ chain, we calculate, Am = 0.09958 ± 0.000015, w = 0.6481 ± 0.00015, no = 287.57 ± 0.05 and R2 = 0.999. 7. Results: Route to equilibrium We turn our attention now to a problem of historical significance. We consider the FPUT chain, presented previously in Section 6.2. We take a different approach from the authors of the original FPUT system by focusing on the strongly nonlinear FPUT system. In particular, we study the purely nonlinear FPUT chain with a seeded LNE. The reasons for this approach are as follows.

Fig. 15. We show a snapshot of the kink profiles for the Morse chain at t = 310 and for the LJ chain at t = 54 in (a) and (b) respectively. The data points are obtained from numerical simulations and the curve shows the fit to Eq. (20).

While the traditional FPUT system with weak nonlinearity has been the subject of many studies, the strongly nonlinear regimes have not been sufficiently probed. This is due to the lack of analytical techniques available for these systems. Hence, we believe that simulations are of utmost importance for studying strongly nonlinear systems. These systems also bear a strong resemblance to granular systems such as those described in Ref. [23]. Phenomena such as formation of SWs, interaction between them and the formation of the quasi-equilibrium (QEQ) phase (described below) have been reported in both granular systems as well as 1D spring–mass systems with strong nonlinearity [21,35,45]. Further, in the purely nonlinear regime these systems show qualitatively different behavior [45,51]. For instance, phenomena such as rogue fluctuations, which are discrete analogs of rogue waves found in oceans [52,53] and unstable nonlinear objects [32] are observed in these systems. However, such phenomena can be suppressed except at very late times due to the presence of interactions between the SWs, the unstable nonlinear objects and the phonons. As demonstrated in the previous section, a seeded LNE eventually delocalizes. As it delocalizes, it does so in a way that as the system evolves sufficiently far into time, its average kinetic and potential energies approach the values predicted by the virial theorem. At early times in the LNE evolution, the motion is very nearly periodic with the frequencies corresponding to those of the unforced and undamped Duffing oscillator [32]. As it evolves and emits energy via SWs and other excitations, the frequency peaks characterizing the dynamics of the LNE broaden as shown in Fig. 16, i.e., the LNE’s periodicities no longer remain well defined. This frequency broadening further accelerates the delocalization process of the LNE. Once delocalization is complete the system enters a state where we find a Gaussian distribution of velocities of

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

147

Fig. 16. Shown here are the direct cosine transforms (DCTs) of the kinetic energy of the LNE particle (particle 50) at two different time intervals in its evolution. In (a) we show the DCT for the time interval t = 0 to t = 100 and in (b) we show the DCT computed for the time interval from t = 2500 to 2600. The parameter values are k1 = 0, k2 = 0, k3 = 1.0 and A = 0.30. At early times, the DCTs show sharper peaks. With time, as the LNE delocalizes, the DCT shows frequency spreading i.e. the LNE is not quite periodic anymore.

the particles, no dependence on initial conditions and no equipartitioning of energy. We have introduced this state as the QEQ state in earlier work [45]. For more recent work on QEQ the reader is referred to Refs. [20,21,34,46]. As the LNE delocalizes, the kinetic energy fluctuations δK of the system begin to decrease since the motion of the LNE particles is not nearly as periodic as at early times. A plot of δK fluctuations during the LNE lifetime is shown in Fig. 17(a). With time, the fluctuations go down further and reach a stable value as shown in Fig. 17(b). This regime is well after the delocalization of the LNE. It is natural to ask whether the system remains in the QEQ phase indefinitely or whether it attains an equipartitioned state. In true equilibrium, ergodicity would hold true and the energy would be equipartitioned. We know that the answer to this question is known for the FPUT system in the weakly nonlinear regime, i.e., the FPUT chain goes to the energy equipartitioned state at late times for weak nonlinearity [36–38]. Our objective hence is to explore whether PULSEDYN can help settle whether the system in the fully nonlinear case (recall k1 = k2 = 0) attains the equipartitoned state. To address this question, we performed calculations to verify the virial theorem and the specific heat of the system. The parameters used are k1 = 0, k2 = 0, k3 = 1.0 and A = 0.3. These are the same parameters as those for the LNE in Fig. 13. We use a time step of 10−5 after LNE delocalization and run the simulation for 1012 time iterations using the velocity-Verlet algorithm. We sample data at every 105 iterations and we have called each such sampling a single time step. We have recorded data for 1.4 × 107 time steps.

Fig. 17. In (a) and (b) we show the kinetic energy fluctuations as a function of time at early time (before LNE delocalization) and at late time respectively. The time averaged and normalized kinetic energy is shown as a function of time deep into QEQ in (c). At late times, notice that the fluctuations stabilize and the kinetic energy average matches the virial theorem prediction.

We find that at late times – the last 15% of the time evolution – the space average and the time average of K are equal to each other, i.e., the system is ergodic. We then check whether the virial theorem results are satisfied. The virial theorem states n ⟨K ⟩ = E. (21) n+2 Here, E is the total energy of the system and n is the exponent in the power law potential. For the purely nonlinear β -FPUT chain, n = 4. Therefore, ⟨K ⟩ = 23 E at late times. During LNE evolution, the kinetic energy of the chain fluctuates with high frequency [32]. After LNE delocalization, the average kinetic energy reaches the value predicted by the virial theorem as shown in Fig. 17(c). However, this is not sufficient evidence for energy equipartitioning since ergodicity and the virial theorem can be both satisfied even in the QEQ state [46]. Further, the velocity distribution in the QEQ state has been found to be Gaussian [20,35,45]. Therefore, for the final state obtained to be truly in equilibrium, the system has to satisfy energy equipartitioning. We test this, by calculating the specific heat, Cv , of the system at late times.

148

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

Calculations of response functions in the microcanonical ensemble of a nonlinear many-body system are challenging [54,55]. To examine whether the FPUT system has reached the equipartitioned state, we rely on a recent similar study carried out for an alignment of grains where the grains repel via a strongly nonlinear algebraic potential [56]. There it has been shown that

[ Cv ≈ kB

n+2

−

2n

1 N

(

n+2 n

+

4(N − 2) 2nN

)] ,

(22)

where, N is the number of particles in the finite sized system. In the above equation, if N → ∞, we recover the result from Tolman’s generalized equipartition theorem due to the equivalence of statistical ensembles in the thermodynamic limit [57] as

( Cv ≈

n+2

)

2n

kB .

(23)

Since the result depends only on the exponent in the potential n, Eq. (22) holds for the purely nonlinear β -FPUT chain when the energy is in the equipartitioned state. For N = 100, the theoretical value for Cv /kB from Eq. (22) is calculated to be 0.7252. A relation also has been established between the Cv of a system in the microcanonical ensemble and the kinetic energy fluctuations in [56,58], Cv =

2kB N

( 1−

(N − 8)⟨1/K 2 ⟩ (N − 4)⟨1/K ⟩2

)−1 .

(24)

From our numerical calculations, we recover a value of 0.7287 at late times using Eq. (24). The numerical value from our simulations is within 0.5% of the theoretical value. The numerical error in the calculation of energies is of the order of 10−9 . Using the methods outlined in Ref. [59] for error propagation, we estimate that the error in Cv is of the order of 10−5 . Therefore, the calculation of Cv is safe from numerical errors up to the fourth decimal place. It must be noted however, that this is an estimate for how numerical errors might affect Cv calculation. The fluctuations in the system kinetic energy themselves will result a fluctuation of Cv estimated from Eq. (24), depending on the time interval chosen for averaging. It is also important to note here that the numerical error is not the same as the fluctuation in calculation of Cv which might arise due to the fluctuations in K about ⟨K ⟩. These fluctuations may cause the calculated value of Cv to show small changes as Nt is reduced. On the other hand, using a large Nt must also be avoided since large values of Nt would cause the averaging interval to overlap with time intervals of evolution of the system in which energy is not yet equipartitioned. For our calculations, we have set Nt = 1 × 106 time steps and we recover reliable results. Therefore, our studies using PULSEDYN provide evidence that the purely nonlinear β -FPUT system goes to the equipartitioned state at late times.

license. This code in its current form can be used to study the Toda chain (an integrable system), and the FPUT, Morse and LennardJones chains (non-integrable systems). Conservative systems can be solved accurately using the velocity Verlet algorithm and the Gear algorithm whereas the driven dissipative system can only be studied using the Gear algorithm. The program has been written with a modular design and the programming style is expected to allow for additions and enhancements to the code by an experienced C++ programmer. To demonstrate the capabilities of PULSEDYN, we have recovered the following results. First, we recovered the exact analytical results on the propagation of a soliton and the interaction between solitons in the Toda chain. Next, we considered the FPUT chain and confirmed the presence of the recurrence phenomenon in the weakly nonlinear regime as studied originally by Fermi et al. [1]. Thereafter we considered the FPUT chain in the less explored strongly nonlinear regime and demonstrated that our results for LNE delocalization and SW interaction are in agreement with known results [32]. We then recovered the approximate forms of SWs in the Morse and the LJ potentials that have been previously investigated [25,26]. Finally, we presented new results regarding the relaxation processes of LNEs in the purely nonlinear FPUT chain. We demonstrated that as the LNE delocalizes, the fluctuations in the system’s kinetic energy K decrease and stabilize at late enough times as expected. At sufficiently late times the energy distribution reaches the values predicted by the virial theorem and the system is found to reach the energy equipartitioned state. The equipartitioning is evident through the excellent agreement between the specific heat calculated from the simulations and the same predicted theoretically. In closing we note that recent work has shown that the dynamics of nonlinear systems can be highly sensitive to the competition between linear and nonlinear forces. While there are analytical tools that can sometimes be used to study weakly nonlinear systems, very little analytical work is possible for most strongly nonlinear systems. And even less can typically be calculated when the linear and nonlinear forces become comparable [19]. We hope that PULSEDYN will help stimulate more studies on the dynamics of strongly nonlinear chains where, perhaps, fundamentally important processes remain to be discovered. Acknowledgments We thank Tyler Barrett, Guo Wen, Kevin VanSlyke, Nathaniel Fuller and Alexandra Westley for their help in testing the PULSEDYN code for various cases. Rohith Kashyap’s advice during the development of the code is gratefully acknowledged.

8. Conclusions We have presented the detailed architecture of the PULSEDYN code and how it can be openly used by anyone to carry out highly accurate dynamical simulations of discrete, nonlinear systems in one dimension. The design is such that both conservative as well as driven, dissipative systems can be studied using PULSEDYN. The program can be used by someone with limited or no coding experience. The discussions presented address the details of how to do long-time simulations of these systems with strict control on round-off errors. Further, the code provides features that can be readily tuned using a parameter file interface. PULSEDYN is distributed as an executable for both Windows and linux platforms and is distributed under the GNU GPL v3

Appendix. Required metadata

A.1. Current executable software version See Table 3. A.2. Current code version See Table 4.

R. Kashyap and S. Sen / Computer Physics Communications 239 (2019) 134–149

149

Table 3 Software metadata. Nr.

(executable) Software metadata description

Please fill in this column

S1

Current software version

v1.1

S2

Permanent link to executables of this version

https://github.com/rahulkashyap7557

S3

Legal Software License

GNU General Public License 3 (GPL)

S4

Computing platform/Operating System

Windows, Linux

S5

Installation requirements

none

S6

If available, link to user manual — if formally published include a reference to the publication in the reference list

https://github.com/rahulkashyap7557

S7

Support email for questions

[email protected]

Table 4 Code metadata. Nr.

Code metadata description

Please fill in this column

C1

Current code version

v1.1

C2

Permanent link to code/repository used of this code version

https://github.com/rahulkashyap7557

C3

Legal Code License

GNU General Public License 3 (GPL)

C4

Code versioning system used

Github

C5

Software code languages, tools, and services used

C++

C6

Compilation requirements, operating environments & dependencies

GNU GCC compiler (Linux) MinGW compiler or Cygwin GCC compiler (Windows)

C7

If available Link to developer documentation/manual

https://github.com/rahulkashyap7557

C8

Support email for questions

[email protected]

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]

E. Fermi, J. Pasta, S. Ulam, Los Alamos Report la-1940, 1955. E. Fermi, Collected Papers, vol. II, Univ. of Chicago Press, 1965. M. Toda, J. Phys. Soc. Japan 22 (1967) 431–436. M. Toda, J. Phys. Soc. Japan 23 (1967) 501–506. T. Schneider, E. Stoll, 1978 Proc. of the Symp. on Nonlinear (Soliton) Structure and Dynamics in Condensed Matter, vol. 1, Springer-Verlag, Berlin, 1978. S. Flach, A.V. Gorbach, Phys. Rep. 467 (2008) 1–116. M.J. Ablowitz, B.M. Herbst, C.M. Schrober, Physica D 87 (1995) 37–47. G.P. Berman, F.M. Izrailev, Chaos 15 (2005) 015104. D.K. Campbell, P. Rosenau, G.M. Zaslavksy, Chaos 15 (2005) 015101. J. Ford, Phys. Rep. 213 (1992) 271–310. N. Saito, H. Hirooka, J. Phys. Soc. Japan 23 (1967) 167–171. G.E. Duvall, R. Manvi, S.C. Lowell, J. Appl. Phys. 40 (1969) 3771–3775. R. Manvi, G.E. Duvall, S.C. Lowell, Int. J. Mech. Sci. 11 (1969) 1–8. R. Reigada, A. Sarmiento, K. Lindenberg, Phys. Rev. E 64 (2001) 066608. R. Reigada, A. Sarmiento, K. Lindenberg, Physica A 305 (2002) 467–485. A.J. Sievers, S. Takeno, Phys. Rev. Lett. 61 (1988) 970–973. M. Sato, A.J. Sievers, Nature 432 (2004) 486–488. G. Bennetin, A. Ponno, J. Stat. Phys. 144 (2011) 793–812. Y. Takato, S. Sen, Europhys. Lett. 100 (2012) 24003. E. Avalos, D. Sun, R.L. Doney, S. Sen, Phys. Rev. E 84 (2011) 046610. S. Sen, T.R. Krishna Mohan, Phys. Rev. E 79 (2009) 036603. S. Sen, J. Hong, J. Bang, E. Avalos, K.L. Doney, Phys. Rep. 462 (2008) 21–63. V. Nesterenko, Dynamics of Heterogenous Materials, Springer-Verlag, New York, 2001. G. Eilenberger, Solitons: Mathematical Methods for Physicists, SpringerVerlag, Berlin, 1981. T.P. Valkering, C. de Lange, J. Phys. A: Math. Gen. 13 (1980) 1607–1621. N. Flytzanis, St. Pnevmatikos, M. Peyrard, J. Phys. A: Math. Gen. 22 (1989) 783–801. C.W. Gear, The Numerical Integration of Ordinary Differential Equations of Various Orders, Technical Report 7126, Argonne National Laboratory, 1966. W.C. Swope, H.C. Anderson, P.H. Berens, K.R. Wilson, J. Chem. Phys. 76 (1982) 637. L. Verlet, Phys. Rev. 159 (1967) 98–103.

[30] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon, Oxford, 1987. [31] M. Toda, Theory of Nonlinear Lattices, Springer-Verlag Berlin, Heidelberg, Berlin, 1981. [32] R. Kashyap, A. Westley, A. Datta, S. Sen, Internat. J. Modern Phys. B 31 (2017) 1742014. [33] H. Zhao, Z. Wen, Y. Zhang, D. Zheng, Phys. Rev. Lett. 94 (2005) 025507. [34] E. Ávalos, S. Sen, Phys. Rev. E 79 (2009) 046607. [35] T.R.K. Mohan, S. Sen, Pramana 64 (2005) 423–431. [36] R. Livi, M. Petini, S. Ruffo, M. Sparpaglione, A. Vulpiani, Phys. Rev. A 31 (1985) 1039. [37] M. Onorato, L. Vozella, D. Proment, Y.V. Lvov, Proc. Natl. Acad. Sci. 112 (2015) 4208. [38] G. Benettin, H. Christodoulidi, A. Ponno, J. Stat. Phys. 152 (2013) 195. [39] W.E. Ferguson, H. Flaschka, D.W. McLaughlin, J. Comput. Phys. 45 (1982) 157. [40] M. Hénon, Phys. Rev. B 9 (1974) 1921. [41] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, New York, 1985. [42] P.M. Morse, Phys. Rev. 34 (1929) 57–64. [43] J.E. Lennard-Jones, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 106 (1924) 463–477. [44] Y. Shen, P.G. Kevrekidis, S. Sen, A. Hoffman, Phys. Rev. E 90 (2014) 022905. [45] S. Sen, T.K. Mohan, J.M. Pfannes, Physica A 342 (2004) 336–343. [46] E. Ávalos, S. Sen, Phys. Rev. E 89 (2014) 053202. [47] J.H. Batteh, J.D. Powell, J. Appl. Phys. 49 (1978) 3933–3940. [48] D.E. Strenzwilk, J. Appl. Phys. 50 (1979) 6767. [49] J. Dancz, S.A. Rice, J. Chem. Phys. 67 (1977) 1418–1426. [50] T.J. Rolfe, S.A. Rice, J. Dancz, J. Chem. Phys. 70 (1979) 26–33. [51] R. Reigada, A. Sarmiento, K. Lindenberg, Chaos 13 (2003) 646–656. [52] D. Han, M. Westley, S. Sen, Phys. Rev. E 90 (2014) 032904. [53] R. Kashyap, S. Sen, Phys. Rev. E (2018) submitted for publication. [54] E. Scalas, A.T. Gabriel, E. Martin, G. Germano, Phys. Rev. E 92 (2015) 022140. [55] J.R. Ray, H.W. Graben, Phys. Rev. A 44 (1991) 6905–6908. [56] M. Przedborski, S. Sen, T.A. Harroun, Phys. Rev. E 95 (2017) 032903. [57] R.C. Tolman, Phys. Rev. 11 (1918) 261–275. [58] H.H. Rugh, J. Phys. A Math. Gen. 31 (1998) 7761–7770. [59] J.R. Taylor, An Introduction to Error Analysis: The Study of Uncertainties in Physics Measurements, Second Ed., University Science Books, Sausalito, 1997.