This document describes a method for the creation of thermally equilibrated molecular configurations for large scale molecular dynamics simulations and its parallel implementation. The method combines relaxation of a stochastic initial configuration by Monte Carlo moves followed by molecular dynamics. Benchmarks are presented for up to 800000 particles

Hloucha, Matthias

Juelich Shared Electronic Resources

FORSCHUNGSZENTRUM JÜLICH GmbHZentralinstitut für Angewandte MathematikD-52425 Jülich, Tel. (02461) 61-6402Interner BerichtThermalization of Large Fluid Systemsconsisting of Lennard-Jones MixturesMatthias HlouchaFZJ-ZAM-IB-9902Februar 1999(letzte Änderung: 4.2.1999)Contents1 Program MCMD 11.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Usage 42.1 Compilation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2 Running Jobs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.3 Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.4 Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.5 Recipe for building a large thermalized configuration . . . . . . . . . . . . . . . . 53 Results 73.1 Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73.2 Benchmarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73.3 Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11iii CONTENTSChapter 1Program MCMDOld Chemists never die – they just reach equilibrium1.1 IntroductionAt the Center for Applied Mathematics, Research Center Ju¨lich, Germany, work has been started toinvestigate massively parallel algorithms for molecular dynamics computer simulations of liquids.As part of this project, a program MCMD was developed in order to construct starting configurationsfor simulations with a large number of particles. The ideas, usage and some benchmarks of MCMDare presented by the present text.There are two common methods to construct a starting point for a subsequent liquid simulation.While some authors suggest to start from a regular lattice [1], others prefer random configurations[2]. Both methods have their advantages and drawbacks. A start from a regular lattice might beproblematic for dipolar fluids at high densities, where melting is inhibited by long range dipolarinteractions. This approach is also expected to fail for mixtures of molecules of significantly dif-ferent size. Random mixing of molecules, however, often leads to unrealistic short intermoleculardistances. The resulting huge repulsive forces would make molecular dynamics simulations eitherinstable or one would have to use a very short time step.The idea of this work is to use a combination of both methods – random mixing and lattice start –together with Monte Carlo and molecular dynamics simulation techniques.1.2 MethodThe principal function of MCMD is illustrated in the Figures 1.1–1.4. First, the required particlenumber and volume is divided into a subsystem, which is allowed to contain not more than 1,000particles. The subsystem is also chosen to be a unit cell of a cubic lattice for the whole system.If it is not possible to divide the particle number exactly into a cubic subsystem, then the particlenumber in the subsystem is slightly increased.The coordinates of particles in the subsystem are randomly distributed in the given volume (seeFigure 1.1). In the picture of an Argon–Krypton mixture at a typical liquid density some unrealisticshort interatomic distances are clearly visible as overlapping spheres.A Monte Carlo simulation in the NV T ensemble is performed next in order to allow a config-urational relaxation of the system (see Figure 1.2) [3]. The displacements of the MC moves areindependent of the potential difference; therefore the simulations are stable, even for such unreal-istic initial configurations. Unbiased molecular dynamics codes would fail, because of the hugeaccelerations.12 CHAPTER 1. PROGRAM MCMDFigure 1.1: Construct a random configuration of a subsystem.Figure 1.2: Monte Carlo simulation to thermalize the subsystem.In a third step, a cubic lattice with a large number of particles is constructed from periodic copiesof the subsystem (see Figure 1.3). The periodic boundary conditions imposed on the subcell avoidconflicts between neighbor cells. At this stage there is no limit for the desired particle number.Additional molecular dynamics simulations are used to randomize and thermalize the large sys-tem (see Figure 1.4). The MD part of the program MCMD is optimized for large particle numbers.The linked cell method is used to reduce the effort for the calculation of pair distances [4]. Thecalculation of the force matrix – which usually comsumes more than 98 percent of the comput-ing time – is parallelized for shared memory and distributed memory parallel computing machines(SGI,SUN,Cray). Parallelization was not a primary goal of this project, but a speedup of 10 is easilyachieved by a simple master–slave approach. The parallelized version of the program is capable ofsimulations with up to one million particles.1.2. METHOD 3Figure 1.3: Build a lattice from the equilibrated subsystem.Figure 1.4: Parallelized molecular dynamics code to thermalize the large system.Chapter 2Usage2.1 CompilationThe program MCMD is written in C. A serial version is running on all standard UNIX platforms. Itwas tested on: IBM AIX, DEC Alpha, SGI Irix, SUN Solaris, Intel based Linux and Cray UNICOS.Optimized and parallelized versions of this program exist for SUN Solaris and SGI Irix sharedmemory architectures and Cray T3E. Shared memory parallelization is done by pragma directives;the version for Cray T3E is based on MPI routines.The compilation is managed by compiler options in the Makefile. They are summarized in Figure2.1. Preprocessor defines are used to build the platform dependent optimized version (SERIAL,MPP, SUNPAR, SGIPAR). Further speedup on SGI Power Challenge machines with Mips R8000processors is obtained by an additional flag: -r8000.2.2 Running JobsThe following lines contain the UNIX commands to start the program on two CPUs:SGI: setenv NUM_THREADS 2; mcmdSUN: setenv PARALLEL 2; mcmdT3E: mpprun -n 2 mcmdDefaultCC = ccCFLAGS = -DSERIALLIBS = -lmCray MPI parallelCFLAGS = -DMPP -h inline0 -h scalar2 -h vector2LIBS = -lm -lmpiSUN shared memory parallelCFLAGS = -DSUNPAR -xO5 -fast -xexplicitpar -xloopinfo -xvparaLIBS = -lmSGI shared memory parallelCFLAGS = -DSGIPAR -O3 -mpLIBS = -lm -mpFigure 2.1: Compiler options42.3. TEST 5control of the program MCMD0 O:MD 1:MC20 n-step blocking of vir., g(r) and energy calc.1500 resolution g(r): number of points15.0 max. g(r) [A]0.4 MC: acceptance ratio4 MD: temp bath 0:strict 1:weak 2: Andersen 3:NVE 4:Nose1 MD: linked cell method; 0=off4 MD: factor between mean density and max. cell density2000.0 MD: max. velocity distribution [m/s]100 MD: n-step waiting before start of acf calculation1 MD: n-step blocking of acf calculation100 MD: resolution acf: number of points1.0 MD: Nose heatbath: (Chi) []1.0e-45 MD: Nose heatbath: (Q) [Js2]1.0 MD: Nose heatbath: (lnS) []Figure 2.2: mh.inp2.3 TestIn order to check the syntax and runtime performance of the programm a test case is implemented.This may be executed by:cp TEST.konfig.xyz konfig.xyzcp TEST.mcmd.inp mcmd.inpcp TEST.mh.inp mh.inpmcmd > outdiff out TEST.outdiff konfig.xyz TEST.sav.konfig2.4 ParametersThe program is controlled by two input files: mcmd.inp and mh.inp. For compatibility reasonsthe first one has the same structure as the input file for the program DMMD, which is under construc-tion at ZAM. The file mcmd.inp contains temperature, volume, number of species in the fluidsystem and their interaction parameters; it also contains simulation control parameters like numberof simulation steps, blocking size, time step and cutoff radius (see Figure 2.3). In mh.inp the typeof the simulation is specified: NV T -MC, NVE-MD or NV T -MD with four different methods forfixing the temperature. Some parameters for the evaluation of system properties are also given inthis file (see Figure 2.2).2.5 Recipe for building a large thermalized configuration1. Monte Carlo run to build and relax a subsystem and to construct a large system from thesubsystem. Necessary steps: Change line 15 in mcmd.inp into false (construct a new configuration). Change line 2 in mh.inp into 1 (use Monte Carlo). Run the program Check the file with the large configuration in konfig.xyz.6 CHAPTER 2. USAGEN=1000 Argon Krypton mixtureparameter V,T3.85408e25 molar volume [A3]300.0 temperature [K]10 MD: time step [fs]10e0 cutoff [A]parameter simulation run500 simulation steps100 blocking0 x0 x0 x.true. restartx xx xparameters for system propertiesx xx xparameters for xyz film filex x0 xParameter Temperature-Scaling0 x0 x0 x400 MD: t_temp_scal [fs]molecule parameters2 n_spec (no. of lennard-jones species)500 500 n_part(i)119.8 167.0 eps(i)3.405 3.633 sig(i)39.948 83.800 mass(i)Ar Kr species_name(i)Figure 2.3: mcmd.inp2. Molecular dynamics for further thermalization of the large system. Change line 15 in mcmd.inp into true (restart from konfig.xyz). Change line 2 in mh.inp into 2 (molecular dynamics). Change line 7 in mh.inp into 0 (strict scaling is about 15 percent faster than Nose´thermostat). Change line 8 in mh.inp into 1 (linked cell method is necessary for good perfor-mance). Run the program on SGI, SUN or Cray machines with up to 8 CPUs.Chapter 3Results3.1 PropertiesSimulations with an Argon–Krypton mixture were performed in order to demonstrate the capa-bilities of MCMD. The potential parameters are given in Table 3.1. The interaction parameters forArgon and Krypton are standard Lennard-Jones values; the Lorentz–Bertholet combination ruleswere used to calculate the Ar–Kr potential parameters. Details of the simulation are listed in Table3.2. Several properties are accessible from simulations with MCMD: internal energy U , pressureP , kinetic temperatures, average displacement per time step, pair distribution functions g(r) forall possible interaction types, velocity autocorrelation functions and velocity distributions for eachspecies [5].Simulated pair distribution functions are shown in Figure 3.1. The first peak of g(r) is shifted withincreasing  to larger distances. The height of the first peak is noticeably influenced by the depthsof the pair potential (compare Table 3.1). The velocity autocorrelation functions for Argon andKrypton are shown in Figure 3.2. The correlation time of Krypton is clearly larger than that ofArgon. The lower mass of Argon leads to – as expected from the Maxwellian velocity distributionfunction – a broader velocity distribution (see Figure 3.3).3.2 BenchmarksThe parameters for the benchmark runs are given in Tables 3.1 and 3.2. The radial cutoff forsimulations with 1,000,000 particles is reduced to 10  10 10 m. The trajectories are calculatedusing the velocity Verlet algorithm [4] with strict temperature scaling is applied in every time step.Computing times per MD time step for systems with one million particles are shown in Table 3.3.These values were obtained from an average interval over 10 time steps. =10 10m =KAr–Ar 3.405 119.8Kr–Kr 3.633 167.0Ar–Kr 3.519 141.4Table 3.1: MD simulation of Ar–Kr mixture. Lennard-Jones pair potential parameters.temperature density number of particles time step cutoff=K =g cm 3 =10 15s =10 10m300 1:605 15; 000  800; 000 10:0 15:0Table 3.2: MD simulation parameters for the Ar–Kr system.78 CHAPTER 3. RESULTS0.0 2.0 4.0 6.0 8.0r /10−10m0.01.02.0g(r)Ar−ArAr−KrKr−KrFigure 3.1: MD simulation of Ar–Kr mixture. Radial pair distribution function g(r) at T = 300 K and = 1; 6 g cm 3.0.0 0.2 0.4 0.6 0.8 1.0t /10−12s0.00.20.40.60.81.0<v(0)v(t)>ArKrFigure 3.2: Simulated velocity autocorrelation functions in a mixture of Argon and Krypton at T = 300 Kand  = 1; 6 g cm 3.3.2. BENCHMARKS 9SUN Enterprise SGI Power Challenge Cray T3Enumber of CPUs 8 8 8time per MD step /s 43,9 45,3 18,5Table 3.3: MD simulation with 1,000,000 particles.PEs MPI code shared memory codetime MD step /s time MD step /s1 14.3 14.32 6.85 7.274 3.90 3.998 2.56 2.3416 2.36 1.5032 3.02 1.0864 4.65 0.88Table 3.4: Simulation with 100,000 particles. Computing times for one MD time step on SGI O2000.Performance measurements on Cray T3E are shown in Figures 3.4 and 3.5. The program exhibitsvery high MFLOPS rates and good scaling on up to eight CPUs. The MFLOPS rate drops no-ticeably with increasing particle number (see values for 8 PEs in Figure 3.4). The master-slaveparallelization of the force matrix calculation leads to more data traffic for larger systems.Speedup data for simulations with 100,000 particles on SGI O2000 and Cray T3E are shown inFigure 3.6. The best speedup of almost 10 is achieved with the shared memory code on SGI O2000using 16 processors. The shared memory code shows better scaling in comparison to the MPI codeon SGI O2000. The scaling of the MPI code is slightly better for Cray T3E than for SGI O2000.Computing times for 100,000 particles are shown Table 3.4.−1500.0 −1000.0 −500.0 0.0 500.0 1000.0 1500.0v /ms−10.0000.0020.0040.0060.008P(v)ArKrFigure 3.3: Velocity distribution in a liquid mixture of Argon and Krypton at T = 300 K and  =1; 6 g cm 3.10 CHAPTER 3. RESULTS0 2 4 6 8 10 12 14 16 18PEs05001000150020002500MFLOPS 15.000 30.000100.000800.000Figure 3.4: MD simulation of Ar–Kr mixture. MFLOPS rates on Cray T3E.0 2 4 6 8 10 12 14 16PEs0.02.04.06.08.010.0Speedup15.00030.000100.000800.000Figure 3.5: MD simulation of Ar–Kr mixture. Speedup on Cray T3E for different particle numbers.3.3. REMARKS 110 4 8 12 16PEs0.02.04.06.08.010.0SpeedupCray T3E, MPISGI O2000, MPISGI O2000, shared memoryFigure 3.6: MD simulation of Ar–Kr mixture with 100,000 particles. Speedup on Cray T3E (MPI code) andSGI Origin2000 (MPI and shared memory code)3.3 RemarksThe method of the program MCMD, to construct a large system from periodic copies of a subsystem,implies the danger of long range correlations. This could be problematic, especially for dipolarsubstances or ions, where electrostatic interactions are important. These correlations can be avoidedby subsequent MC runs – with appropiate boundaries – for each subcell of the large system.The question when equilibrium is reached has been the subject to many discussions. In principleevery accessible property and its fluctuation could be used as a check for thermalization. Internalenergy, pressure and kinetic temperatures are printed out during a run of MCMD. They should bechecked, together with the structural and dynamical functions, which are given at the end of asimulation.Bibliography[1] F. Vesely. Computerexperimente an Flu¨ssigkeitsmodellen. Physik Verlag, Weinheim, 1978.[2] G. Corongiu. Molecular dynamics simulation for liquid water using a polarizable and flexiblepotential. Int. J. Quant. Chem., 42:1209–1235, 1992.[3] J. A. Barker and D. Henderson. What is “liquid”? Understanding the states of matter.. Rev.Mod. Phys., 48:587–671, 1976.[4] D. Frenkel and B. Smit. Understanding molecular simulation. Academic Press, Amsterdam,1996.[5] M. P. Allen and D. J. Tildesley. Computer simulation of liquids. Clarendon Press, Oxford,1987.12

Thermalization of Large Fluid Systems consisting of Lennard-Jones Mixtures

http://juser.fz-juelich.de/record/155122/files/FZJ-2014-04307.pdf

Thermalization of Large Fluid Systems consisting of Lennard-Jones Mixtures

Abstract

Similar works

Full text

Available Versions

Juelich Shared Electronic Resources