Home

DL MESO USER MANUAL

1. ikr2arln 1 we Tij lt TO Tmax U rij ae 10 60 00 Tij gt ro T Tmazx 3 Marko Siggia Worm Like Chain WLC 28 kgT 1 1 y ij ph Tij Tig ij lt Tmax U riz ae a A a ae 10 61 00 Tij gt Tmar 4 Morse potential bond 35 U rij De 1 exp 8 ri ro 10 62 where Tij Ti T ro is an equilibrium bond length rmaz the maximum specified bond length or extension k is a spring force constant A the persistence length of a wormlike chain De the potential well depth and the potential width The coefficients for these bonds as used in DL_MESO_DPD and stated in the FIELD file are summarized in Table 10 4 10 6 BOND INTERACTIONS BETWEEN PARTICLES 135 Table 10 4 Bond coefficients used in DL_MESO_DPD Bond type a b c d 1 Harmonic K ro 2 FENE K ro Liat 3 WLC Ap Tmax 4 Morse De To B The force on particle due to a bond potential is obtained from the general formula 5 1 o al il ET Orij Tij with the force acting on particle j equal to the negative of this and the virial contribution from the stretching bond given by W Fij Fi 10 64 with only one contribution per bondt 10 6 2 Bond angles DL_MESO_DPD includes two methods for modelling potentials and forces between three bonded particles due to the angle formed between them 0 x The potentials are given as follows 1 Harmonic K U bijk 3 Oije 00 10 65
2. e Function Outputs concentration of solute iprop and compressible fluid velocity at each lattice point for system in Structured Grid Legacy VTK format e Dependencies fOutputLegacyVTIKCB2D fOutputLegacyVTKCB3D e Arguments filename input array of characters e Comments The default output file name is 1bout vtk and the concentration of solute 0 is output by default These can be changed by specifying an output file name and solute number up to 1bsy ns 1 when calling the routine fOutputLegacy VTKCBIncom e Header records int fOutputLegacyVTKCBIncom const char filename lbout int iprop 0 e Function Outputs concentration of solute iprop and incompressible fluid velocity at each lattice point for system in Structured Grid Legacy VTK format e Dependencies fOutputLegacyVTKCB2DIncom fOutputLegacyVTKCB3DIncom e Arguments filename input array of characters e Comments The default output file name is lbout q and the concentration of solute 0 is output by default These can be changed by specifying an output file name and solute number up to 1bsy ns 1 when calling the routine fOutputLegacyVTKT e Header records int fOutputLegacyVTKT const char filename 1bout e Function Outputs macroscopic temperature and compressible fluid velocity at each lattice point for system in Structured Grid Legacy VTK format e Dependencies fOutputLegacyVTKT2D fOutputLegacyVTKT3D 6 2 DLMESO_LBE SUBROUTINES AND
3. e manybody_module Contains subroutines for calculating density dependent forces between particles including the calculation of localized densities e bond module Contains book keeping and force calculation routines for bonds angles and dihedrals e domain_module Contains subroutines to construct parallel link cells import and deport particles in and export particle data to domain boundary halos 139 140 CHAPTER 11 DL_MESO_DPD PACKAGE REFERENCE e start_module Contains subroutines to initialize and restart DPD calculations e config module Contains subroutines to read in input files with system and molecule bond data to zero all parameters and accumulators for statistical data and for parallel version only determine 3D domain decomposition e field module Contains subroutines to calculate pairwise forces between particles e integrate_module Contains subroutines to integrate the equations of motion using the Velocity Verlet scheme and apply various thermostats and barostats e statistics module Contains subroutines to calculate and write out statistical and trajectory data e run_module Contains program loops for different integrators and barostats Most of the above modules and the main program file dlmesodpd 90 are identical for both serial and parallel versions of DL_MESO_DPD The filenames for the serial versions of comms_module and domain_module end with _ser but are referred to in the code by their
4. 1 1 1 1 2 2 1 Cia 1 1 1 1 1 1 3 10 12 4 9 13 18 Speed vector 12 13 14 15 16 17 18 4 1 LATTICE MODELS D3Q19 11 11 4 1 1 1 11 4 1 1 2 ooo oo ooo Oo 1 30 12 0 0 0 0 0 0 0 30 r LTs jne TF 9 bulk CHAPTER 4 DL_MESO_LBE BASIC DEFINITION p p 19 e amp 4 llp po is tiy J eea wept dy Pda z de qa 3Jx dy jy a 3ty iz de qsa ie 3Jz enn La 2j j2 32 TT PO ca 2 2 p am po 252 jy Ja Psw tute med Wax 2 2 apap po is 92 Pry daly 7 Po Pye Jylz ed Po pe dada mo n ma 0 m 0 0 38 07 Fo Vy Fy v2F 11 vs Fs vyFy 004 E Fx 2F Fy 3Fy Fz S 2 2u Fy vyFy uz Fz 207Fy vy Fy vz F 2 2u fy UzFs vyFy v F Uy Fy vy Fo Uy EF uzFy v2 Fy UL Fs 0 0 0 T 1 1 1 1 1 1 bulk 52 1 84 1 84 1 84 T 154 TF 94 T Tp Tp 16 516 516 a aay 31 4 1 LATTICE MODELS D3Q27 Weight factor 4 al 3 PISAS Ko O NIN IN I ol jo SHI IN elOl mim al M TFT a ha re Speed vector Ci z Ci y Cia 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 32 CHAPTER 4 D
5. DL_MESO_LBE can write standard Plot3D Structured Grid Legacy VTK and Structured Grid XML VTK data format files in binary format for parallel calculations and ANSI for serial These may be modified by the user as required The utility lbeplot3dgather in the LBE utility directory can combine Plot3D files generated in parallel while Parallel Structured Grid XML VTK files 1btout pvts that refer to the files from each processor can be created using the utility lbevtkgather further details can be found in Appendix B or the README file in the same directory Plot3D format Output grid position lbout xyz NL NY NZ T0 0 0 Tngz 1 ny 1 nz 1 Y0 0 0 Ynz 1 ny 1 nz 1 Z0 0 0 Zng 1 ny 1 nz 1 where nx is the total number of grid points in x direction ny is the total number of grid points in y direction nz is the total number of grid points in z direction and 2 j k Yi j k Zi j k is the Cartesian coordinate of grid point i j k Output macroscopic quantities 1bout q NX NY NZ c 1 0 Re t P0 0 0 Pna 1 ny 1 nz 1 Uo 0 0 Uiicinyonsai Vo 0 0 tee Var thy imesti Wo 0 0 Wrz 1 ny 1 nz 1 D0 0 0 gt SB Ona 1my 1 n2 1 where c is the speed of sound for the lattice Re the Reynolds number for the flow t the time step p the density U the x component of velocity V the y component of velocity W the z component of velocity and the space boundary condition property The 1
6. The boundary condition number can be rather confusing and difficult to understand The GUI in DL_MESO therefore includes a translator which interprets a word as its corresponding integer number The word is made 4 2 DATA STRUCTURE 33 Table 4 2 Boundary condition category letter meaning V Constant Velocity P Constant Pressure Density C Constant Solute Composition T Constant Temperature B Neumann Boundary Condition Solute Composition or Temperature PS Planar Surface CC Concave Corner CE Concave Edge T Normal Vector Pointing to Top D Normal Vector Pointing Downwards L Normal Vector Pointing to Left R Normal Vector Pointing to Right F Normal Vector Pointing to Front B Normal Vector Pointing to Back up of defined letters as listed in Table 4 2 The boundary conditions with combinations of type and orientation are listed in Table 4 3 The letters are in the order of 1 2 3 Fluid property constant speed or constant pressure Solute property constant composition or Neumann boundary Temperature property isothermal constant or heat bath Neumann boundary Geometric property planar surface concave corner or concave edge Boundary orientation one letter for planar surface two letters for concave corners or three letters for concave edges For example a shearing planar surface facing down the y axis with constant composition and temperature i e isothermal is represent
7. The boundary width for running the parallel version of DL_MESO_LBE is set using this box The serial version by default automatically resets the boundary width to zero The sound speed c is the real life quantity for the first main fluid The kinetic viscosity v is the real life quantity for the first main fluid and together with the sound speed allows conversions between lattice and real units the time step and lattice spacing are given by oe V e V3v At AEE and Ag E The specified top boundary speed in lattice units i e at the maximum y value for the grid This respectively and similar properties only need to be specified if the boundary includes a fixed velocity any specified value will be ignored for periodic bounce back and fixed density conditions Note that the z component will be greyed out for two dimensional systems The specified bottom boundary speed in lattice units i e at the minimum y value for the grid The specified left boundary speed in lattice units i e at the minimum x value for the grid The specified right boundary speed in lattice units ie at the maximum zx value for the grid The specified front boundary speed in lattice units i e at the maximum z value for the grid For two dimensional systems this can be omitted and is thus greyed out The specified back boundary speed in lattice units i e at the minimum z value for the grid For two dimensional systems this can
8. critigroup output integer reference e Comments The totalgrid grid points are distributed among totalgroup processes so that the first critigroup processes have indigrid grid points and the others have indigrid 1 fCppMod e Header records two cases int fCppMod int a int b long fCppMod long a long b e Function Ensure that a is in a range between 0 and b 1 so that the value beyond the maximum value equals the minimum and vice versa e Dependencies None 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 55 e Arguments a input integer long integer b input integer long integer fCppMod output integer long integer e Comments fCppMod O when a b or fCppMod b 1 when a 1 This function is useful for periodic boundary conditions fPrintLine e Header records int fPrintLine e Function Prints a line of 72 characters e Dependencies None fPrintDoubleLine e Header records int fPrintDoubleLine e Function Prints a line of 72 characters e Dependencies None fRandom e Header records double fRandom e Function Creates a double precision random number between 1 and 1 e Dependencies None e Arguments fRandom output double precision e Comments There is a built in seed in the function which is only activated when the function is initially called fBigEndian e Header records int fBigEndian e Function Determines endianness for machine returns 1 for big end
9. e Header records SUBROUTINE exportdensity nlimit mdir mp begin final shove Parallel version SUBROUTINE exportdensity nlimit mdir begin final Serial version 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 161 e Function Exports particle data local densities to neighbouring domain as boundary halo for calculation of many body DPD interaction forces e Dependencies error global_sca_and Parallel version only msg_receive_unblocked Parallel version only msg_receive_sca_blocked Parallel version only msg_send_blocked Parallel version only msg_send_sca_blocked Parallel version only msg wait Parallel version only e Arguments nlimit input output integer mdir input integer mp input integer begin input real KIND dp final input real KIND dp shove input real KIND dp deportdata e Header records SUBROUTINE deportdata nlimit e Function Applies deport of particles from boundary halos to neighbouring domains and or periodic boundary conditions e Dependencies deport Parallel version only e Arguments nlimit input output integer importdata e Header records SUBROUTINE importdata nlimit e Function Applies import of particle forces from boundary halos of neighbouring domains and or across periodic boundary conditions e Dependencies import Parallel version only e Arguments nlimit input output integer e Comments This subroutine is applicable for integrators thermostats
10. system temperature at current step stpttp real KIND dp system bond potential energy at current step stpbe real KIND dp system angle potential energy at current step stpae real KIND dp system dihedral energy at current step stpde real KIND dp system electrostatic energy at current step stpee real KIND dp system mean bond length at current step stpbdl real KIND dp system mean bond angle at current step stpang real KIND dp system mean bond dihedral at current step stpdhd real KIND dp rolling average system potential energy ravpe real KIND dp rolling average system virial ravvir real KIND dp rolling average system kinetic energy ravtke real KIND dp rolling average system total energy ravte real KIND dp rolling average system pressure ravprs real KIND dp rolling average system volume ravvlm real KIND dp rolling average system temperature ravttp real KIND dp rolling average system bond potential energy ravbe real KIND dp rolling average system angle potential energy ravae real KIND dp rolling average system dihedral potential energy ravde real KIND dp rolling average system electrostatic potential energy ravee real KIND dp stage number for data stack nav integer data stack for potential energy stkpe real KIND dp A mxstak data stack for virial stkvir real KIND dp A mxstak data stack for kinetic energy stktke real KIND dp A mxstak data stack for volume stkvlm real KIND dp A mxstak data stack for bond potential ene
11. 1 when calling the routine fOutputVTKT e Header records int fOutputVTKT const char filename lbout e Function Outputs macroscopic temperature and compressible fluid velocity at each lattice point for system in Structured Grid XML VTK format e Dependencies fOutputVTKT2D fOutputVTKT3D 76 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE e Arguments filename input array of characters e Comments The default output file name is lbout vts This can be changed by specifying an output file name when calling the routine fOutputVTKTIncom e Header records int fOutputVTKTIncom const char filename lbout e Function Outputs macroscopic temperature and incompressible fluid velocity at each lattice point for system in Structured Grid XML VTK format e Dependencies fOutputVTKT2DIncom fOutputVTKT3DIncom e Arguments filename input array of characters e Comments The default output file name is lbout vts This can be changed by specifying an output file name when calling the routine Other output routines Subroutines with names fsQutput are suitable for serial running and produce output files for entire systems Unlike the routines listed above these omit any domain boundary lattice points used in calculations Notes regarding q vtk and vts files e lboutxz is the q vtk or vts file at the xth saved step during serial running e lboutyatx is the q vtk or vts file at the xth saved step written
12. 2 Harmonic cosine U Oij 5 cos Ox cos 6o 10 66 where is an angle force constant and o an equilibrium bond angle Table 10 5 gives the assignments used by DL_MESO_DPD and the FIELD file for the bond angle coefficients with angles given in degrees and cosines automatically calculated from angles Table 10 5 Bond angle coefficients used in DL_MESO_DPD Angle type a ble d 1 Harmonic kK A 2 Harmonic cosine k 09 The angle across particles j and k can be determined from the bond vectors rj T Fj and Tj Tk T Oije cos uta 10 67 Tijfkj The most general form for the bond angle potential is given thus U 8 56 Tij Thi A Oije S riz S rki S rik 10 68 4This expression is also used for the virial contribution from the standard DPD pairwise forces i e Equation 8 3 again only applying a single contribution per particle pair 136 CHAPTER 10 DL MESO_DPD FEATURES with A 0 as a purely angular function and S r a screening or truncation function The force on particle in dimension a is thus given by 0 ES ago Oise Mian Mes ai o riz S rki S rik Bra 051 rij A Oije S rri S vin ei 005 28 ris Tij Tij Toj fa A ijn S riz S rik Sex e ree rk I J 8 A ijr S rig S 145 Bee Ses ES rin 10 70 ik Oik with 6g 1 if a b and dq 0 if a b
13. 8 1 Introduction ecc a a ee aa 8 2 Outline of Method sari a atp dat 8 3 Equation of state and dynamic properties 8 4 Derivation of Equilibrium a 8 5 Summary of Dissipative Particle Dynamics DL_MESO_DPD Basic Definition gTr Data structure ed oh ays er Sek Be E aR ie i 9 2 The Parameters and Their Functions DL_MESO_DPD Features 10 1 Domain decomposition and linked list cell calculations 10 2 Thermostats and integration algorithms 10 3 Barostats eee cuse a ee ee ee bee ca Be ee ees 10 4 Particle particle interactions 2 e a 02000000004 10 5 Long ranged Electrostatic Coulombic Potentials 10 6 Bond interactions between particles 000 10 7 Surface interactions ssi t studa 6 ae BP ban Ot pe ats DL_MESO_DPD Package Reference LD le OVerviGw a A Oe RE ein ee DA Pee eg 11 2 DL_MESO_DPD Subroutines and Functions 11 3 DL MESO_DPD Data Files DL_MESO DPD Examples 12 1 Mixture Sereda de dd De A A Ed alate 12 2 Ph seSeparabiOn eta a A A at 12 3 APTA a ae eee eee De ean ae A o da 12 4 Polyelectrolyte si a Seg a asp eg Sm pod Nae wee ES 12 5 AmphiphileMesophases 0 000 00000 eee eee 12 6 Mixture Par a gerei poids done did od ok AOE aed At Ba ae Mee a 12 7 VesicleFormation as elt 6 id Pana is a ee E Manual compiliation and running of DL_MESO Ay DE MES
14. If running DL MESO_DPD in serial the modules comms_module and domain_module should be replaced by comms_module_ser and domain_module_ser respectively To simplify the process a makefile may be created either in the DPD directory or in the working directory to automatically compile the modules and build the executable Examples of these for running in the DPD directory may be found in the DPD makefiles directory and modified by the user the compiler after FC and flags after FFLAGS may need changing depending on the Fortran90 compiler available If invoking from the working directory the modules for DL_MESO_DPD should either be preceded by the path i e DPD in the list of compile sources or the directive VPATH DPD can be used before the source list the latter strategy is used by the DL_MESO GUI when creating makefiles DL_MESO_DPD can be compiled by the command make if the makefile is called Makefile or if it has a custom name e g Makefile custom by the command A 2 DL MESO DPD 185 e make f Makefile custom The example makefiles will produce an executable with the name dpd exe which can be copied to the working directory if necessary and run at the command line after the required input files CONTROL and where appli cable MOLECULE and FIELD are created or copied into the same directory A CONFIG file can also be placed in the same directory if an initial configuration is required for a new simulation export fil
15. Po Tp with the option of keeping the system isotropic by using equal values of P for all dimensions This may be switched off if imposing constant surface tensions at given planar surfaces Rescaling of the simulation volume takes place in the first stage while calculation of the rescaling factor ug takes place in the second stage after force integration 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 171 11 2 15 statistics_module This module requires the modules constants variables start_module comms_module and numeric_container to be loaded beforehand printout e Header records SUBROUTINE printout time lbegin e Function Writes summary of simulation at the current time step to OUTPUT file e Dependencies None e Arguments time input real KIND dp lbegin input logical e Comments The logical 1begin indicates that the column titles are to be printed before the data corout e Header records SUBROUTINE corout time e Function Writes statistical data to CORREL data file e Dependencies None e Arguments time input real KIND dp histout e Header records SUBROUTINE histout time e Function Writes trajectory data including virial to HISTORY data files e Dependencies None e Arguments time input real KIND dp 172 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE result e Header records SUBROUTINE result e Function Writes final summary of simulation to OUTPUT file e Dependencies rev
16. e Header records two cases 3D and 2D int fGetCoord long tpos int amp xpos int amp ypos int amp zpos int fGetCoord long tpos int amp xpos int amp ypos e Function Calculates the Cartesian coordinate of a grid point from its position in a one dimensional array e Dependencies None e Arguments xpos output ypos output zpos output tpos input fGetOneMassSite integer reference integer reference integer reference long integer e Header records three cases double fGetOneMassSite double startpos double fGetOneMassSite int fpos long tpos double fGetOneMassSite int fpos int xpos int ypos int zpos e Function Calculates the mass density of one of the fluids at a grid point e Dependencies None e Arguments startpos fpos tpos xpos ypos Zpos input double pointer input integer input long integer input integer input integer input integer fGetOneMassSite output double precision e Comments Mass density is calculated according to the definition p gt fi In the first case startpos is the start point for the summation of particle distribution functions and must be assigned correctly The second case carries out the same calculation for the fpos th fluid and tpos th grid point while the third carries it out for the fpos th fluid at the grid point indicated by the Cartesian coordinate xpos ypos zpos The latter two are more readable to the user but a bit slower to c
17. e The system clock is checked timchk e The first stage of the required integrator is used to update the motion of the particles and their positions e Masses and particle molecule names are reassigned to the particles e The parallel link cell structure is set up and the pairwise forces calculated plcfor_ e The time taken to calculate the forces is determined timchk e The second stage of the integrator is applied to calculate the velocities at the end of the time step This may include e g recalculation of dissipative forces or resizing of the system for a barostat e Statistical properties for the system are calculated and during equilibration the particle velocities are rescaled for the specified temperature statis e After every nsbpo time steps the system clock is consulted timchk and statistical data is printed to the OUTPUT file printout If equilibration has come to an end i e nstep nseql this is also reported equilout e Ifrequested by the user after every iscorr time steps statistical data is written to the CORREL file corout and after every nsbsv time steps trajectory data is saved to the HISTORY or HISTORY file s e The step time is calculated timchk and if the allocated time has expired the job is closed down otherwise restart data is saved revive After all the time steps have been calculated or the allocated time has elapsed e The final calculation summary is printed to the OUTPUT file
18. pressible fluids e Dependencies Many for different lattice schemes e Arguments tpos input long integer prop input long integer uwall input array of doubles e Comments Planar surface calculations are based on 21 concave edges and corners assume zero speed at the boundary The fluid velocity at the lattice point given by the array uwal1 is required for this boundary condition 82 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE fFixedTemperature e Header records int fFixedTemperature long tpos int prop double uwall e Function Calculates the particle distribution function at a fixed temperature boundary for compressible and incom pressible fluids e Dependencies Many for different lattice schemes e Arguments tpos input long integer prop input long integer uwall input array of doubles e Comments Planar surface calculations are based on 21 concave edges and corners assume zero speed at the boundary The fluid velocity at the lattice point given by the array uwall is required for this boundary condition fPostCollBoundary e Header records int fPostCollBoundary e Function Calculates the particle distribution function at different boundaries after the collision step prior to prop agation for compressible fluids e Dependencies Many for different lattice schemes e Comments Algorithms for other boundary conditions can be added by the user although care should be taken as to when they are
19. 1 0 1 0 1 0 1 0 1 0 p P ee 20 3 j jy es wep M2 432 M9 go jz eq as ace Poh om vis jody 0 6 vz Fy VyFy 6 0 Fs ty Fy Y II a F 2 Ur Fy Uy Fy Un Py vy Fy T 1 mh ea a Gia 1 sa 1 54 T TF 1 1 Ax v ron 5 Ari 2 At 27 4 1 LATTICE MODELS D3Q15 Wi aoo 72 Weight factor 1 3 8 10 4 7 11 14 Speed vector Ci z Ciy Cix 10 11 12 13 14 1 4 1 4 1 4 4 4 4 16 1 1 1 1 1 1 4 1 1 1 1 1 E A 0 heel Ht 0 28 1 7 CHAPTER 4 DL_MESO_LBE BASIC DEFINITION p Cn eea 2p Gett H eo wep FH ja jy j2 qe de dy Jy Gy y Jz Jz qs rie 3ne4 2j2 j2 j pes jody pet dude pet dede ae J 0 2 Ug Ey Vy Fy vzF 2 10 vz Fr vy Fy vF Fo iF Fy 2 3Fy S F FF 2 2u Fr Uy Fy v Fz 2 vy Fy vzF Uy Fy vy Fy Vy uzFy vz Fy Vrt 0 T 1 1 1 1 1 1 bulk 5251 84 1 84 1 84 T Tp TF Tf Tf 514 V 1 2 Ar At 2 Tf bulk 9 29 1 I 2 2 1 2 2 1 11 4 4 1 1 11 4 4 1 2 11 4 4 4 8 1 36 Ci z hiena Wi 1 1 2 2 1 1 1 Weight factor Ciy
20. 1 The type of thermostat to be used in the simulation 2 Thermostat parameters for the Stoyanov Groot thermostat currently the only type that requires an additional parameter a global coupling parameter for the Nos Hoover part a is required Figure 2 6 d gives the layout for the barostat pop up window 1 The type of barostat to be used in the simulation 2 Barostat parameters for the Langevin barostat a barostat relaxation time 7 and piston drag coefficient yp are required while the Berendsen barostat requires the compressibility relaxation ratio 2 3 This check box determines whether or not an isotropic system i e one where pressure acts uniformly in all dimensions should be modelled If unchecked the barostat will act differently in each dimension and the shape of the system will change over time The parameters for electrostatics can be given in the pop up window shown in Figure 2 6 e 1 The type of electrostatics to be used in the simulation 2 Electrostatic parameters for the Ewald sum with Slater type smearing the system coupling constant T Ewald real space convergence a and charge smearing 8 coefficients need to be specified 3 The Ewald sum method also requires a reciprocal space k vector range If non periodic boundaries are to be used the parameters for surfaces can be entered in the appropriate pop up window Figure 2 6 f 16 CHAPTER 2 THE DL_MESO GUI DL_MESO LBE
21. In the absence of screening terms the above formula reduces to 0 FA A 0 10 71 are Gist 10 71 1 S a A sin ijk olijk ik fe Tj ro ro Tkj 02 02 e Sek cos 0 51 02 02 3 Sp vk 5 10 72 Tijf kj Tijfkj Tij Tkj The contribution to the virial from the angle is given by W a A F Tkj Fa 10 73 which is equal to zero for bond angle potentials without screening terms 47 10 6 3 Bond dihedrals Three potential models for bond dihedrals along particles 7 j k and l are provided in DL_MESO_DPD as follows 1 Cosine torsion U di x1 A 1 cos moizr 9 10 74 2 Harmonic K U Pijki 5 Dist do 10 75 3 Harmonic cosine K U ijt 5 cos ijki cos bo 10 76 where A and are dihedral force constants m is the multiplicity 6 the dihedral at minimum potential and go an equilibrium bond dihedral The assignment of coefficients used in DL_MESO_DPD and the FIELD file is given by Table 10 6 Table 10 6 Bond dihedral coefficients used in DL_ MESO_DPD Dihedral type alblcld 1 Cosine torsion A m 2 Harmonic kK do 3 Harmonic cosine k do The dihedral angle across all four particles or between planes ij and kl is given by dig cos B Fij Fito Tkl 10 77 10 7 SURFACE INTERACTIONS 137 where a 3 E Fij X Tjk Tir X Te Fij X Polym X Fil B Tij Tik Tki 10 78 whi
22. Reads in molecule and bond data from MOLECULE and FIELD files Dependencies MOLECULE file FIELD file getdble getint getword error global_sca_and Header records SUBROUTINE zero Function Sets step counters initial time parameters system parameters accumulators for statistical properties and long range potential corrections to zero initializes random number generators 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 167 e Dependencies duni mtrnd forgen e Header records SUBROUTINE forgen e Function Sets up arrays for potential and force parameters e Dependencies error global_sca_and elecgen e Header records SUBROUTINE elecgen e Function Sets up electrostatic parameters for self interaction and charged system corrections e Dependencies error global_sca_and 11 2 13 field module This module requires the modules constants variables bond_module manybody_module ewald_module surface module numeric_container and domain_module to be loaded beforehand Different versions of each subroutine are available in this module for different integrators and or thermostats DPD thermostat with standard molecular dynamics Velocity Verlet integration MD VV 54 mdvv DPD thermostat with DPD Velocity Verlet integration DPD VV 3 dpdvv e Lowe Andersen thermostat 27 lowe Peters thermostat 38 peters Stoyanov Groot thermostat 49 stoyanov forces_ e Header records SUBROUTINE
23. gt qi qchg real KIND dp Ewald real space convergence parameter a alphaew real KIND dp 27k gT At 2This non global value is determined according to the number of particles and link cells required not required for DPD thermostat MD VV or DPD VV 3Only required for Stoyanov Groot thermostat This incorporates the time step for Velocity Verlet integration and is thus equal to 9 2 THE PARAMETERS AND THEIR FUNCTIONS 121 Table 9 4 DL MESO_DPD Parameters continued function parameter data type notation reciprocal of real space convergence parameter ralphaew real KIND dp Slater charge smearing coefficient 8 betaew real KIND dp maximum reciprocal vector size in x dimension KIIF kmax1 integer gt 14 maximum reciprocal vector size in y dimension ky kmax2 integer maximum reciprocal vector size in z dimension kg kmax3 integer Ewald self interaction correction engsic real KIND dp charged system correction qfixv real KIND dp bond interaction parameter a aabond real KIND dp A mxbonddef bond interaction parameter b bbbond real KIND dp A mxbonddef bond interaction parameter c ccbond real KIND dp A mxbonddef bond interaction parameter d ddbond real KIND dp A mxbonddef angle interaction parameter a aaang real KIND dp A mxbonddef angle interaction parameter b bbang real KIND dp A mxbonddef angle interaction parameter c ccang real KIND dp A mxbonddef angle interacti
24. int main int argc char argv bigend fBigEndian fStartMPI argc argv fDefineSystem lbin sys fDefineDomain fStartDLMESO fMemoryAllocation fDefineNeighbours O fDefineMessage fPrintSystemInfo fAllReady fPrintDomainInfo fAllReady fInputParameters lbin sys fReadSpaceParameter lbin spa fGetModel fBoundNonBlockCommunication fMarkBoundArea fInitializeSystem fReadInitialState lbin init fAllReady fOutputInfo f OutputGrid Only needed for Plot3D output files timetotal fCheckTimeMPI for int i 0 i lt lbtotstep i fNonBlockCommunication if i lbsave 0 fOutputVTK Substitute with required file type and or output qVersion if 1bdm rank 0 printf d Xt i fPrintSystemMass fPrintSystemMomentum fInteractionForceShanChen fForceNonBlockCommunication CollisionBGK fPostCollBoundary fPropagationCombSwap fPostPropBoundary timetotal fCheckTimeMPI fFreeMemory fFinishDLMESO fCloseMPI return 0 51 52 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE 6 2 2 IbpMODEL D2Q9 e Header records int D2Q9 e Function Assign the weight factor lbw speed vector 1bv index for opposite speed vector lbopv MRT transformation matrices 1btr and 1btrinv and tuneable parameters 1lbmrts and 1bmrtw for the D2Q9 lattice model e Dependencies None D3Q15 e Hea
25. nents This method is based on calculating pairwise interaction potentials by using an effective mass for each component 4 which is a function of density and most frequently defined as ya Z 1 exp p 5 14 where p is the density of component a This can be changed by the user by modifying the subroutine fCalcPotential_ShanChen Defining gap as the interaction coefficient between components a and b the overall force on component a due to interactions with other components is defined as ES 4 8 X gab X with E i i 5 15 b i 5 4 DIFFUSION AND HEAT TRANSFER 47 and any suitable forcing algorithm can be used to apply this force on a Lattice Boltzmann Equation simulation For a particular interaction the resulting equation of state 41 is defined as 1 P pes 5 5Gav p 5 16 and the interfacial tension between the phases 42 is given as 4 2 oo 2 eagapCs Ax i dy ab mnnm 1 Oab 2 Fags de dz 5 7 where e is a lattice dependent constant and z distance along the normal from the phase interface 5 4 Diffusion and heat transfer In a similar fashion to multiple fluid systems the Lattice Boltzmann Equation method can be applied to problems involving diffusion and or heat transfers by using additional distribution functions for each solute and or temperature 20 59 For a system consisting of a number of dilute solutes along with a bulk fluid the governing equation for each solu
26. result e The duration of the calculation run is determined and printed timchk e All remaining output channels for OUTPUT HISTORY files are closed e For parallel running only MPI communications are closed down exitcomms 11 2 2 numeric_container This package contains general purpose functions which may be replaced with any suitable functions in Fortran90 standard libraries as well as bookkeeping subroutines for the global local particle number list duni e Header records REAL KIND dp FUNCTION duni idnode e Function Creates a double precision random number between 0 and 1 e Dependencies None e Arguments idnode input integer duni output real KIND dp e Comments The random number generator is an implementation of the Universal Random Number Generator 30 The processor name idnode is used as a seed which is activated the first time the function is called 142 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE mtrnd e Header records REAL KIND dp FUNCTION mtrnd idnode e Function Creates a double precision random number between 0 and 1 e Dependencies None e Arguments idnode input integer mtrnd output real KIND dp e Comments The random number generator is an implementation of the Mersenne Twister random number generator 33 The processor name idnode is used as a seed which is activated the first time the function is called gaussmp e Header records REAL KIND dp FUNCTION gaussmp idn
27. z component of interparticle distance dzz real KIND dp A mxpair number of particle pairs in thermostat pair list npair integer system wide virial corrections vrlsg real KIND dp A nsyst barostat type btype integer barostat target pressure Po prszero real KIND dp thermostat parameter a atherm real KIND dp thermostat parameter b btherm real KIND dp thermostat parameter c ctherm real KIND dp thermostat parameter d dtherm real KIND dp x component of piston velocity upx real KIND dp y component of piston velocity upy real KIND dp z component of piston velocity upz real KIND dp x component of piston velocity at previous timestep upx1 real KIND dp y component of piston velocity at previous timestep upy1 real KIND dp z component of piston velocity at previous timestep upzi real KIND dp piston mass Wg psmass real KIND dp x component of piston force fpx real KIND dp y component of piston force fpy real KIND dp z component of piston force fpz real KIND dp instantaneous virial ivrl real KIND dp A 3 Langevin random force parameter op sigmalang real KIND dp electrostatic algorithm type etype integer electrostatic interaction parameter a aelec real KIND dp electrostatic interaction parameter b belec real KIND dp electrostatic interaction parameter c celec real KIND dp electrostatic interaction parameter d delec real KIND dp electrostatic coupling parameter T gammaelec real KIND dp total system charge gt
28. 0 between the lattice speed of sound and flow Reynolds number represents the freestream angle of attack The density of the fluid may be replaced with its concentration or scalar temperature 6 3 DL MESO_LBE DATA FILES 103 Legacy VTK format lbout vtk Legacy VTK files written by DL_MESO_LBE include the lattice dimensions numbers of grid points in each direction the Cartesian coordinates of the grid points in real life units and the following data e A scalar property fluid density solute concentration or scalar temperature e g P0 0 0 Pna 1 ny 1 nz 1 e Fluid velocity Uo 0 0 Vo o 0 Wo 0 0 Una 1 ny 1 nz 1 Vne l ny 1l nz 1 Wha 1 ny 1 nz 1 e The space boundary condition property 0 0 0 Qnesinye igi Structured Grid XML VTK format lbout vts Structured Grid VTK files written by DL_MESO LBE include the lattice dimensions numbers of grid points in each direction the Cartesian coordinates of the grid points in real life units and the same data as in Legacy VTK files Output files produced in serial include the data between XML tags e g lt DataArray gt while those produced in parallel use the lt DataArray gt tags to refer to the starting point for the data in a stream of binary numbers inside an lt AppendedData gt tag The latter represent the data in each subdomain and should be retained when plotting the entire system since the parallel VTK format only links to the piece VTK files Cha
29. 15 Hans C Andersen Molecular dynamics simulations at constant pressure and or temperature Journal of Chemical Physics 72 4 2384 2393 1980 H J C Berendsen J P M Postma W F van Gunsteren A DiNola and J R Haak Molecular dynamics with coupling to an external bath Journal of Chemical Physics 81 8 3684 3690 1984 Gerhard Besold Ilpo Vattulainen Mikko Karttunen and James M Polson Towards better integrators for dissipative particle dynamics simulations Physical Review E 62 6 R7611 R7614 December 2000 P L Bhatnagar E R Gross and M Krook A model for collision processes in gases I Small amplitude processes in charged and neutral one component systems Physical Review 94 3 511 525 May 1954 Shiyi Chen and Gary D Doolen Lattice Boltzmann method for fluid flows Annual Review of Fluid Dynamics 30 1 329 364 1998 Paul J Dellar Bulk and shear viscosities in lattice Boltzmann equations Physical Review E 64 3 031203 August 2001 Dominique d Humieres Irina Ginzburg Manfred Krafczyk Pierre Lallemand and Li Shi Luo Multiple relaxation time lattice Boltzmann models in three dimensions Philosophical Transactions of the Royal Society of London Series A Mathematical Physical and Engineering Sciences 360 1792 437 451 2002 Burkhard Diinweg and Wolfgang Paul Brownian dynamics simulations without Gaussian random numbers International Journal of Modern Physics C 2 3 817 827 1991 U Fris
30. 2 o ti f r exp 28 3 4 vib where 9 kgT m kg is the Boltzmann constant T is temperature m is molar mass D is the space dimension is the thermal velocity and w the macroscopic velocity When l lt v Equation 3 4 can be expanded into pola pow e En Ea 270 20 B o 26 2 3 5 For a microscopic quantity Y the associated macroscopic quantity Y is calculated by w ue suas 3 6 Let E v20c where is a rescaled thermal velocity the macroscopic velocity can be similarly rescaled to v20u Equations 3 5 and 3 6 can thus be combined to give v feyt A 1 2 2 0 2 2 8 u de 3 7 Using Gaussian quadrature Equation 3 7 changes into V20p be ae ae T eae 142 7 2 2 4 u 3 8 Let a 3 9 270 and fe wip i 2 42a ul 3 10 a 3 4 STRUCTURAL RELAXATION AND MACROSCOPIC EQUATIONS 23 The value of w c can be obtained from Gauss Hermite integration Equation 3 10 is the equilibrium particle distribution function in the discrete regime w is called the weight factor for speed vector c Equation 3 10 can also be written as 3 t 7 9 u 3u de c2 2c4 2c fi wip eh where c V30 2T is the modulus of the basic lattice vector For water at 20 C c 367 8 m s m 3 4 Structural Relaxation and Macroscopic Equations The Lattice Boltzmann method often uses the BGK Bhatnager Gross and
31. DL_MESO_LBE Parameters continued 39 function parameter data type notation position for sending boundary message position for receiving boundary message position for sending force message position for receiving force message message type message type message type message type message type message type message type message type message type message type message type message type message type message type message type endianness of system total calculation time output file number lbnb bspos lbnb brpos lbnb fspos lbnb frpos 1bmsg2x 1bmsg2y 1bmsg3x 1bmsg3y 1bmsg3z 1bbmsg2x 1bbmsg2y 1bbmsg3x 1bbmsg3y 1bbmsg3z lbfmsg2x lbfmsg2y 1b msg3x lbfmsg3y lbfmsg3z bigend totaltime qVersion unsigned long unsigned long unsigned long unsigned long MPI_Datatype MPI_Datatype MPI Datatype MPI Datatype MPL Datatype MPL Datatype MPL Datatype MPI Datatype MPI Datatype MPL Datatype MPL Datatype MPL Datatype MPI Datatype MPI Datatype MPI_Datatype int double int Chapter 5 DL MESO LBE Features 5 1 Collision and Propagation Algorithms The collision and propagation routines are a fundamental part of Lattice Boltzmann Equation calculations Implementation involves the calculation of post collisional values for the distribution functions at each lattice point at time t For the generalized form of the Lattice Boltzmann Equation with the collision o
32. DPD The latter is greyed out if not in use If required the electrostatic cutoff re for short range electrostatic interactions and the surface cutoff zc for interactions between particles and solid walls The size of each time step At for integrating the equations of motion and the total steps required for the DPD simulation The number of time steps required to equilibrate the system equilibration steps and the number of time steps to store system variables for rolling averages stack interval The numbers of time steps between rescaling of particle velocities to the desired system temperature during equilibration temp scaling interval and printing summaries in the OUTPUT file print interval The former can be set to zero if no temperature rescaling is required The numbers of time steps between saves of trajectory data to HISTORY files for later visualization save interval and between outputs of statistical data system energy potential energies pressure temperature etc to a plottable CORREL file These values can respectively be set to zero if no trajectory or plot files are required The job time is the maximum real time that can be spent carrying out the DPD simulation the close time gives the time needed to write restart files and shut down the calculation in a controlled manner The restart key for the simulation this can either be set to none for a new simulation a full restart to continue a previ
33. FUNCTIONS 73 e Arguments filename input array of characters e Comments The default output file name is lbout vtk This can be changed by specifying an output file name when calling the routine fOutputLegacy VTKTIncom e Header records int fOutputLegacyVTKTIncom const char filename 1bout e Function Outputs macroscopic temperature and incompressible fluid velocity at each lattice point for system in Structured Grid Legacy VTK format e Dependencies fOutputLegacyVTKT2DIncom fOutputLegacyVIKT3DIncom e Arguments filename input array of characters e Comments The default output file name is lbout vtk This can be changed by specifying an output file name when calling the routine 6 2 8 IbpIOVTK fOutputVTK e Header records int fOutputVTK const char filename lbout int iprop 0 e Function Outputs macroscopic mass density and velocity speeds along x y and z directions at each lattice point of compressible fluid iprop for system in Structured Grid XML VTK format e Dependencies fOutputVTK2D fOutputVTK3D e Arguments filename input array of characters e Comments The default output file name is lbout vts and the density of fluid 0 is output by default This can be changed by specifying an output file name and fluid number up to 1bsy nf 1 when calling the routine fOutputVTKIncom e Header records int fOutputVTKIncom const char filename l1bout int iprop 0 74 CHAPTER 6 DL_MESO_LB
34. For a double precision scalar the alternative SUBROUTINE msg_send_sca_blocked msgtag buf length pe is available msg_wait e Header records SUBROUTINE msg_wait request e Function Causes process to wait for an unblocked message Dependencies None e Arguments request input integer mynode e Header records INTEGER FUNCTION mynode e Function Returns name of process e Dependencies None e Arguments mynode output integer e Comments The serial version of this function returns 0 numnodes e Header records INTEGER FUNCTION numnodes e Function Returns total number of processes available e Dependencies None e Arguments numnodes output integer e Comments The serial version of this function returns 1 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 147 timchk e Header records SUBROUTINE timchk ktim time e Function Determines the time elapsed since the start of the calculation run and if ktim gt 0 prints the time to the OUTPUT file e Dependencies None e Arguments ktim input integer time output real KIND dp e Comments The serial version of DL_ MESO_DPD uses the generic SYSTEM_CLOCK call while the parallel version uses MPI_wtime 11 2 4 error module This module is used to print error messages and shut down DL_MESO_DPD in a controlled manner It requires the modules constants and comms_module to be loaded beforehand error e Header records
35. Licence Agreement provided that 2 1 1 the Licensee may not distribute any of the DL_MESO Software or any Derived Work to any third party or share their use with any third party whether free of charge or otherwise and the Licensee may not sub license the use of any of the DL_MESO Software to any third party unless in each case that third party has complied with clause 2 3 below 2 1 2 the Licensee may not use the DL_MESO Software on behalf of or for the benefit of any third party including without limitation using it to provide bureau outsourcing or application services or facilities management services party unless that third party has complied with clause 2 3 below and 2 1 3 the DL_MESO Software and any Derived Work may be used by the Licensee and its employees and registered students for Academic Purposes only 2 2 If the Licensee wishes to use the DL_MESO Software or any Derived Work in any way except for Academic Purposes or wishes to distribute or make the DL_MESO Software or any Derived Work 2 3 2 4 195 available to any third party for non Academic Purposes it must obtain a commercial licence from STFC STFC may refuse to grant the Licensee a commercial licence If STFC agrees to grant a commercial licence that licence will be on such terms and conditions as STFC sees fit If the Licensee wishes to carry out any collaboration for Academic Purposes with any third party and that third party needs to
36. Peters thermostat itype 3 The Peters thermostat 38 is a modification of the Lowe Andersen thermostat that reduces to standard DPD as the time step tends to zero After integrating all forces using the Velocity Verlet scheme all particle pairs have their velocities modified in a random order using 1 B H ass 645 Tig At bigi VAR i 10 9 1 Bj Tj ass i Ty At bygi WAL s 10 10 where the coefficients as and bij are chosen so that aij At di 2kpT as a2 J i J 2 To ensure that the thermostat both reduces to the DPD thermostat as the time step reduces to zero and is not restricted by the choice of time step the coefficients are chosen as follows lij yijw rij At iran kBT hi 2yijw rij At bij y A 1 exp e 10 12 A correction to the virial of ai ij Tij bijCij v At r is also applied for each particle pair This thermostat can be used with larger time steps than the standard DPD thermostat but with similarly low system viscosities As for the Lowe Andersen thermostat its implementation in parallel running uses a replicated data strategy to carry out the velocity corrections which requires additional memory per processing unit for storing the velocities of all particles in the system and the data required to modify them The efficiency of the Peters thermostat therefore depends upon the total number of particles in the system since all particle pairs ar
37. a z is the distance between the particle and the wall and ze the surface repulsion range The repulsion is applied at the same time as all pairwise forces using the routine hardreflect Since the walls are assumed to be non porous no interactions across them are included and thus boundary halos adjacent to walls are not used The repulsive force magnitudes between the boundary and each species are specified in the CONTROL file at the end of the line with the keyword SURFACE in order of definition Chapter 11 DL MESO DPD Package Reference 11 1 Overview DL_MESO_DPD consists of seventeen Fortran90 modules which should be compiled in the following order prior to the main program itself e constants Contains the constants and parameters required by DL_MESO_DPD e variables Contains the variables and arrays required by DL_MESO_DPD e numeric_container Contains random number generators and a general purpose scalar sum function e comms_module Contains all subroutines necessary for parallel computation e error module Contains subroutines to print error messages and close down DL_MESO_DPD in a controlled manner e parse_utils Contains functions to read in text data from input files e surface_module Contains subroutines for applying boundary conditions at system planes e g solid walls e ewald module Contains subroutines for calculating forces due to electrostatic interactions using Ewald summation based methods
38. a calculational precision of e is required for the Ewald sum the required value of a is equal to wale The long range term for the Ewald sum requires the reciprocal vector k which for an orthogonal simulation box of dimensions L x Ly x Lz is given by Qrky Qrke Ly 27k3 Lz w II where k ka and kg are integers positive and negative from zero to large values specified by the user in the CONTROL file using the keyword KMAX as kj k3 and ky respectively for x y and z dimensions the maximum reciprocal vector ae using the maximum k values can also be defined The long range potential energy term for the entire system is given by kK k2 2 U OV 5 p2 5 qj exp ik 7 10 55 kA0 j where is the imaginary constant y 1 Differentiation of the potential gives the long range electrostatic force on a particle ABI Tq en a exp E r Ep X ikexp ik 7 Y an exp ik Fa 10 56 k0 n The self energy correction term is constant for a given system for all time steps and is equal to Ta U 4 10 57 472 i 134 CHAPTER 10 DL MESO DPD FEATURES while the system charge correction is 2 cc E Vie aay a 10 58 No additional forces are required for these terms In DL MESO DPD the module ewald module contains the most important routines for calculating electrostatic interactions using the algorithm described above The real space short range component is calculated in rou tine ewald_re
39. and a b c and d are system wide parameters whose meanings for each thermostat are given in Table 10 1 Table 10 1 Thermostat coefficients used in DL_MESO_DPD SJ o O a itype 0 MD VV 1 DPD VV 2 Lowe Andersen 3 Peters 4 Stoyanov Groot H2 AR 2 1 i 1 10 2 1 DPD thermostat with standard Velocity Verlet integration MD VV itype 0 This algorithm uses the drag dissipative and random forces FR and pR respectively as described in Chapter 8 as the system thermostat These forces are combined with all other forces between particles pairwise 1This is the first referenced particle in stretching bonds and the second for bond angles and dihedrals 2These values can also be specified using the keyword GAMMA in the CONTROL file the maximum of the two values is used 10 2 THERMOSTATS AND INTEGRATION ALGORITHMS 125 conservative standard and or density dependent bonding electrostatic planar surface external body forces and integrated using the standard Velocity Verlet integrator The combination of the DPD thermostat with the standard MD type VV algorithm is the simplest and least time consuming thermostating algorithm available in DL_MESO_DPD If no integrator thermostat type is selected in the CONTROL file DL MESO_DPD will use this algorithm by default The drag force does however depend upon particle velocities and is therefore only approximated using the mid step values
40. applied in each time step the conditions invoked in this routine are applied to post collisional distribution functions before propagation takes place fPostCollBoundaryIncom e Header records int fPostCollBoundaryIncom e Function Calculates the particle distribution function at different boundaries after the collision step prior to prop agation for incompressible fluids e Dependencies Many for different lattice schemes e Comments Algorithms for other boundary conditions can be added by the user although care should be taken as to when they are applied in each time step the conditions invoked in this routine are applied to post collisional distribution functions before propagation takes place 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 83 fPostPropBoundary Header records int fPostPropBoundary Function Calculates the particle distribution function at different boundaries after propagation prior to the next collision step for compressible fluids Dependencies Many for different lattice schemes Comments Algorithms for other boundary conditions can be added by the user although care should be taken as to when they are applied in each time step the conditions invoked in this routine apply to distribution functions after propagation fPostPropBoundaryIncom Header records int fPostPropBoundaryIncom Function Calculates the particle distribution function at different boundaries after propagati
41. are a number of subdirectories e LBE containing the LBE source code e DPD containing the DPD source code 1 5 DISCLAIMER 3 e JAVA containing the GUI source code e MAN containing the DL_MESO user manual e DEMO containing test cases for DL_MESO e WORK an example working directory 1 5 Disclaimer Neither STFC CCP5 nor any of the authors of the DL_MESO package guarantee that the package is free from error Neither do they accept responsibility for any loss or damage that results from its use 1 6 Copyright STFC Daresbury Laboratory 2011 1 7 Authors Dr Michael Seaton and Prof William Smith Computational Science amp Engineering Department STFC Daresbury Laboratory Warrington WA4 4AD United Kingdom 1 8 Suggestions and Bug Reports We encourage users to send suggestions for improvements and new features for DL_MESO including bug reports and subroutines as well as any additional test cases that demonstrate its features All of these should be sent to michael seaton stfc ac uk Chapter 2 The DL MESO GUI 2 1 Getting Started with the DL_MESO GUI The DL_MESO GUI offers a convenient way of using the DL_MESO package although it is not an essential tool for those who prefer command line operation Appendix A provides details on compiling the DL_MESO program codes manually Working with the GUI requires the availability of Java tools particularly the javac compiler and the java runner for Java 2 v
42. available for the D3Q15 and D3Q19 lattices there is currently no MRT scheme available for D3Q27 The MRT schemes without source terms can be applied using fCollisionMRT while the schemes with Guo like moment transformed source terms can be invoked using fCollisionMRTGuo 5 1 2 Propagation algorithms The simplest implementation the two lattice algorithm involves the use of a temporary array 1bft to copy post collisional distribution functions to their new positions which are subsequently copied back to the main distribution function array 1bf While this particular method is clear easy to understand and can be applied throughout the system s lattice points in any order its drawbacks include the use of two loops for propagation 44 CHAPTER 5 DL_MESO_LBE FEATURES and array copying two large arrays for the distribution functions at each lattice node and significant amounts of time expended in memory access If the user wishes to use this method the routine Propagation is available An alternative more efficient implementation of propagation is the swap algorithm 34 whereby this process is carried out by systematic swapping of pairs of collided distribution function values To make this easier to implement the lattice links are organised so that the conjugate link j to link i e is equal to a m for i 1 m Looping i between 1 and m the post collisional distribution functions for each lattice point f 7 are
43. be changed by specifying an output file name and solute number up to 1bsy ns 1 when calling the routine fOutputQCBIncom e Header records int fOutputQCBIncom const char filename 1bout int iprop 0 e Function Outputs concentration of solute iprop and incompressible fluid velocity at each lattice point for system in Plot3D format e Dependencies fOutputQCB2DIncom fOutputQCB3DIncom e Arguments filename input array of characters e Comments The default output file name is lbout q and the concentration of solute 0 is output by default These can be changed by specifying an output file name and solute number up to 1bsy ns 1 when calling the routine fOutputQT e Header records int fOutputQT const char filename lbout e Function Outputs macroscopic temperature and compressible fluid velocity at each lattice point for system in Plot3D format e Dependencies fOutputQT2D fOutputQT3D 70 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE e Arguments filename input array of characters e Comments The default output file name is lbout q This can be changed by specifying an output file name when calling the routine fOutputQTIncom e Header records int fOutputQTIncom const char filename lbout e Function Outputs macroscopic temperature and incompressible fluid velocity at each lattice point for system in Plot3D format e Dependencies fOutputQT2DIncom fOutputQT3DIncom e Arguments filenam
44. be omitted and is thus greyed out The noise magnitude only has an effect for initializing multiple phase simulations DL_MESO_LBE may include either a phase field parameter or no phase field parameter for systems with multiple phases the Shan Chen 43 algorithm does not require it and this option is currently disabled The number of fluids phases can be increased if a multiple fluid system is to be studied up to 6 fluids may be modelled in DL_MESO_LBE The parameters and boundary conditions for the fluid s must then be set by clicking the set fluid parameters button see below for more details 8 CHAPTER 2 THE DL_MESO GUI 15 The number of solutes needs changing if solute parameters are required if the number of solutes is greater than zero and up to 6 the number of fluids in the above row must be set to 0 or 1 If used the parameters and boundary conditions for the solutes must be set by clicking the set solute parameters button see below for more details 16 The using temperature scalar box may be clicked yes if thermal systems are to be studied If checked the thermal parameters must be set by clicking the set thermal parameters button see below for more details Once all the data in this window and any pop up windows for fluid solute and thermal parameters are filled in the SAVE button should be clicked to write the 1bin sys file 2 2 1 1 Fluid solute and thermal parameters Examples of the pop up windows for
45. by processor y during parallel running e lbout files produced using multiple processors will require gathering or simultaneous plotting see Ap pendix B for more details 6 2 9 IbpBOUND fNextStep e Header records three cases 3D 2D and coordinate long fNextStep int q int xpos int ypos int zpos long fNextStep int q int xpos int ypos long fNextStep int dx int dy int dz long tpos e Function Finds particle position at the next time step when currently at xpos ypos zpos or tpos and moving along direction q or dx dy dz 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 77 e Dependencies fCppMod e Arguments q input integer dx input integer dy input integer dz input integer xpos input integer ypos input integer zpos input integer tpos input long integer fNextStep output long integer fBounceBackF e Header records int fBounceBackF long tpos e Function Performs an on grid bounce back for the fluid distribution function at the tpos th grid point e Dependencies None e Arguments tpos input long integer e Comments This bounce back boundary condition is carried out using f lbopv i f 1bv i i e populations are exchanged with conjugate values fBounceBackC e Header records int fBounceBackC long tpos e Function Performs an on grid bounce back for the solute distribution function at the tpos th grid point e Dependencies None e Arguments tpos input long intege
46. concentrations at the right boundary The solute concentrations at the front boundary Greyed out for two dimensional systems The solute concentrations at the back boundary Greyed out for two dimensional systems The relaxation time Ts for each solute representing diffusivities If selected for inclusion the required thermal properties are 1 2 The heat relaxation time 7 for the system which represents the thermal diffusivity The initial T K and initial dT dt initial temperature and rate of change in temperature related to heat transfers to from the system for the entire system The temperatures and rates of temperature change at the top and bottom boundaries The temperatures and rates of temperature change at the left and right boundaries The temperatures and rates of temperature change at the front and back boundaries These are greyed out for two dimensional systems After filling in all the required values clicking the relevant save button SAVE F SAVE C or SAVE T will store the data in preparation for writing to the lbin sys input file The cancel buttons CANCEL F CANCEL C and CANCEL T will close the pop ups without saving any values 10 CHAPTER 2 THE DL_MESO GUI g LBE DPD SPH Manual Help atticainaltemanin top periodic boundary condition 4 o uh WN 4 bottom periodic boundary condition eft periodic boundary condition right periodic boundary condition front p
47. double speed double startpos int fGetOneMomentSite double speed int fpos long tpos int fGetOneMomentSite double speed int fpos int xpos int ypos int zpos e Function Calculates the momentum of one of the fluids at the grid point e Dependencies None e Arguments startpos input double pointer fpos input integer tpos input long integer xpos input integer ypos input integer Zpos input integer speed output double precision array e Comments The calculation is based on pa gt fi ia The second and third cases are more readable than the first but also slightly slower fGet TotMomentSite e Header records int fGetTotMomentSite double speed double startpos e Function Calculates the momentum of all fluids at the grid point e Dependencies None e Arguments startpos input double pointer speed output double precision array 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 63 fGet TotMoment Domain e Header records int GetTotMomentDomain double momentum e Function Calculates the momentum of all fluids in the domain e Dependencies fGetTotMomentSite e Arguments momentum output double precision array e Comments This function is mainly used to verify that the domain momentum along each axis is conserved fGetOneDirecSpeedSite e Header records three cases float fGetOneDirecSpeedSite int dire double startpos float fGetOneDirecSpeedSite int dire long tpos float fGetOne
48. double pointer rho input double precision Comments The equilibrium distribution function calculated here is 3 8 9 d 3u c2 i 2c4 2c f wip 14 which is only suitable for square lattices Other lattice models e g FHP 9 would require modification or alternative versions of this routine fGetEquilibriumC Header records int fGetEquilibriumC double feq double v double rho Function Calculates equilibrium distribution function for solute Dependencies None Arguments feq output double pointer v input double pointer rho input double precision Comments The equilibrium distribution function calculated here is 3 U gt wC f a c using the solute concentration C and the velocity of the bulk fluid 20 This subroutine can be changed for other Lattice Boltzmann models e g free energy model 51 88 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE fGetEquilibriumT e Header records int fGetEquilibriumT double feq double v double tem e Function Calculates equilibrium distribution function for temperature e Dependencies None e Arguments feq output double pointer v input double pointer tem input double precision e Comments The equilibrium distribution function calculated here is het wT h i aoe using the velocity of the bulk fluid 20 This subroutine can be changed for other Lattice Boltzmann models e g 16 fGetEquilibriumIncom e Header records int
49. e 1bpGET Contains subroutines to calculate macroscopic quantities e g macroscopic density speed and momentum e lbpI0 Contains subroutines to read parameters and write numerical results for plotting and visualization 1bpI0 Contains subroutines to read input files calculate and write summaries 1bpI0Plot3D Contains subroutines to write calculation output files in Plot3D format lbpIOLegacyVTK Contains subroutines to write calculation output files in Legacy VTK structured grid format lbpIOVTK Contains subroutines to write calculation output files in XML VTK structured grid format e 1bpBOUND Contains subroutines for boundary conditions e g calculating the particle distribution function in a shear boundary e 1bpSUB Contains the most important subroutines for Lattice Boltzmann calculations e g particle propagation and site collision e 1bpMPI Contains all subroutines necessary for parallel computation This package can be left out if the user uses only a single processor workstation or a Windows PC It is recommended that DL_MESO users put self defined subroutines into a package called 1bpUSER so upgrades of DL_MESO will not interfere with their contributions 49 50 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE 6 2 DL_MESO_LBE Subroutines and Functions 6 2 1 main There are two primary versions of the main DL_MESO_LBE program serial slbe cpp and parallel plbe cpp The default listings for th
50. e Dependencies contract_bndtbl contract_angtbl contract_dhdtbl error global_sca_and msg_receive_unblocked msg_receive_sca_blocked msg_send_blocked msg_send_sca_blocked msg_wait 158 e Arguments nlimit mdir mp begin final shove input output input input input input input e Comments Only exists in parallel version of DL_MESO_DPD import e Header records CHAPTER 11 integer integer integer real KIND dp real KIND dp real KIND dp DL_MESO_DPD PACKAGE REFERENCE SUBROUTINE import nlimit mdir mp begin final shove e Function Imports particle properties forces virials from boundary halo for integration schemes that only require one set of forces to be calculated e Dependencies error global_sca_and msg receive unblocked msg receive_sca blocked msg_send_blocked msg_send_sca_blocked msg_wait e Arguments nlimit mdir mp begin final shove input output input input input input input e Comments integer integer integer real KIND dp real KIND dp real KIND dp Only exists in parallel version of DL_MESO_DPD This routine is suitable for force integrators using the standard MD Velocity Verlet algorithm including the Lowe Andersen and Peters thermostats import_dpdvv e Header records SUBROUTINE import_dpdvv nlimit mdir mp begin final shove first e Function Imports particle properties forces virials from boundary halo for
51. export nlimit mdir begin final shove Serial version 160 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE e Function Exports particle data positions velocity to neighbouring domain as boundary halo e Dependencies error global_sca_and Parallel version only msg_receive_unblocked Parallel version only msg_receive_sca_blocked Parallel version only msg_send_blocked Parallel version only msg send_sca_blocked Parallel version only msg wait Parallel version only e Arguments nlimit input output integer mdir input integer mp input integer begin input real KIND dp final input real KIND dp shove input real KIND dp exportvelocity e Header records SUBROUTINE exportvelocity nlimit mdir mp begin final shove Parallel version SUBROUTINE exportvelocity nlimit mdir begin final Serial version e Function Exports particle velocities to neighbouring domain as boundary halo for recalculation of dissipative forces as required for DPD Velocity Verlet algorithm e Dependencies error global_sca_and Parallel version only msg_receive_unblocked Parallel version only msg_receive_sca_blocked Parallel version only msg send_blocked Parallel version only msg_send_sca_blocked Parallel version only msg_wait Parallel version only e Arguments nlimit input output integer mdir input integer mp input integer begin input real KIND dp final input real KIND dp shove input real KIND dp exportdensity
52. for DPD calculations involving many body density dependent interactions 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 163 exportdata e Header records SUBROUTINE exportdata nlimit e Function Applies export of particle properties positions velocities to boundary halos of neighbouring domains and or across periodic boundary conditions e Dependencies export e Arguments nlimit input output integer e Comments Non periodic boundaries do not require boundary halos these are excluded from the export of particle properties to prevent interactions between particles near opposite boundaries exportvelocitydata e Header records SUBROUTINE exportvelocitydata nlimit e Function Applies export of particle velocities to boundary halos of neighbouring domains and or across periodic boundary conditions for recalculation of dissipative forces e Dependencies exportvelocity e Arguments nlimit input output integer e Comments Non periodic boundaries do not require boundary halos these are excluded from the export of particle velocities exportdensitydata e Header records SUBROUTINE exportdensitydata nlimit e Function Applies export of local densities for particles to boundary halos of neighbouring domains and or across periodic boundary conditions for calculation of many body DPD interaction forces e Dependencies exportdensity e Arguments nlimit input output integer e Comments Non period
53. for further information fStartMPI e Header records int fStartMPI int argc char argv e Function Starts Message Passing Interface MPI e Dependencies None 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS fCloseMPI e Header records int fCloseMPI O e Function Closes Message Passing Interface MPI e Dependencies None fGlobal Value e Header records six cases int fGlobalValue double vqua int nnum int fGlobalValue int vqua int nnum int vtot int fGlobalValue int vqua int nnum int fGlobalValue long int vqua int nnum int fGlobalValue long int vqua int nnum long int vtot e Function Sums values from all processes and distributes the sum fGlobalProduct e Header records two cases int fGlobalProduct double vqua int nnum int fGlobalProduct int vqua int nnum e Function Multiplies together values from all processors and distributes the product fArrangeProcessors e Header records int fArrangeProcessors e Function Arrange processors according to system dimensions e Comments Calculations are based on lbdm xdim 1bdm ydim J1bdm zdim lbsy nx lbsy ny lbsy nz lbdm xdim x lbdm ydim x lbdm zdim lbdm size fDefineDomain e Header records int fDefineDomain e Function Determines domain parameters for system calculation 99 100 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE fDefineMessage e Header records int fDefineMessage e Fu
54. integration schemes that require two separate sets of forces to be calculated e Dependencies error global_sca_and msg_receive_unblocked 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 159 msg_receive_sca_blocked msg_send_blocked msg_send_sca_blocked msg_wait e Arguments nlimit input output integer mdir input integer mp input integer begin input real KIND dp final input real KIND dp shove input real KIND dp first input logical e Comments Only exists in parallel version of DL_MESO_DPD Setting first to true imports both sets of forces constant and variable while setting to false imports the variable force set and virials This routine is suitable for integrators where types of forces need to be separated and or recalculated e g the DPD Velocity Verlet scheme and Stoyanov Groot thermostat importdensity e Header records SUBROUTINE importdensity nlimit mdir mp begin final e Function Imports local densities from boundary halo e Dependencies error global_sca_and msg receive unblocked msg receive_sca blocked msg_send_blocked msg_send_sca_blocked msg_wait e Arguments nlimit input output integer mdir input integer mp input integer begin input real KIND dp final input real KIND dp shove input real KIND dp e Comments Only exists in parallel version of DL_MESO_DPD export e Header records SUBROUTINE export nlimit mdir mp begin final shove Parallel version SUBROUTINE
55. integrator DPD MD VV is the default but recalculation of dissipative forces at the end of each step DPD DPD VV 10 the Lowe Andersen 27 Peters 38 and Stoyanov Groot 49 thermostats can also be selected Values of y for the DPD and Peters thermostats and T for the Lowe Andersen and Stoyanov Groot thermostats can be specified elsewhere for the entire system individual species or between unlike species but an additional parameter for the Stoyanov Groot thermostat should be set by clicking on set thermostat see below for more details 2 3 DISSIPATIVE PARTICLE DYNAMICS AND THE DL_MESO GUI 13 L_MESO LBE DPD SPH Manual Help Exit E o E job header CONTROL file generated by DL MESO GUI Ij Dissipative Particle Dynamics PA xefi250 0 Xx E 2 temperature fio pressure 3 dissipative factor 1 0 boundary halo 0 0 4 interaction cutoff 1 0 many body cutoff 1 0 5 electrostatic cutoff k surface cutoff N 6 time step 0 01 total steps ooo 7 equilibration steps 000 stack interval 10 8 temp scaling interval fo print interval 10 9 save interval 10 plot file interval 10 10 job time s 3600 0 close time fioo 11 number of species 1 Figure 2 5 Define DPD System 16 The system barostat no barostat is used by default but Langevin 22 and Berendsen 2 barostats are available in combination with all five thermostats If either barostat is selected its paramete
56. interaction length does not have to equal the interaction cut off radius re which applies for dissipative and random force interactions For all thermostats the drag coefficient y or collision frequency I is by default considered to be constant for all species and interactions between species Individual values for the drag coefficient or collision frequency between the same and different species Yag or Pag can be specified in the CONTROL file to override this default Four types of pairwise interactions between particles are available in DL_MESO_DPD Lennard Jones Weeks Chandler Andersen Groot Warren standard DPD and many body density dependent DPD The keyword KTYPE in the CONTROL file determines which interactions are used 10 4 1 Lennard Jones ktype 0 The Lennard Jones potential 23 is a mathematically simple model that approximates interactions both attrac tive and repulsive between pairs of neutral atoms or molecules U rij 4 2 2 10 32 where is the depth of the potential well and oj is the finite distance at which the potential is zero between particles and j in DL_MESO_DPD these are represented by A and B respectively This potential and its related force are calculated for all interparticle distances r less than the interaction cutoff radius re 130 CHAPTER 10 DL MESO_DPD FEATURES Long range system wide corrections to the potential and virial are required MA gt a a U oy
57. of processors along x axis npy number of processors along y axis npz number of processors along z axis volm volume of system dimx size of system in x dimension dimy size of system in y dimension dimz size of system in z dimension sidex size of domain in x dimension sidey size of domain in y dimension sidez size of domain in z dimension delx absolute x coordinate of domain origin dely absolute y coordinate of domain origin delz absolute z coordinate of domain origin nlx number of link cells in domain along z axis nly number of link cells in domain along y axis nlz number of link cells in domain along z axis wdthx link cell size in x dimension wdthy link cell size in y dimension wdthz link cell size in z dimension Table 9 3 Neighbour information parameter meaning map k processor name of neighbour k k 1 left neighbour k 2 right neighbour k 3 lower neighbour k 4 upper neighbour k 5 back neighbour k 6 front neighbour can also be found in pot i and vr1 i The particles modelled by a particular processor have both local and global identity numbers which are stored in the integer arrays loc i and lab i respectively Numbers representing species and molecule types are stored in arrays 1tp i and 1tm i which are used to assign particle masses and the names of particles and molecules respectively to the arrays weight i atmnam i and molnam i Bonded particles are listed by global identity numb
58. or not The barostat can be selected in the CONTROL file using the keyword BARO btype is an integer specifying the ther mostat integration algorithm and a b c and d are system wide parameters whose meanings for each thermostat are given in Table 10 2 Table 10 2 Barostat coefficients used in DL MESO_DPD btype ab cld 1 Langevin Tp yp 2 Berendsen 10 3 1 Langevin barostat btype 1 The governing equation for the Langevin barostat on an orthorhombic simulation cell 22 is the force exerted by the piston expressed as the time derivative of its momentum Pg a WgUg a 1 Dga V Pa Po X miv YpPg a OpLp a 10 17 N 4 where Ny is the number of degrees of freedom for a three dimensional box containing N moving particles Ny 3 N 1 yp and op y 27 W kpT are respectively the drag and random coefficients for the piston and Cp a is a Gaussian random number for dimension a this is set to the same value for all three dimensions if operating isotropically When both y and 9 are set to zero the Langevin barostat reduces to the extended system method The subsequent simulation cell size La can be determined by _ PgaLla Wg La ML 10 18 with the barostat mass W chosen to be equal to iN kpTt where 7 is the characteristic barostat time and 10 should be set equal to between 2 and 12 Yp Yp 128 CHAPTER 10 DL MESO DPD FEATURES The velocitie
59. particle A value unlikei integer A mxpot species number j for unlike particle A value unlikej integer A mxpot integrator thermostat type itype integer thermostat parameter a atherm real KIND dp thermostat parameter b btherm real KIND dp thermostat parameter c ctherm real KIND dp thermostat parameter d dtherm real KIND dp dissipative coefficient y collision frequency T gamma real KIND dp A mxpot random force parameter o sigma real KIND dp A mxpot Lowe Andersen probability of velocity rescaling P At gammala real KIND dp A mxpot Peters dissipative coefficient y gammap real KIND dp A mxpot Stoyanov Groot probability of velocity rescaling P At gammasg real KIND dp A mxpot Stoyanov Groot Nos Hoover coupling parameter a alphasg real KIND dp instantaneous temperature nominator accumulator itpnm real KIND dp instantaneous temperature denominator accumulator itpdm real KIND dp thermostat pair list particle i pairlsti integer A mxpair thermostat pair list particle j pairlstj integer A mxpair thermostat pair list species of particle i pairspei integer A mxpair thermostat pair list species of particle j pairspei integer A mxpair thermostat pair list Maxwell distributed velocity v2 veleq real KIND dp A mxpair thermostat pair list x component of interparticle distance dxx real KIND dp A mxpair thermostat pair list y component of interparticle distance dyy real KIND dp A mxpair thermostat pair list
60. specified value 86 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE fStartDLMESO e Header records int fStartDLMESO O e Function Announces start of DL_MESO_LBE run e Dependencies None e Comments If preferred the call to this routine can be commented out fFinishDLMESO e Header records int fFinishDLMESO e Function Announces end of DL_MESO_LBE run prints simulation time and efficiency measure Millions of Lattice Updates Per Second e Dependencies None e Comments If preferred the call to this routine can be commented out fGet Model e Header records int GetModel e Function Initializes vectors 1bv lbw and lbopv for lattice model e Dependencies D2Q9 D3Q15 D3Q19 D3Q27 e Comments Parameters are specified according to requested space dimension and number of discrete velocities fMarkBoundArea e Header records int fMarkBoundArea e Function Denotes where boundary areas for message passing and or periodic boundary conditions are located e Dependencies int fMarkBoundArea3D int fMarkBoundArea2D 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 87 Comments Only used when boundary areas are used primarily for parallel computing fGetEquilibriumF Header records int fGetEquilibriumF double feq double v double rho Function Calculates equilibrium distribution function for compressible fluid Dependencies None Arguments feq output double pointer v input
61. standard names The modules manybody_module and bond_module may be modified by users to incorporate alternative many body DPD schemes or additional bond angle and dihedral models For anything else however we recommend that DL_MESO users put any self defined subroutines and functions into a module file called user_module 90 so future upgrades of DL_MESO will not interfere with their contributions 11 2 DL_MESO_DPD Subroutines and Functions 11 2 1 Main program dlmesodpd Both the serial and parallel versions of the program operate in a similar way Before DPD calculations start the following tasks are carried out e For parallel running only MPI is started up initcomms and the node properties are determined e An I O channel for the general OUTPUT file is opened e The system clock is consulted for a start time timchk e The starting banner for DL_MESO_DPD is printed e Initial values are set zero e System and bond data are read in and initialized sysdef e The starting configuration is set up start either from scratch or specified by the user using a CONFIG file e The system clock is consulted again for the start of the DPD calculation cycle timchk A loop for DPD calculations is then called from run_module depending on the integrator and barostat selected by the user Each step of DPD calculations involves the following 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 141 e The step counter nstep is increased
62. thermal model for the lattice Boltzmann method in incompressible limit Journal of Computational Physics 146 1 282 300 1998 Xiaoyi He and Li Shi Luo Lattice Boltzmann model for the incompressible Navier Stokes equation Journal of Statistical Physics 88 3 4 927 944 1997 F J Higuera and J Jim nez Boltzmann approach to lattice gas simulations EPL Europhysics Letters 9 7 663 668 1989 R W Hockney and J W Eastwood Computer simulation using particles McGraw Hill International 1981 Takaji Inamuro Masato Yoshino Hiroshi Inoue Riki Mizuno and Fumimaru Ogino A lattice Boltzmann method for a binary miscible fluid mixture and its application to a heat transfer problem Journal of Computational Physics 179 1 201 215 2002 Takaji Inamuro Masato Yoshino and Fumimaru Ogino A non slip boundary condition for lattice Boltz mann simulations Physics of Fluids 7 12 2928 2930 1995 Ask F Jakobsen Constant pressure and constant surface tension simulations in dissipative particle dy namics Journal of Chemical Physics 122 12 124901 2005 J E Jones On the determination of molecular fields II From the equation of state of a gas Proceedings of the Royal Society of London Series A 106 738 463 477 October 1924 Simon Jury Peter Bladon Mike Cates Sujata Krishna Maarten Hagen Noel Ruddock and Patrick Warren Simulation of amphiphilic mesophases using dissipative particle dynamics Physical Chemistry Chemic
63. those expressly set out in this Licence Agreement The Licensee waives any claim for breach of or any right to rescind this Licence Agreement in respect of any representation which is not an express provision of this Licence Agreement However this clause does not exclude any liability which STFC may have to the Licensee or any right which the Licensee may have to rescind this Licence Agreement in respect of any fraudulent misrepresentation or fraudulent concealment before the signing of this Licence Agreement Amendments No variation of or amendment to this Licence Agreement will be effective unless it is made in writing and Signed by each party s representative Third parties No one who is not a party to this Licence Agreement has any right to prevent the amendment of this Licence Agreement or its termination and no one except a party to this Licence Agreement may enforce any benefit conferred by this 6 8 199 Licence Agreement unless this Licence Agreement expressly provides otherwise Governing law This Licence Agreement is governed by and is to be construed in accordance with English law The English Courts will have exclusive jurisdiction to deal with any dispute which has arisen or may arise out of or in connection with this Licence Agreement except that STFC may bring proceedings against the Licensee or for an injunction in any jurisdiction Bibliography 10 11 12 13 14
64. to 1bsy nf 1 when calling the routine fOutputLegacy VTKCA e Header records int fOutputLegacyVTKCA const char filename lbout int iprop 0 e Function Outputs mass fraction and speeds along x y and z directions of compressible fluid iprop at each lattice point for system in Structured Grid Legacy VTK format e Dependencies fOutputLegacyVTKCA2D fOutputLegacyVTKCA3D e Arguments filename input array of characters e Comments The default output file name is lbout vtk and the mass fraction of fluid O is output by default These can be changed by specifying an output file name and fluid number up to 1bsy nf 1 when calling the routine fOutputLegacyVTKCAIncom e Header records int fOutputLegacyVTKCAIncom const char filename lbout int iprop 0 e Function Outputs mass fraction and speeds along x y and z directions of incompressible fluid iprop at each lattice point for system in Structured Grid Legacy VTK format e Dependencies fOutputLegacyVTKCA2DIncom fOutputLegacyVIKCA3DIncom e Arguments filename input array of characters e Comments The default output file name is lbout vtk and the mass fraction of fluid 0 is output by default These can be changed by specifying an output file name and fluid number up to 1bsy nf 1 when calling the routine 72 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE fOutputLegacyVTKCB e Header records int fOutputLegacyVTKCB const char filename lbout int iprop 0
65. too many processes due to either excessive communications or inability to insert molecules into the system 179 180 CHAPTER 12 DL MESO DPD EXAMPLES a t 1 b t 10 t 25 t 50 Figure 12 2 Visualizations of DPD PhaseSeparation test case et for particle type blue for type 2 12 3 Aggregate This serial simulation example consists of 3000 unbonded particles and 30 molecules of 10 particles each with harmonic bonds between them The unbonded particles and molecules are made up of different species with a higher energy parameter for unlike particle interactions This causes the molecules to aggregate which can be seen in Figure 12 3 a snapshot of the simulation x Figure 12 3 Visualization of system at t 2010 from DPD Aggregate test case 12 4 Polyelectrolyte This serial simulation example consists of a slightly hydrophobic polyelectrolyte molecule of 50 particles each with a relative charge of 0 5 immersed in a salt solution of concentration 0 14M 11 The salt solution consists of 9900 neutral water particles 75 cationic salt particles net charge 1 75 anions of charge 1 and 25 counterions of charge 1 to keep the system neutral A similar simulation is included with the polyelectrolyte replaced with a neutral polymer of the same number of particles and the counterions replaced with water CONTROL neutral this should be renamed to CONTROL to run the simulation and used with the same FIELD and MOLECULE files Fi
66. types can be used by default in DL_MESO_DPD exportconfig exportconfig is a utility written in Fortran90 to produce a configuration file in DL_POLY format CONFIG from DL_MESO_DPD restart files export which can be used as a starting point for new simulations including restarting simulations on different numbers of processes Since a limited amount of data is included in restart files the CONTROL file for the simulation is needed to provide some additional information The source code for this utility exportconfig f90 can be used for export or export files created by both the serial and parallel versions of DL_MESO_DPD If the available Fortran90 compiler is invoked by the command 90 the executable config exe can be produced by typing e 90 o config exe exportconfig f90 and either run at the command line or by using the GUI in Process DPD Data after entering the number of processes used in the required field and selecting the required CONFIG file key in the pulldown list This utility can be run with two command line arguments the first indicating the number of processes used to generate the restart data and the second denoting the CONFIG file key levcfg 0 positions only 1 positions and velocities 2 positions velocities and force e g if 16 processes were used and the particle positions and velocities are required either of the following commands can be used e export exe 16 1 e export exe 16 1 If no command l
67. use the DL_MESO Software in connection with that collaboration or if the Licensee wishes to make the DL_MESO Software available online to any third party for Academic Purposes the Licensee must direct that third party to the DL_MESO Website That third party may use the DL_MESO Software and any Derived Work whether obtained from STFC or from the Licensee only if it has completed the registration process on the DL_MESO Website and agreed to the terms and conditions of the Licence Agreement for the use of the DL_MESO Software for Academic Purposes that then appear on the DL_MESO Website This Licence Agreement allows the Licensee to use only the release or version of the DL_MESO Software downloaded by the Licensee from the DL_MESO Website immediately after the Licensee agrees to the terms and conditions of this Licence Agreement the Licensee must acquire a new licence for any future release or version of the DL_MESO Software that it wishes to use The Licensee will not tamper with or remove any copyright or other proprietary notice or any disclaimer that appears on or in any part of the DL_MESO Software and will reproduce the same in all copies of any of the DL_MESO Software and in all Derived Works 3 WARRANTIES AND LIABILITY 3 1 3 2 The DL_MESO Software is provided for Academic Purposes free of charge Therefore STFC gives no warranty and makes no representation in relation to the DL_MESO Software or any assistance or advice th
68. using single sets of particle forces i e standard Velocity Verlet Lowe Andersen Peters 162 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE importdata_dpdvv e Header records SUBROUTINE importdata_dpdvv nlimit e Function Applies import of particle forces constant and variable from boundary halos of neighbouring domains and or across periodic boundary conditions e Dependencies import_dpdvv Parallel version only e Arguments nlimit input output integer e Comments This subroutine is applicable for integrators thermostats using two sets of particle forces i e DPD Velocity Verlet Stoyanov Groot importforcedata e Header records SUBROUTINE importforcedata nlimit e Function Applies import of variable particle forces and virials from boundary halos of neighbouring domains and or across periodic boundary conditions e Dependencies import_dpdvv Parallel version only e Arguments nlimit input output integer e Comments This subroutine is used for importing recalculated forces for particles while using the DPD Velocity Verlet algorithm importdensitydata e Header records SUBROUTINE importdensitydata nlimit e Function Applies import of local densities for particles from boundary halos of neighbouring domains and or across periodic boundary conditions e Dependencies importdensity Parallel version only e Arguments nlimit input output integer e Comments This subroutine is essential
69. 3 VCTPSL VCTPSB VCTCCDLB VCTCCTLF VCTCETR VCTCEDR VCTCEDF VCTCELB VBBPST VBBPSR VBBCCTRB VBBCCDRB VBBCCDLF VBBCETL VBBCETF VBBCERF VBBCEDB VCBPSD VCBPSF VCBCCTLB VCBCCTRE VCBCCDRF VCBCEDL VCBCELF VCBCETB VCBCERB VBTPSL VBTPSB VBTCCDLB VBTCCTLF VBTCETR VBTCEDR VBTCEDF VBTCELB PCTPST PCTPSR PCTCCTRB PCTCCDRB PCTCCDLF PCTCETL PCTCETF PCTCERF PCTCEDB PBTPSD PBTPSF PBTCCTLB PBTCCTRF PBTCCDRF 123 126 129 132 143 146 149 152 221 224 227 230 233 244 247 250 253 322 325 328 331 334 345 348 351 304 423 426 429 432 443 446 449 452 521 524 527 530 533 544 547 550 553 622 625 628 631 634 4 2 DATA STRUCTURE 35 Table 4 3 Notation of boundary condition continued PBTCETR 643 PBTCETL 644 PBTCEDL 645 PBTCEDR 646 PBTCETF 647 PBTCELF 648 PBTCEDF 649 PBTCERF 650 PBTCETB 651 PBTCELB 652 PBTCEDB 653 PBTCERB 654 PCBPST 721 PCBPSD 722 PCBPSL 723 PCBPSR 724 PCBPSF 725 PCBPSB 726 PCBCCTRB 727 PCBCCTLB 728 PCBCCDLB 729 PCBCCDRB 730 PCBCCTRF 731 PCBCCTLF 732 PCBCCDLF 733 PCBCCDRF 734 PCBCETR 743 PCBCETL 744 PCBCEDL 745 PCBCEDR 746 PCBCETF 747 PCBCELF 748 PCBCEDF 749 PCBCERF 750 PCBCETB 751 PCBCELB 752 PCBCEDB 753 PCBCERB 754 PBBPST 821 PBBPSD 822 PBBPSL 823 PBBPSR 824 PBBPSF 825 PBBPSB 826 PBBCCTRB 827 PBBCCTLB 828 PBBCCDLB 829 PBBCCDRB 830 PBBCCTRF 831 PBBCCTLF 832 PB
70. ATTICE BOLTZMANN AND THE DL_MESO GUI 9 Body force z axis the z component of body force on each fluid Greyed out for two dimensional systems The initial fluid densities are applied throughout the system and used to initialize LBE calculations The constant fluid density pp for incompressible systems The fluid densities at the top boundary The fluid densities at the bottom boundary The fluid densities at the left boundary The fluid densities at the right boundary The fluid densities at the front boundary Greyed out for two dimensional systems The fluid densities at the back boundary Greyed out for two dimensional systems The relaxation time 7f for each fluid these values should be greater than 0 5 to give non zero kinetic viscosities The bulk relaxation time Tf bulk for each fluid these values are only used in multiple relaxation time schemes to define the bulk viscosity and again should be greater than 0 5 If more than one phase is being modelled non zero interaction parameters between the fluids are required gap for the Shan Chen algorithm Solutes require the following data in the following row numbers Ls 2 3 4 The initial concentrations of the solutes throughout the system as used for initialization The solute concentrations at the top boundary The solute concentrations at the bottom boundary The solute concentrations at the left boundary The solute
71. BCCDLF 833 PBBCCDRF 834 PBBCETR 843 PBBCETL 844 PBBCEDL 845 PBBCEDR 846 PBBCETF 847 PBBCELF 848 PBBCEDF 849 PBBCERF 850 PBBCETB 851 PBBCELB 852 PBBCEDB 853 PPBCERB 854 4 2 3 Storage of running information The Lattice Boltzmann component of DL_MESO defines three structures to store the system information The parameters in these structures are listed in Tables 4 4 4 5 and 4 6 Table 4 4 System information parameter meaning lbsy nd space dimension lbsy nq number of discrete speeds lbsy nf number of fluids lbsy nc number of solutes lbsy nt number of temperature scalars 0 or 1 lbsy pf phase field order parameter lbsy nx number of grid points in x direction lbsy ny number of grid points in y direction lbsy nz number of grid points in z direction 1bsy nz 1 when 1bsy nd 2 1Referring to the physical system being simulated rather than the computer system 36 CHAPTER 4 DL_MESO_LBE BASIC DEFINITION Table 4 5 Domain information parameter meaning lbdm rank name of the processor lbdm size number of processors lbdm bwid domain boundary width set to zero for serial running lbdm xcor x coordinate of the processor lbdm ycor y coordinate of the processor lbdm zcor z coordinate of the processor lbdm zcor 0 when lbsy nd 2 lbdm xdim number of processors along x axis lbdm ydim number of processors along y axis lbdm zdim number of processors alon
72. Comments Collisions are carried out using multiple relaxation time MRT schemes coupled with Guo like forcing terms 39 8 ei a fi T 5 TF Mf 1 315 5 fCalcPotential_ShanChen e Header records int CalcPotential_ShanChen 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 95 e Function Calculates the interaction potential 1 exp p as suggested by the Shan Chen model e Dependencies fGetOneMassSite e Comments The Shan Chen model is detailed in 43 44 Alternative mesoscale interaction potentials can be introduced here fCalcInteraction_ShanChen e Header records int fCalcInteraction ShanChen int xpos int ypos int zpos e Function Calculates particle interaction forces according to Shan Chen model e Arguments xpos input integer ypos input integer zpos input integer e Comments Further details can be found in 43 44 Similarly named routines to calculate interaction forces for alternative mesoscale algorithms can be added by the user fInteractionForceShanChen e Header records int fInteractionForceShanChen e Function Calculates interaction forces for all fluids based on the Shan Chen model 43 44 e Dependencies fCalcPotential_ShanChen fCalcInteraction_ShanChen e Comments Alternative mesoscale interactions can be applied using similar routines fCollisionBGK e Header records int fCollisionBGKO e Function Collision steps for all compressible fluids us
73. D subdirectory These can be divided into problems that should by default be run using the serial and parallel versions of DL MESO_DPD the latter have been tested with 96 processing units but should also run on smaller and larger numbers They consist of the following cases 12 1 Mixture_Ser This serial simulation example consists of 1000 particles with 3 species Species 1 has a population of 333 species 2 has a population of 333 and species 3 has a population of 334 particles All particle types have identical sizes and masses but different energy parameters using the default values for unlike particle parameters Figure 12 1 gives a snapshot of the system at the end of the simulation demonstrating mixing between the three particle types represented by different colours Figure 12 1 Visualization of system at final time step from DPD Mixture Ser test case 12 2 PhaseSeparation This serial simulation example consists of 3000 particles with 2 species Species 1 has a population of 1500 and species 2 has a population of 1500 particles Both particle types have identical sizes masses and like like energy parameters but the unlike energy parameter has been set to a larger value to produce phase separation which can clearly be seen in Figure 12 2 An AVI video file of the first half of the simulation can be found in the subdirectory 1These can also be run using the parallel version of DL_MESO_DPD but care should be taken not to use
74. DL_MESO USER MANUAL M A Seaton and W Smith STFC Daresbury Laboratory Daresbury Warrington Cheshire WA4 4AD United Kingdom Version 2 4 May 2011 Contents Contents DL_MESO General Information Td Description nse sad a Gaye ce hak SE ae ae ee a Ace OTB ES ES 1 2 Functionality xe ows top eb ere A A hod A ae ee Sag eee A BS BORA to ARequirements st et se te lee oles te See ay ee Be eh hk A 1 4 The DL_ MESO Directory Structure 2 0 0 00 2 00000000 0000000000004 db Disclaimer sen io ides shed hh Ate ese eee BR Bs Pa ee Pe tay hed Ged Be ee ee LG Copyright hap 2AM ARAN A ed eaten a Bh Bote EE Pah AA Be oe a Te Authors a ataie cat ta we Be A den Se ee Ae ey Se ee AEA 1 8 Suggestions and Bug Reports a The DL_MESO GUI 2 1 Getting Started with the DL MESO GUI 0 00 00 0200 00 00000 2 2 Lattice Boltzmann and the DL MESO GUI 02 0 0 0 0 0000 0 2 3 Dissipative Particle Dynamics and the DL MESO G l o o o o 2 4 Compiling and running DL MESO 0 000000000000 220 o Notes s gue ek aah Gays a aye ae eee BR Be Beet a eee tht Syd ad Sie ee a Lattice Boltzmann Equation LBE The Lattice Boltzmann Equation Basic Theory 3 1 Introduction e ee ee O A ee EE oe ee Oe a ee et 3 24 Basics Dennitions a ee ood Be oe ee hyena Gok Bee Soe EK eS 3 3 Derivation of Equilibrium 3 4 Structural Relaxation and Macroscopic Equations ooa a 00200 eee
75. DPD SPH Manual Help Exit Dissipative Particle Dynamics operating system 1 Prior compilation of molecule generate utility required interaction bond x add to FIELD 3 type hamoni a parameters a 0 0 c 5 b 0 0 d molecule name MOLI 6 molecule population o C bonds 7 angles dihedrals global Figure 2 7 Set DPD Molecules The type of surface interactions or boundary conditions to be applied Wall directions if the checkbox for a particular dimension is ticked the boundary condition will be applied to the surfaces orthogonal to the specified axis Surface parameters for the hard surfaces with soft repulsions the surface repulsion for each species Awail i is required After filling in all the required values clicking the relevant save button SAVE SP SAVE P SAVE T SAVE B SAVE E or SAVE SF will store the data in preparation for writing to the CONTROL file The cancel buttons CANCEL SP CANCEL P CANCEL T CANCEL B CANCEL E and CANCEL SF will close the pop ups without saving any values 2 3 2 Defining DPD Molecules Figure 2 7 shows the Set DPD Molecules panel with the sections for data entry labelled in red numbering The following data are required 1 The operating system is required before launching a command line terminal and running the molecule generation utility molecule generate cpp to create molecules for the DPD sim
76. DirecSpeedSite int dire int xpos int ypos int zpos e Function Calculates the grid speed for all compressible fluids along direction dire 0 for x 1 for y and 2 for z e Dependencies None e Arguments startpos tpos xpos ypos Zpos dire fGetOneDirecSpeedSite e Comments input input input input input input output double pointer long integer integer integer integer double pointer floating point Mainly used to output grid speed The second and third cases are more readable than the first but also slightly slower fGetOneDirecSpeedIncomSite e Header records three cases float fGetOneDirecSpeedIncomSite int dire double startpos float fGetOneDirecSpeedIncomSite int dire long tpos float fGetOneDirecSpeedIncomSite int dire int xpos int ypos int zpos e Function Calculates the grid speed for all incompressible fluids along direction dire 0 for x 1 for y and 2 for z e Dependencies None 64 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE e Arguments startpos input double pointer tpos input long integer xpos input integer ypos input integer Zpos input integer dire input double pointer fGetOneDirecSpeedIncomSite output floating point e Comments Mainly used to output grid speed The second and third cases are more readable than the first but also slightly slower 6 2 5 IbpIO fDefineSystem e Header records int fDefineSystem const char filename lbin sys e F
77. E PACKAGE REFERENCE Function Outputs macroscopic mass density and velocity speeds along x y and z directions at each lattice point of incompressible fluid iprop for system in Structured Grid XML VTK format Dependencies fOutputVTK2DIncom fOutputVTK3DIncom Arguments filename input array of characters Comments The default output file name is lbout vts and the density of fluid 0 is output by default This can be changed by specifying an output file name and fluid number up to 1bsy nf 1 when calling the routine fOutputVTKCA Header records int fOutputVTKCA const char filename 1lbout int iprop 0 Function Outputs mass fraction and speeds along x y and z directions of compressible fluid iprop at each lattice point for system in Structured Grid XML VTK format Dependencies fOutputVTKCA2D fOutputVTKCA3D Arguments filename input array of characters Comments The default output file name is lbout vts and the mass fraction of fluid 0 is output by default These can be changed by specifying an output file name and fluid number up to 1bsy nf 1 when calling the routine fOutputVTKCAIncom Header records int fOutputVTKCAIncom const char filename l1bout int iprop 0 Function Outputs mass fraction and speeds along x y and z directions of incompressible fluid iprop at each lattice point for system in Structured Grid XML VTK format Dependencies fOutputVTKCA2DIncom fOutputVTKCA3DIncom Argument
78. Krook approximation 4 to describe the structural relaxation The single particle distribution function evolves to the equilibrium state via fil At t At fi Z t a fs Bt 22 3 12 where Ty is called the relaxation time and is related to the kinetic viscosity of fluid To derive the macroscopic equations the left hand side of Equation 3 12 can be expanded as o O At mea fi At t At fi Zt E 0 eiaa fi Zt 3 13 m 1 Expanding the instantaneous particle distribution function around its equilibrium and retaining only the first order gives fi 5 8 f 74 t Tf Ot eia a Zt O 0 3 14 2 Substituting Equations 3 13 and 3 14 into the left hand side of Equation 3 12 gives the second order differential equation for the equilibrium distribution LE L 8 eia a fet ws Os Ciada 17 O 0 i Tf A At where wf Tf gt Summing Equation 3 15 over i and ignoring the second order deriviative we obtain 0 ip OaPUa wf OB onus On 5 fenci 0 0 3 16 Summing Equation 3 16 times e over i we obtain 0 ipua Og 5 fi eiaeig wr Oy a 5 fi eiaCiy 08 dE ieee 0 6 3 17 Equation 3 17 shows that the second term in Equation 3 16 is of the third order in the derivative Therefore we have the continuity equation to the second order of the derivative Defining the third and fourth order moments SS FP eiaeip Pag p
79. L_MESO_LBE BASIC DEFINITION Table 4 1 Boundary condition category value meaning 0 liquid 10 domain boundary 11 inside solid 12 on grid bounce back boundary 13 mid link bounce back boundary 100 199 constant speed composition and temperature boundary 200 299 constant speed Neumann composition and temperature boundary 300 399 constant speed and composition Neumann temperature boundary 400 499 constant speed and temperature Neumann composition boundary 500 599 constant pressure composition and temperature boundary 600 699 constant pressure Neumann composition and temperature boundary 700 799 constant pressure and composition Neumann temperature boundary 800 899 constant pressure and temperature Neumann composition boundary 4 2 Data structure 4 2 1 Storage of particle distribution functions For a system with a square lattice the total number of grid points lbsy nx x lbsy ny x lbsy nz where lbsy nx lbsy ny and lbsy nz are the numbers of grid points along the x y and z axes respectively The grid points are arranged in a serial order of 9000 9001 900 1bsy nz gt 9010 9011 JO 1bsy ny lbsy nz J100 9101 Jlbsy nx lbsy ny lbsy nz At each grid point DL_MESO_LBE arranges the particle distribution functions in order of fluid functions solute functions temperature functions and phase field order parameter For example for a D2Q9 lattice with two fluids scalar temperat
80. Licence Agreement will take effect and the Licence Period will start when the Licensee has agreed to the terms and conditions of this Licence Agreement and downloaded the DL_MESO Software from the DL_MESO Website This Licence Agreement will terminate immediately and automatically if 5 2 1 the Licensee is in breach of this Licence Agreement or 5 2 2 the Licensee becomes insolvent or if an order is made or a resolution is passed for its winding up except voluntarily for the purpose of solvent amalgamation or reconstruction or if an administrator administrative receiver or receiver is appointed over the whole or any part of its assets or if it makes any arrangement with its creditors The Licensee s right to use the DL_MESO Software will cease immediately on the termination of this Licence Agreement and the Licensee will destroy all copies of the DL_MESO Software that it or any of its employees or students then holds Clauses 1 2 2 3 4 5 3 5 4 5 5 and 6 will survive the expiry of the Licence Period and the termination of this Licence Agreement and will continue indefinitely STFC may stop providing any assistance or advice in relation to or any corrections new releases or versions of the DL_MESO and 198 APPENDIX C DL_MESO LICENCE AGREEMENT ACADEMIC PURPOSES may stop updating or publishing the DL_MESO Website at any time 6 GENERAL 6 1 6 2 6 3 6 4 6 6 6 7 Headings The headings i
81. N Table 9 4 DL MESO_DPD Parameters continued function parameter data type notation Cartesian coordinate x for particle XXX real KIND dp A maxdim Cartesian coordinate y for particle yyy real KIND dp A maxdim Cartesian coordinate z for particle zzz real KIND dp A maxdim particle mass weight real KIND dp A maxdim virial of particle vrl real KIND dp A maxdim potential energy of particle pot real KIND dp A maxdim potential energy accumulator pe real KIND dp virial accumulator vir real KIND dp kinetic energy accumulator tke real KIND dp bond potential energy accumulator be real KIND dp angle potential energy accumulator ae real KIND dp dihedral potential energy accumulator de real KIND dp electrostatic potential energy accumulator ee real KIND dp bond length accumulator bdlng real KIND dp bond angle accumulator bdang real KIND dp bond dihedral accumulator bddhd real KIND dp average system potential energy avepe real KIND dp average system virial avevir real KIND dp average system kinetic energy avetke real KIND dp average system total energy avete real KIND dp average system pressure aveprs real KIND dp average system volume avevlm real KIND dp average system temperature avettp real KIND dp average system bond potential energy avebe real KIND dp average system angle potential energy aveae real KIND dp average system dihedral potential energy avede real KIND dp average system electros
82. ND dp A mxpot localized densities rhomb real KIND dp A maxdim mxspe x component of external body force per particle bdfrcx real KIND dp A mxspe y component of external body force per particle bdfrcy real KIND dp A mxspe z component of external body force per particle bdfrcz real KIND dp A mxspe species name for particle atmnam character 8 A maxdim molecule name for particle molnam character 8 A maxdim particle global identity number lab integer A maxdim particle species number ltp integer A maxdim particle molecule type number ltm integer A maxdim particle link cell population number Tet integer A maxdim particle link cell number link integer A maxdim particle local domain cell identity number loc integer A maxdim particle molecule type number ltm integer A maxdim force x component on particle xx real KIND dp A maxdim force y component on particle fyy real KIND dp A maxdim force z component on particle fzz real KIND dp A maxdim variable force z component on particle fvx real KIND dp A maxdim variable force y component on particle fvy real KIND dp A maxdim variable force z component on particle fvz real KIND dp A maxdim velocity xz component of particle VXX real KIND dp A maxdim velocity y component of particle vyy real KIND dp A maxdim velocity z component of particle vzz real KIND dp A maxdim Requirement only applicable if Ewald sum electrostatics are used 122 CHAPTER 9 DL_MESO_DPD BASIC DEFINITIO
83. O L BE cg ispi soa rd blk a WO ad ae ae CONTENTS 109 111 179 CONTENTS 111 A 2 DL_MESO_DPD A ion ie dh AS SD AT Re E A E A ad a 184 B DL_MESO Utilities 187 Bil DEMESO LBE cd Gel aed A A Ag Sea A ig ae ae es e 2 a 187 B2 DE MESOEDED teta cakes Be hap A A ee Ah de Ma Se hea te ale da aie A 188 C DL_MESO Licence Agreement Academic Purposes 193 Bibliography 201 Acknowledgements DL_MESO was developed under the auspices of the Engineering and Physical Sciences Research Council EP SRC for the EPSRC s Collaborative Computational Project for the Computer Simulation of Condensed Phases CCP5 The members of the CCP5 DL_MESO consortium were David M Heyes University of Surrey Chris M Care Sheffield Hallam University Peter V Coveney University College London David Emerson STFC Daresbury Laboratory Rob English North East Wales Institute Andrea Ferrante Unilever Port Sunlight lan Halliday Sheffield Hallam University John Harding University of Sheffield Sebastian Reich Imperial College Bill Smith STFC Daresbury Laboratory Patrick B Warren Unilever Port Sunlight Julia Yeomans Oxford University Many other people have given advice and encouragement in the development of DL_MESO We gratefully acknowledge the support of the following people Maurice Leslie Richard Wain Alexandre Dupuis Jonathan Chin Michael Dupin Weiming Liu and John Purton Particular thanks go to Rongshan Qin at Imperial College as th
84. REFERENCE deportdata error exportvelocitydata DPD VV only dragforces_dpdvv DPD VV only importforcedata DPD VV only global_sca_sca int global_sca_sca_dble global_sum_dble global_sum_int Lowe Peters Stoyanov e Arguments stage input integer e Comments The Langevin barostat is configured to apply the following piston force rate of change of momentum in dimension a A Dga V Pa Po N 2 MiV YpPg a Oplp a with the option of keeping the system isotropic by using equal values of Py and Cp for all dimensions This may be switched off if imposing constant surface tensions at given planar surfaces Rescaling of the simulation volume takes place in the first stage coupled with the integration of particle forces while iteration of the barostat force to achieve convergence of particle velocities to machine precision takes place in the second stage verlet_ _berend e Header records SUBROUTINE verlet_ _berend stage e Function Solves the equations of motion using the Velocity Verlet scheme coupled with the Berendsen barostat 2 e Dependencies hardreflect deportdata error exportvelocitydata DPD VV only dragforces_dpdvv DPD VV only importforcedata DPD VV only global_sca_sca_int global_sum_dble global_sum_int Lowe Peters Stoyanov e Arguments stage input integer e Comments The Berendsen barostat rescales particle positions and system sizes in dimension a by the factor Ha 1 Pap
85. Run DPD Program button which activates a panel that allows you to select the required submission command and then submit the job You may need to create a suitable run script in your working directory beforehand if running the job in parallel 18 2 5 CHAPTER 2 THE DL_MESO GUI Collecting data from multiple processors and processing it for visualization is possible using the Gather LBE Data and Process DPD Data buttons Note that the utilities need to be compiled in the working directory prior to use details on this and their functions can be found in Appendix B The results of LBE and DPD calculations may be plotted using the Plot LBE Results and Plot DPD Results buttons which allows the user to select plotting and visualization applications including those not available in the pull down lists Note that these need to be already installed on the workstation in use before being invoked by the GUI if they require running from a command line tick the run in terminal box before launching the application Notes There are some inactive buttons reserved for later use The GUI does not produce initial state files lbin init for LBE CONFIG for DPD prior to simulations Click EXIT to close down the GUI Part I Lattice Boltzmann Equation LBE Chapter 3 The Lattice Boltzmann Equation Basic Theory 3 1 Introduction The Lattice Boltzmann method utilizes fully discretized space time and velocity to describe the evolution o
86. S 145 e Arguments iii input output array of integers non input integer e Comments For integer scalars the alternative SUBROUTINE global_sca_sum_int iii is available msg_receive_blocked e Header records SUBROUTINE msg_receive_blocked msgtag buf length e Function In a blocking call receives data in the form of a double precision array buf O e Dependencies None e Arguments buf output array of real KIND dp msgtag input integer length input integer e Comments For a double precision scalar the alternative SUBROUTINE msg_receive_sca_blocked msgtag buf length is available msg_receive_unblocked e Header records INTEGER FUNCTION msg_receive_unblocked msgtag buf length e Function In a non blocking call receives data in the form of a double precision array buf e Dependencies None e Arguments buf output array of real KIND dp msgtag input integer length input integer msg_receive_unblocked output integer e Comments Since no scalars are received in non blocking calls no scalar version of this function exists msg_send_blocked e Header records SUBROUTINE msg_send_blocked msgtag buf length pe e Function In a blocking call send data from a double precision array buf e Dependencies None 146 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE e Arguments buf input array of real KIND dp msgtag input integer length input integer pe input integer e Comments
87. SUBROUTINE error idnode iode value e Function Prints user friendly error message in OUTPUT file and closes down DL_MESO_DPD e Dependencies abortcomms e Arguments idnode input integer iode input integer value input integer e Comments The error code iode closes down DL_MESO_DPD when positive negative values can be used to print warning messages for non fatal problems 11 2 5 parse_utils getword e Header records CHARACTER 16 FUNCTION getword txt n e Function Obtains the nth word from a line of text txt separated by spaces or commas 148 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE e Dependencies None e Arguments txt input character 160 n input integer getword output character 16 parseint e Header records INTEGER KIND 1i FUNCTION parseint word e Function Reads the integer contained in the string word e Dependencies None e Arguments word input Ccharacter 16 parseint output integer KIND li parsedble e Header records REAL KIND dp FUNCTION parsedble word e Function Reads the double precision number contained in the string word e Dependencies None e Arguments word input character 16 parsedble output real KIND dp getint e Header records INTEGER KIND 1i FUNCTION getint txt n e Function Reads the nth word of the string txt to obtain an integer e Dependencies parseint e Arguments txt input character 160 n input integer getin
88. ULE file includes the species and relative Cartesian coordinates for the particles with the origin at the centre of the cube in all molecule types up to eight by default each molecule type is formatted thus MOLECULE n cubesize 1 T1 Y 21 S2 T2 Y2 22 Sn n Yn Zn where n is the number of particles in the molecule default maximum mxbonds of 500 cubesize is the size of cube within which the molecule lies and used to insert it into the system which can be no bigger than the smallest dimension for the entire system while the species number and Cartesian coordinates x y z of particle are given by si i yi and z respectively If the user is likely to change the order or number of species in the CONTROL file the names of species as defined in the CONTROL file can be used in place of s The FIELD file contains both the properties of the bonds angles and dihedrals required and the bond pairs bond angle triples and bond dihedral quadruples for each molecule type The bond angle dihedral properties are listed first in the following formats BOND bondtype abcd ANGLE angletype abcd DIHEDRAL dihedraltype abcd SIf any of the values are omitted they are given the values of Ree 11 3 DL_MESO_DPD DATA FILES 175 where the bond angle dihedral types are integers as listed in Section 10 6 and a b cand d are the corresponding parameters double precision floating point numbers as per Tables 10 4 10 5 and 10 6 The relati
89. Y 2 bap NaN peop ES Es 10 33 a B gt a 12012 20 lr a a W oy gt 2 bap IN 99 5 10 34 a Broa where dag is the Kronecker delta 1 when a 8 0 when a 4 and Na the total number of particles of species a These corrections multiplied by the volume are stored in the array clr to eliminate the need to adjust them if the system volume is changed by a barostat 10 4 2 Weeks Chandler Andersen ktype 1 The Weeks Chandler Andersen potential 57 is a modification of the Lennard Jones potential to produce purely repulsive short range interactions ne O U rij 4 e H ij 10 35 ij ij 1 This interaction is applied for interparticle distances up to 260 which should be less than or equal to the interaction cutoff radius re No long range corrections to potential energy or virials are required 10 4 3 Standard DPD ktype 2 The Groot Warren standard form of DPD 13 uses the following purely repulsive soft potential 1 U rij 5 This conservative interaction is applied for interparticle distances up to Bij which should be less than or equal to Te 10 4 4 Many body DPD ktype 3 The conservative force in standard DPD depends only upon the species interacting and the interparticle sepa ration which yields a quadratic equation of state Many body DPD 37 53 is a method of providing alternative thermodynamic behaviours to DPD particles by making conservative forc
90. _LBE This is based upon substituting the unknown distribution functions for a boundary point with local equilibrium values but using an adjusted solute concentration or temperature to produce the correct value for the property at that point This concentration temperature is obtained by substituting the equilibrium distribution function for the unknown populations into the sum of distribution functions equal to the required concentration temperature For example for a top edge with a specified temperature Ty using the D2Q9 lattice scheme TCEDF the adjustment temperature T is given as 6 T 3 To ho ha ha he hr ha w y and the missing populations i 3 4 5 in this case are given by n wr 1 3 E c where Vu is the known or calculated velocity for all fluids at the same boundary DL_MESO_LBE includes implementations of the Inamuro boundary condition for all four lattice schemes Only planar surfaces can be dealt with using this method as for constant pressure boundaries concave edges and corners use the equilibrium distribution function for the required concentration temperature and zero velocity On grid bounce back is used for Neumann zero gradient conditions of solute concentrations and or tempera ture to keep this type of boundary condition entirely local 5 3 Mesoscale interactions DL_MESO_LBE uses the Shan Chen model 43 44 to model interactions between multiple phases and compo
91. a t 10 31 The remainder of the Velocity Verlet algorithm is unchanged although the scaling factor for the beginning of the next time step can be calculated at this point using Equation 10 29 No iteration is required for this barostat 10 4 Particle particle interactions Pairwise particle interaction parameters in DL MESO DPD are specified in the CONTROL file for each species Aa and are stored in the array aaa i Other properties that are also required for each species are its name up to 8 characters stored in namspe i the number of unbonded particles of that species nspec i the particle mass amass i radius B bbb i and charge q ccc i The interaction parameter values and lengthscales are stored in the array vvv k 1 4 in preparation for DPD calculations taking advantage of the fact that Aag Aga by converting the species numbers and j into a potential number k 14 i By default the conservative interaction parameter between different particle species a and is set to the geometric mean i e Aag VA As This choice can be overridden by specifying a replacement value for Aag in the CONTROL file using the keyword UNLIKE By default the interaction lengthscales between species a and is set to the arithmetic mean i e Bag 4 Ba Bg Like the interaction parameters the default lengthscale can be overridden by specifying a re placement value for Byg in the CONTROL file It should be noted that the
92. ade between them in DL_MESO_DPD 10 7 Surface interactions The default boundaries for a simulation box are periodic i e particles leaving the system are replaced at the opposite face with the same velocity Certain systems may require alternative boundary conditions and DL_MESO_DPD can include these at the system boundaries The keyword SURFACE in the CONTROL file can be used to specify the type of surface interaction srftype and which surface s are to be included srfx srfy srfz each set to 0 for periodic boundaries and 1 or greater for other types as well as any necessary parameters The module surface_module includes routines to apply boundary conditions at those surfaces and to identify the processing units that include surfaces surfacenodes 138 CHAPTER 10 DL MESO DPD FEATURES 10 7 1 Hard reflecting boundaries srftype 1 This boundary condition is applied by using a combination of specular free slip reflection at the system boundaries and a soft short range wall repulsion to reduce density oscillations 56 The reflection is achieved by moving any particle leaving the system at a particular boundary back into it and inverting the velocity component normal to the wall but maintaining the tangential momentum which is achieved with the routine surfacebounce The soft short range wall repulsion is given by baa Unautz Awall a 5 5 z lt Ze 10 85 e where Awall a is the repulsive force magnitude with species
93. al using a link cell algorithm with a cutoff radius of re a larger boundary halo than for standard pairwise interactions is therefore required The reciprocal space long range component is calculated using the scheme described by 46 parallelized by distribution over atomic sites which requires global summations but is more efficient in terms of memory usage than distribution of k vectors The routine ewald_reciprocal also adds the self interaction and charged system corrections which are calculated using the routine elecgen in config_module 10 6 Bond interactions between particles Molecules of particles bonded together can be included in calculations using a MOLECULE file to define the relative positions of particles and a FIELD file to define the properties and topologies of the bonds angles and dihedrals between them These data are used in the start subroutine to add the specified numbers of molecules into the system before DPD calculations commence and to create tables listing the particles that are included in bonds angles and dihedrals 10 6 1 Stretching bonds DL_MESO_DPD can model four forms of bond potential and corresponding force between specified particles all of which are functions of the distance between them The available bond forms between particles and j are as follows 1 Harmonic Hookean Fraenkel bond U rij 5 rij ro 10 59 2 Finitely Extendible Non linear Elastic FENE bond
94. al Physics 1 9 2051 2056 1999 J M V A Koelman and P J Hoogerbrugge Dynamic simulations of hard sphere suspensions under steady shear EPL Europhysics Letters 21 3 363 368 1993 Pierre Lallemand and Li Shi Luo Theory of the lattice Boltzmann method Dispersion dissipation isotropy Galilean invariance and stability Physical Review E 61 6 6546 6562 June 2000 C P Lowe An alternative approach to dissipative particle dynamics EPL Europhysics Letters 47 2 145 151 July 1999 John F Marko and Eric D Siggia Stretching DNA Macromolecules 28 26 8759 8770 1995 G Marsaglia and T A Bray A convenient method for generating normal variables SIAM Review 6 3 260 264 July 1964 George Marsaglia Arif Zaman and Wai Wan Tsang Toward a universal random number generator Statistics and Probability Letters 9 1 35 39 January 1990 C A Marsh G Backx and M H Ernst Static and dynamic properties of dissipative particle dynamics Physical Review E 56 2 1676 1691 August 1997 Nicos S Martys and Hudong Chen Simulation of multicomponent fluids in complex three dimensional geometries by the lattice Boltzmann method Physical Review E 53 1 743 750 January 1996 Makoto Matsumoto and Takuji Nishimura Mersenne twister a 623 dimensionally equidistributed uniform pseudo random number generator ACM Transactions on Modeling and Computer Simulation 8 1 3 30 1998 Keijo Mattila Jari Hyv luoma Tuomo Ross
95. arry out 58 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE fGet Tot MassSite e Header records two cases double fGetTotMassSite double startpos double fGetTotMassSite long tpos e Function Calculates the mass density of all fluids at a grid point e Dependencies None e Arguments startpos input double pointer tpos input long integer fGetTotMassSite output double precision e Comments Mass density is calculated according to the definition p 5 fi The second case carries out the same calculation as the first but using tpos for the grid point fGetOneMassDomain e Header records double fGetOneMassDomain int fpos e Function Calculates the total mass of fpos th fluid in the domain e Dependencies None e Arguments fpos input integer fGetOneMassDomain output double precision e Comments The total mass of the domain does not include boundary areas used for message passing or nodes used to apply boundary conditions fGet TotMassDomain e Header records double fGetTotMassDomain e Function Calculates the total mass of all fluids in the domain e Dependencies None e Arguments fGetTotMassDomain output double precision e Comments The total mass of the domain does not include boundary areas used for message passing or nodes used to apply boundary conditions 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 59 fGetFracSite e Header records three cases double fGetFracSite int fpos double sta
96. at STFC may give in connection with the DL_MESO Software The Licensee will indemnify STFC against any and all claims arising as a result of the Licensee having made any of the DL_MESO Software or any Derived Work available to any third party Before using any of the DL_MESO Software the Licensee will check that the DL_MESO Software does not contain any Harmful Element STFC does not warrant that the DL_MESO Software will run without interruption or be error free or be free from any Harmful Element STFC is not obliged to provide any support or error correction service assistance or advice in relation to the DL_MESO Software but the Licensee may access any error corrections and online assistance that STFC chooses to make 196 3 3 3 4 3 6 3 7 3 8 APPENDIX C DL_MESO LICENCE AGREEMENT ACADEMIC PURPOSES available on the DL_MESO Website from time to time If STFC does provide that sort of service assistance or advice subject to clause 3 7 STFC will not be liable for any loss or damage suffered by the Licensee as a result STFC will not be liable to the Licensee to the extent that any loss or damage is caused by the Licensee s failure to implement or the Licensee s delay in implementing any correction or advice in relation to the DL_MESO Software that STFC has made available on the DL_MESO Website or by the Licensee s failure to acquire a licence of and to implement any new release or version of the DL_MESO So
97. ate_module This module requires preloading of the constants variables domain_module surface_module comms_module error_module numeric_container and field_module modules Like the field_module different versions of the subroutine mdvv dpdvv lowe peters stoyanov are available for each thermostat integrator verlet_ e Header records SUBROUTINE verlet_ stage e Function Solves the equations of motion using the Velocity Verlet scheme 54 e Dependencies hardreflect deportdata error exportvelocitydata DPD VV only dragforces_dpdvv DPD VV only importforcedata DPD VV only global_sca_sca_int global_sca_sca_dble global_sum_dble Lowe Peters Stoyanov global_sum_int Lowe Peters Stoyanov e Arguments stage input integer e Comments All versions of this subroutine use the standard Velocity Verlet algorithm to integrate the equations di F dt mi dri dt Both the MD VV and DPD VV algorithms use dissipative and random forces as the thermostat with the DPD VV algorithm repeating the calculation of dissipative forces at the end of the time step The Lowe Andersen Peters and Stoyanov Groot thermostats are applied after all other forces are integrated verlet_ _lang e Header records SUBROUTINE verlet_ _lang stage e Function Solves the equations of motion using the Velocity Verlet scheme coupled with a Langevin barostat 22 e Dependencies hardreflect 170 CHAPTER 11 DL MESO DPD PACKAGE
98. brigt double system cooling rate lbsysdt double top boundary cooling rate lbtopdt double bottom boundary cooling rate lbbotdt double front boundary cooling rate lbfrodt double 38 CHAPTER 4 DL_MESO_LBE BASIC DEFINITION Table 4 7 DL_MESO_LBE Parameters continued function parameter data type notation back boundary cooling rate lbbakdt double left boundary cooling rate lblefdt double right boundary cooling rate lbrigdt double fluid inverse relaxation time 7 t lbtf double A lbsy nf bulk fluid inverse relaxation time T7 hulk 1btfbulk double A lbsy nf solute inverse relaxation time 7 lbtc double A lbsy nc temperature inverse relaxation time 7 1btt double A lbsy nt particle interaction function gab lbg double A lbsy nf lbsy nf body force lbbdforce double A 3 lbsy nf interaction force lbinterforce double A 3x 1bsy nf 1bdm touter distribution function 1bf double A lbsitelength lbdm touter temporary function lbft double A 1bdm touter lbsy nf l1bsy nc equilibrium distribution lbfeq double A lbsy nq speed vector for the model lbv int A 3x 1bsy nq index for opposing speed lbopv int A lbsy nq weight factor of speed vector lbw double A lbsy nq MRT transformation matrix lbtr double A lbsy ngq lbsy nq MRT inverse transformation matrix lbtrinv double A lbsy nq lbsy nq MRT tuneable collision parameters lbmrts double A3 MRT tuneable equilibrium parameters lbmrtw double A3 number of pa
99. bsteer int lor 0 system molar weight lbsmw double noise intensity lbnoise double system molar volume lbsmv double initial system velocity lbiniv double A3 top boundary velocity lbtopv double A3 bottom boundary velocity lbbotv double A3 front boundary velocity lbfrov double A3 back boundary velocity lbbakv double A3 left boundary velocity lblefv double A3 right boundary velocity lbrigv double A3 constant incompressible fluid density po lbincp double A lbsy nf initial fluid density lbinip double A lbsy nf top boundary fluid density lbtopp double A lbsy nf bottom boundary fluid density lbbotp double A lbsy nf front boundary fluid density lbfrop double A lbsy nf back boundary fluid density lbbakp double A lbsy nf left boundary fluid density lblefp double A lbsy nf right boundary fluid density lbrigp double A lbsy nf initial composition lbinic double A lbsy nc top boundary composition lbtopc double A lbsy nc bottom boundary composition lbbotc double A lbsy nc front boundary composition lbfroc double A lbsy nc back boundary composition lbbakc double A lbsy nc left boundary composition lblefc double A lbsy nc right boundary composition lbrigc double A lbsy nc initial temperature lbinit double top boundary temperature lbtopt double bottom boundary temperature lbbott double front boundary temperature lbfrot double back boundary temperature lbbakt double left boundary temperature lbleft double right boundary temperature l
100. buttons will produce the Lattice Boltzmann and Dissipative Particle Dynamics panels respectively which will guide you through setting up input files modifying and compiling the program code running the simulation and gathering the results files 5 6 CHAPTER 2 THE DL_MESO GUI DI MESO DL_MESO Welcome to DL MESO R S Qin W Smith amp M A Seaton Copyright STFC 2011 Click one of the above buttons to start Figure 2 1 DL MESO GUI on startup for plotting and visualization The SPH button is for Smoothed Particle Hydrodynamic simulations which will be included in future versions of DL MESO clicking on this button will currently produce a warning message This user manual can be read in Adobe Acrobat Reader if installed by clicking the Manual button while Help will advise you to visit the DL_MESO website at www ccp5 ac uk DL_MESO 2 2 Lattice Boltzmann and the DL_MESO GUI To access the LBE facilities in the DL_MESO GUI proceed as follows e Click the LBE button to get the LBE panel e Click the Define LBE System button and supply the required information The file 1bin sys will be created by the step e Click the Set LBE Space button to define the simulation space The file 1bin spa will be created by this step 2 2 1 Defining the System Figure 2 2 shows the Define LBE System panel with the rows for data entry labelled in red numbering The required data are as follows 1 The requi
101. c Numbers of bonded particles come from the MOLECULE file and therefore do not need to be included in CONTROL 4Parameters after the interaction length B are not essential for basic DPD calculations although C is additionally required for charged systems all omitted parameters will be set to zero with the exception of drag coefficients or frequencies which will be set to the system wide value 11 3 DL_MESO_DPD DATA FILES 173 Table 11 1 Keywords for CONTROL file with single parameters keyword meaning valid values default value VOL initial volume of system real number gt 0 TEMP specified system temperature real number gt 0 PRESS specified system pressure real number gt 0 TSTEP duration of calculation timestep real number gt 0 RCUT interaction cutoff radius real number gt 0 RMBCUT many body interaction cutoff radius real number gt 0 rcut RHALO size of boundary halo real number gt 0 Largest required RELEC short range electrostatic cutoff radius real number gt 0 recut SRFZCUT surface interaction cutoff distance real number gt 0 reut GAMMA system wide drag coefficient collision frequency real number gt 0 TIMJOB maximum time available to run calculation real number gt 0 TCLOSE time required to shut down calculation real number gt 0 KRES calculation restart parameter integer 0 3 KTYPE potential energy function selection integer 0 3 NSYST total number of unbonded particles inte
102. ch B Hasslacher and Y Pomeau Lattice gas automata for the Navier Stokes equation Physical Review Letters 56 14 1505 1508 April 1986 J B Gibson K Chen and S Chynoweth The equilibrium of a velocity Verlet type algorithm for DPD with finite time steps International Journal of Modern Physics C 10 1 241 261 February 1999 Minerva Gonz lez Melchor Estela Mayoral Maria Eugenia Vel zquez and Jos Alejandre Electro static interactions in dissipative particle dynamics using the Ewald sums Journal of Chemical Physics 125 22 224107 2006 R D Groot Electrostatic interactions in dissipative particle dynamics simulation of polyelectrolytes and anionic surfactants Journal of Chemical Physics 118 24 11265 11277 2003 Robert D Groot and Patrick B Warren Dissipative particle dynamics Bridging the gap between atomistic and mesoscopic simulation Journal of Chemical Physics 107 11 4423 4435 1997 Zhaoli Guo Chuguang Zheng and Baochang Shi Discrete lattice effects on the forcing term in the lattice Boltzmann method Physical Review E Statistical Nonlinear and Soft Matter Physics 65 4 046308 April 2002 Cecil Hastings Approximations for digital computers Princeton University Press Princeton NJ 1955 201 202 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 BIBLIOGRAPHY Xiaoyi He Shiyi Chen and Gary D Doolen A novel
103. ch gives a negative value for ij i if the vector Fij X Fix Tjk X Tki is in the same direction as the bond vector Fi and positive if in the opposite direction The force on particle acting in the a direction is given by O Fe U i 10 79 ara 9 jkl U hiu B Piss Fino u 10 80 i Tij Vik T x sin Qijkl O ijki ees ore Hea Using the following definition a 5 F 0 sap 046 o p the derivative of B fij Tjk Tkl is given by O 1 ps 3 e ary B Pigs Paks Tet Fij X Fi Fin X Tel Org Fig X Tir Fie X Fha COS Pi jk 1 o 1 o ror 10 81 2 e x Fel Ore Fi x Foal PT x Tal Ove g lje X Fial ae with Ore Fig X The Fie X Th v FpaP ely Sek 6a FyaPe Sex dej Tik Fij Piel Ser Sex Frei 005 Sei rir FigF ela Sew 005 Finn ei 5ez Arve Pig Thl g Ser Sex 10 82 a ae ed are malig X Fel 217 FieP jel Sey See FigF je 005 dex ork Fij Tis Sex ej Pisin Sei 5e 10 83 a ee ae are a Pik X Pel wey FinP jel Oe Sen FinTral ej der 215 Furl Sek 005 FinFei Sen Ser 10 84 It can be shown both algebraically virial 47 and thermodynamically that the dihedral makes no contribution to the Improper dihedrals which limit the geometry of molecules can be applied using the same procedure as standard dihedrals and no distinction is m
104. ck the names of any self defined variables whenever the package is updated to reduce the possibility of duplications causing unexpected errors The notation column in Table 9 4 gives the restrictions applicable on the parameters A indicates an array of data followed by the number of elements in the array For example A maxdim means the parameter is actually an array with maxdim elements numbered from 1 to maxdim gt 1 means the number must be greater or equal to one while for a Boolean parameter T or F means its value can either be true or false An asterisk in the data type for the array indicates that it is allocatable while two asterisks denote an allocatable array that is deallocated before DPD calculations start 9 2 THE PARAMETERS AND THEIR FUNCTIONS Table 9 4 DL_MESO_DPD Parameters 119 function parameter data type notation kind parameter for double precision numbers dp integer kind parameter for long integers di integer maximum number of particles maxdim integer maximum message buffer size maxbuf integer maximum data stack size mxstak integer maximum number of interaction parameters mxprm integer maximum number of species mxspe integer maximum number of interaction pairings mxpot integer maximum number of molecule definitions mxmoldef integer maximum particle size for molecule mxmolsize integer maximum number of bonds per molecule mxbonds integer restart f
105. cle distribution function at a fixed speed boundary for compressible fluids e Dependencies Many for different lattice schemes e Arguments tpos input long integer prop input long integer uwall output array of doubles e Comments Planar surface calculations are based on 60 concave edges and corners use equilibrium boundary con ditions with the density on the edge and at the grid point assumed to be equal to the values at their nearest neighbours in the bulk fluid The array uwal1 is the velocity at the grid point for all fluids which is subsequently used for solute concentration and temperature boundary conditions fFixedSpeedFluidIncom e Header records int fFixedSpeedFluidIncom long tpos int prop double uwall e Function Calculates the particle distribution function at a fixed speed boundary for incompressible fluids e Dependencies Many for different lattice schemes e Arguments tpos input long integer prop input long integer uwall output array of doubles e Comments Planar surface calculations are based on 60 concave edges and corners use equilibrium boundary con ditions with the density on the edge and at the grid point assumed to be equal to the values at their nearest neighbours in the bulk fluid The array uwal1 is the velocity at the grid point for all fluids which is subsequently used for solute concentration and temperature boundary conditions fFixedDensityFluid Header records int fFixedDensityFluid long tpos int prop doub
106. cles DL MESO DPD divides up the particles and total system volume volm between the processing units available At any given time each process holds nbeads particles Each process also has nbonds bonds nangles bond angles and ndiheds bond dihedrals to deal with If bonds are dealt with locally only the bonds associated with the subdomain are calculated by each process otherwise all processes hold information on all bonds The Cartesian coordinates velocities and forces for the particles are each held in sets of three double precision arrays for x y and z components Particle positions relative to the volume modelled by the individual processor thus not absolute positions unless the serial version of DL_MESO_DPD is used are held in arrays xxx i yyy i and zzz i for particle i Particle velocities are held in vxx i vyy i and vzz i Two sets of arrays for the net forces acting on the particles are available fxx i fyy i and fzz i for forces that remain constant over each time step and fvx i fvy i and fvz i for forces that may vary during the time step e g drag forces for DPD Velocity Verlet integration The potential energy and virial for individual particles 9 1 DATA STRUCTURE 117 Table 9 2 Domain information parameter meaning idnode name of the processor nodes number of processors idx x coordinate of the processor idy y coordinate of the processor idz z coordinate of the processor npx number
107. compi_ vtk with the composition of species in the cells starting from species number 1 e molcompi_ vtk with the composition of molecule type i in the cells starting with number 0 for all unbonded particles The scalar properties including compositions may be considered to act across the entire volumes of the cells while the velocities are representative for the cell centres Appendix C DL MESO Licence Agreement Academic Purposes 1 DEFINITIONS AND INTERPRETATION 1 1 In this Licence Agreement the following expressions have the meanings set opposite Academic Purposes a Derived Work the DL_MESO Software the DL_MESO Website a Harmful Element fundamental or basic research or academic teaching including any fundamental research that is funded by any public or charitable body but not any purpose that generates revenue as opposed to grant income for the Licensee or any third party Any research that is wholly or partially sponsored by any profit making organisation or that is carried out for the benefit of any profit making organisation is not an Academic Purpose any modification of or enhancement or improvement to any of the DL_MESO Software and any software or other work developed or derived from any of the DL_MESO Software the release and version of the DL_MESO Software downloaded by the Licensee from the DL_MESO Website immediately after the Licensee agrees to the terms and conditions of th
108. ct other and type the command for the required compiler in the neighbouring box The Create Makefile button needs to be clicked first to create a makefile in the working directory which automates compilation and may be edited by the user Clicking the COMPILE button will invoke the makefile to compile the code and a message box will signal its completion If the compilation fails you may need to edit the code An editing panel is available for this purpose using either the Change LBE Code or Change DPD code buttons Its function is similar to the compilation panel in operation with a choice of text editors including one packaged with the DL_MESO GUL The files in the LBE code that can be edited include the parallel and serial main files the lattice model file the boundary condition file the core routines for LBE calculations the file for user defined routines the main head file and the head file for user defined routines Others can be edited by selecting other and typing the name of the file in the neighbouring box The files in the DPD code that may be edited include the main program constants global variables and the modules configuration start for system initialization field for force calculations bond interactions many body DPD surfaces and statistics Other code files can be edited by selecting other and typing the name in the neighbouring box Running the LBE DPD code is made possible through the Run LBE Program or
109. cules Journal of Chemical Physics 116 13 5842 5849 2002 M Yoshino and T Inamuro Lattice Boltzmann simulations for flow and heat mass transfer problems in a three dimensional porous structure International Journal for Numerical Methods in Fluids 43 2 183 198 2003 Qisu Zou and Xiaoyi He On pressure and velocity boundary conditions for the lattice Boltzmann BGK model Physics of Fluids 9 6 1591 1598 June 1997
110. d J M Yeomans Lattice Boltzmann simulation of nonideal fluids Physical Review Letters 75 5 830 833 July 1995 I T Todorov and W Smith The DL POLY_4 user manual STFC STFC Daresbury Laboratory Dares bury Warrington Cheshire WA4 4AD United Kingdom version 4 01 0 edition October 2010 S Y Trofimov E L F Nies and M A J Michels Thermodynamic consistency in dissipative parti cle dynamics simulations of strongly nonideal liquids and liquid mixtures Journal of Chemical Physics 117 20 9383 9394 2002 204 54 55 56 57 58 59 60 BIBLIOGRAPHY Loup Verlet Computer experiments on classical fluids I Thermodynamical properties of Lennard Jones molecules Physical Review 159 1 98 103 July 1967 P B Warren Vapor liquid coexistence in many body dissipative particle dynamics Physical Review E 68 6 066702 December 2003 P B Warren P Prinsen and M A J Michels The physics of surfactant dissolution Philosophical Transactions of the Royal Society of London Series A Mathematical Physical and Engineering Sciences 361 1805 665 676 2003 John D Weeks David Chandler and Hans C Andersen Role of repulsive forces in determining the equilibrium structure of simple liquids Journal of Chemical Physics 54 12 5237 5247 1971 Satoru Yamamoto Yutaka Maruyama and Shi aki Hyodo Dissipative particle dynamics study of sponta neous vesicle formation of amphiphilic mole
111. d points using the swap algorithm e Dependencies fSwapPair e Comments This routine is fundamental to Lattice Boltzmann calculations and should not be modified Propagation is carried out by systematic swapping of post collisional values for the distribution function initially at each grid point and then between them in two separate loops as described by 34 and in section 5 1 This version can be used for either serial or parallel calculations and no boundary layer is necessary this is the default propagation routine for serial calculations fPropagationCombSwap e Header records int fPropagationCombSwap e Function Moves lattice particles distribution functions to neighbouring grid points using the swap algorithm e Dependencies fSwapPair e Comments This routine is fundamental to Lattice Boltzmann calculations and should not be modified Propagation is carried out by systematic swapping of post collisional values for the distribution function initially at each grid point and then between them in the same loop as described by 34 and in section 5 1 This version can only be used for calculations with non zero boundary layers this is the default propagation routine for parallel calculations 6 2 11 IbpMPI This package is only required for parallel running and does not require detailed knowledge for its use Sev eral subroutines in this package are not described here interested users should consult the code
112. d temperatures from data file lbin init and calculates initial distribution functions for incompressible fluids Dependencies lbin init data file fReadInitialState2DIncom fReadInitialState3DIncom Arguments filename input array of characters Comments The default file name is 1bin init this can be changed by the user if the initial state data file has a different name This routine will replace the default values of all properties at specified points and should be called after fInitializeSystemIncon fSetoffSteer e Header records int fSetoffSteer e Function Creates a file called notsteer to prevent DL_MESO_LBE from creating new 1bin sys and 1bin spa files which occurs if notsteer is missing e Dependencies None e Comments If the user has changed the input datafiles the code will run with new parameters fCheckSteer e Header records int fCheckSteer 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 67 e Function Checks for the existence of notsteer files if found then reads lbin sys and lbin spa files e Dependencies fInputParameters lbin sys fReadSpaceParameter lbin spa 6 2 6 IbpIOPlot3D fOutputGrid e Header records int fOutputGrid3D const char filename 1bout e Function Outputs grid positions for system in Plot3D format e Dependencies fOutputGrid2D fOutputGrid3D e Arguments filename input array of characters e Comments The default output file name is lbou
113. ded five lines of information of which two are mandatory have to be included e The simulation name 80 characters e The CONFIG file key levcfg integer the periodic boundary key imcon integer the number of particles in the file integer optional and the configuration energy real optional e The x y and z components for the x axis vector real optional This file is formatted identically to CONFIG files used in DL_POLY 48 52 except that the origin is set as the back bottom left corner of the simulation volume instead of the centre 176 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE e The x y and z components for the y axis vector real optional e The x y and z components for the z axis vector real optional The file key levcfg is set depending on the information available for each particle 0 for positions only 1 for positions and velocities or 2 for positions velocities and forces If particle velocities are not specified these are generated at random to produce a distribution corresponding to the required system temperature while unknown forces are set to zero The simulation name number of particles in the file and configuration energy are not read by DL MESO_DPD and can thus be ignored although the line for the simulation name must remain The axes vectors are not used by DL_MESO_DPD and the lines specifying them are skipped over if the value of imcon is set to a value greater than zero if set to ze
114. der records int D3Q15 e Function Assign the weight factor lbw speed vector 1bv index for opposite speed vector lbopv MRT transformation matrices 1btr and 1btrinv and tuneable parameters 1bmrts and 1bmrtw for the D3Q15 lattice model e Dependencies None D3Q19 e Header records int D3Q19 e Function Assign the weight factor lbw speed vector 1bv index for opposite speed vector lbopv MRT transformation matrices 1btr and 1btrinv and tuneable parameters 1bmrts and 1bmrtw for the D3Q19 lattice model e Dependencies None D3Q27 e Header records int D3Q27 e Function Assign the weight factor lbw speed vector 1bv and the index for opposite speed vector lbopv for the D3Q27 lattice model e Dependencies None e Comments No Multiple Relaxation Time MRT scheme is currently available for this lattice 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 53 6 2 3 IbpBASIC This package contains general purpose functions which are not directly related to the Lattice Boltzmann Equa tion method These may be replaced with any suitable functions in C standard libraries fCppAbs e Header records template lt class T gt T fCppAbs T a return a lt 0 a a e Function Calculate absolute value of numerical number a e Dependencies None e Arguments a input any datatype fCppAbs output same datatype as a fSwapPair e Header records template lt class T gt void fSwapPair T amp a T amp b e Functi
115. e Particle Dynamics included with version 2 0 and later 1 2 Functionality The following is a list of the features that DL_MESO currently supports Users are reminded that we are interested in hearing what other features could be usefully incorporated We obviously have ideas of our own and CCP5 strongly influences developments but other input would be welcome nevertheless 1 2 1 Lattice Boltzmann Equation DL_MESO_LBE can simulate lattice gas systems using the Lattice Boltzmann Equation LBE The following properties and features are currently available e Multiple fluid components solutes and coupled heat transfers 59 e Collisions Bhatnagar Gross Krook BGK single relaxation time 4 or Multiple Relaxation Time MRT 26 7 e Boundary conditions Periodic bounce back including stationary objects constant pressure velocity at planar surfaces 60 2 CHAPTER 1 DL MESO GENERAL INFORMATION e Mesoscale interactions Shan Chen pseudopotential method 43 44 1 2 2 Dissipative Particle Dynamics DL_MESO_DPD can model DPD particles beads with soft potential fields by default the standard Groot Warren potential 13 although Lennard Jones 23 and Weeks Chandler Andersen 57 potentials are also avail able along with thermostating dissipative and random forces The following properties and features are currently available e Choice of integrators thermostats standard Velocity Verlet DPD Velocity Verlet 10 Lowe And
116. e force is equal to m Tij rig Tay FS Ay 1 By or pj 1 2 10 46 g ay 1 2 85 0 1 2 2 10 46 132 CHAPTER 10 DL MESO_DPD FEATURES In the routine provided the terms with A for both force and potential are calculated as though they are standard DPD By setting A lt 0 and Bij gt 0 a vapour liquid mixture can be modelled and its equation of state for a single component is given as p pkgT aAp 2aBr4 p cp d where a 0 101 c and d are numerical offsets Parameters A and B for each species and interacting pair of species can be specified in the CONTROL file using the keyword MANY 10 5 Long ranged Electrostatic Coulombic Potentials Compared to other interactions in DPD electrostatic interactions act over considerably longer ranges which can also include periodic images of the system The governing equation for finding the electric potential is the Poisson equation shown here in dimensionless form 12 V PVY Tp 10 47 where y is the electric potential p the charge density concentration of cations minus concentration of anions per unit volume p 7 the local polarizability relative to a reference medium e g water and I the coupling constant for the reference medium The latter is given by e2 kel E0Erre with e as the electron charge 9 the dielectric constant of a vacuum and e the relative permittivity of the 1 reference medium For water at room tempe
117. e input array of characters e Comments The default output file name is lbout q This can be changed by specifying an output file name when calling the routine 6 2 7 IbplOLegacyVTK fOutputLegacy VTK e Header records int fOutputLegacyVTK const char filename lbout int iprop 0 e Function Outputs macroscopic mass density and velocity speeds along x y and z directions at each lattice point of compressible fluid iprop for system in Structured Grid Legacy VTK format e Dependencies fOutputLegacyVTK2D fOutputLegacyVTK3D e Arguments filename input array of characters e Comments The default output file name is 1bout q and the density of fluid 0 is output by default This can be changed by specifying an output file name and fluid number up to 1bsy nf 1 when calling the routine fOutputLegacyVTKIncom e Header records int fOutputLegacyVTKIncom const char filename lbout int iprop 0 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 71 e Function Outputs macroscopic mass density and velocity speeds along z y and z directions at each lattice point of incompressible fluid iprop for system in Structured Grid Legacy VTK format e Dependencies fOutputLegacyVTK2DIncom fOutputLegacyVTK3DIncom e Arguments filename input array of characters e Comments The default output file name is lbout vtk and the density of fluid 0 is output by default This can be changed by specifying an output file name and fluid number up
118. e modified calculation times for this thermostat may be comparable to those for the Lowe Andersen thermostat when At x 1 10 2 5 Stoyanov Groot thermostat itype 4 The Stoyanov Groot thermostat 49 is a combination of the Lowe Andersen thermostat and a Galilean invariant Nos Hoover thermostat which acts locally and on pairs of particles During force calculations after the first Velocity Verlet stage the choice to use either the Lowe Andersen or Nos Hoover thermostats for each particle pair is made at random the Lowe Andersen thermostat is selected with a probability of TAt The system temperature is also determined in terms of relative velocities for all particle pairs i e Dang WT rag hitii 3 pare YT rij where 1 7 rij is a smearing function for the temperature chosen to reduce to zero when r gt Te by default this kpT 10 13 is set as wT rij 1 for r lt re For all particle pairs that are to be subjected to the Nos Hoover thermostat an additional thermostating force is included gt BIEN ce ae E aw rig 1 p Oy gt 55 i 10 14 with a as a system wide thermostat coupling parameter and w rij as a switching function which by default is equivalent to w r 1 for standard DPD A virial correction of FI particle pair All other particle pairs are thermostated using the standard Lowe Andersen method Tij is also made for each This thermostat can produce a wide range of
119. e original author of the Lattice Boltzmann Equa tion source code DL_MESO _LBE and the DL_MESO graphical user interface and Richard Anderson at STFC Daresbury Laboratory for his contributions to the Dissipative Particle Dynamics source code DL_MESO_DPD Chapter 1 DL_MESO General Information 1 1 Description DL_MESO is a general purpose mesoscopic simulation package developed at Daresbury Laboratory by Dr Michael Seaton under the auspices of the Engineering and Physical Sciences Research Council EPSRC for the EPSRC s Collaborative Computational Project for the Computer Simulation of Condensed Phases CCP5 The package is the property of the Science and Technology Facilities Council STFC DL_MESO is issued free under licence to academic institutions pursuing scientific research of a non commercial nature All recipients of the code must first agree to the terms and conditions of the licence and register with us to be kept aware of new developments and discovered bugs Commercial organisations interested in acquiring the package should approach the Computational Science and Engineering Department STFC Daresbury Laboratory in the first instance Daresbury Laboratory is the sole centre for distribution of the package Under no account is it to be redistributed to third parties without consent of the owners DL_MESO contains two mesoscale simulation methods e Lattice Boltzmann Equation included with version 1 0 and later e Dissipativ
120. e the thermostat integrator type and its parameters see Table 10 1 THERM itype abcd 5 5Setting any of these parameters to 0 0 will make DL MESO DPD use the default values i e Aij Aj A Bij 3 Bi B and qij or Ij takes the system wide value 174 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE e the barostat type and its parameters see Table 10 2 BARO btype a bcd e the electrostatics method and its parameters see Table 10 3 ELEC etype a bcd e maximum number of periodic images in each dimension for Ewald sums KMAX ke ket ky e non periodic surface type applicable boundaries and where applicable surface parameters SURFACE srftype srfx srfy srfz etc e constant external body forces per particle for every species BODY F Fy F Note that not all parameters have to be included in the CONTROL file if they are not required for the given type of bond thermostat barostat electrostatics surface etc Superfluous parameters for particular systems e g specified pressure for constant volume simulations can be omitted from the CONTROL file without causing runtime problems If the user wishes to include new parameters and or keywords in the CONTROL file modifications to the parameter recognition loop in the datsys subroutine config_module will be required Define molecules MOLECULE and FIELD A single MOLECULE file and a single FIELD file are required for including molecules in DPD calculations The MOLEC
121. eader records int fCollisionMRTIncom e Function Collision steps for all incompressible fluids using MRT model e Dependencies fSiteCollisionMRTIncom e Comments This routine is fundamental to Lattice Boltzmann calculations and should not be modified fCollisionMRTGuo e Header records int fCollisionMRTGuo O e Function Collision steps for all compressible fluids using MRT model with Guo like forcing terms e Dependencies fSiteCollisionMRTGuo e Comments This routine is fundamental to Lattice Boltzmann calculations and should not be modified fCollisionMRTGuolncom e Header records int CollisionMRTGuoIncom e Function Collision steps for all incompressible fluids using MRT model with Guo like forcing terms e Dependencies f SiteCollisionMRTGuoIncom e Comments This routine is fundamental to Lattice Boltzmann calculations and should not be modified fPropagation e Header records int fPropagation e Function Moves lattice particles distribution functions to neighbouring grid points using the two lattice algorithm e Dependencies None e Comments This routine is fundamental to Lattice Boltzmann calculations and should not be modified This is the least efficient propagation routine available 98 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE fPropagationSwap e Header records int fPropagationSwap e Function Moves lattice particles distribution functions to neighbouring gri
122. ed as VCBPSD and translated as 322 34 CHAPTER 4 DL_MESO_LBE BASIC DEFINITION Table 4 3 Notation of boundary condition VCTPST VCTPSR VCTCCTRB VCTCCDRB VCTCCDLF VCTCETL VCTCETF VCTCERF VCTCEDB VBBPSD VBBPSF VBBCCTLB VBBCCTRE VBBCCDRF VBBCEDL VBBCELF VBBCETB VBBCERB VCBPSL VCBPSB VCBCCDLB VCBCCTLF VCBCETR VCBCEDR VCBCEDF VCBCELB VBTPST VBTPSR VBTCCTRB VBTCCDRB VBTCCDLF VBTCETL VBTCETF VBTCERF VBTCEDB PCTPSD PCTPSF PCTCCTLB PCTCCTRF PCTCCDRF PCTCEDL PCTCELF PCTCETB PCTCERB PBTPSL PBTPSB PBTCCDLB PBTCCTLF 121 124 127 130 133 144 147 150 153 222 225 228 231 234 245 248 251 254 323 326 329 332 343 346 349 352 421 424 427 430 433 444 447 450 453 522 525 528 531 534 545 548 551 554 623 626 629 632 VCTPSD VCTPSF VCTCCTLB VCTCCTRF VCTCCDRF VCTCEDL VCTCELF VCTCETB VCTCERB VBBPSL VBBPSB VBBCCDLB VBBCCTLF VBBCETR VBBCEDR VBBCEDF VBBCELB VCBPST VCBPSR VCBCCTRB VCBCCDRB VCBCCDLF VCBCETL VCBCETF VCBCERF VCBCEDB VBTPSD VBTPSF VBTCCTLB VBTCCTRE VBTCCDRF VBTCEDL VBTCELF VBTCETB VBTCERB PCTPSL PCTPSB PCTCCDLB PCTCCTLF PCTCETR PCTCEDR PCTCEDF PCTCELB PBTPST PBTPSR PBTCCTRB PBTCCDRB PBTCCDLF 122 125 128 131 134 145 148 151 154 223 226 229 232 243 246 249 252 321 324 327 330 333 344 347 350 353 422 425 428 431 434 445 448 451 454 523 526 529 532 543 546 549 552 621 624 627 630 63
123. ed using the pull down list This pull down list will be greyed out for two dimensional systems The back boundary condition can be selected using the pull down list This pull down list will be greyed out for two dimensional systems Solid obstacles can be added to the calculation space by selecting the bounce back on grid or mid grid and obstacle types in the pull down lists entering its location on the grid and if necessary entering its size and clicking add obstacle e A single point will be located at x0 y0 z0 e A sphere is centred at x0 y0 z0 and has radius r e A two dimensional circle is centred at x0 yO and has radius r e A block has a vertex at x0 y0 z0 and has size dx dy dz both z0 and dz can be omitted for two dimensional blocks Note that lattice points well within an obstacle are set as blank sites i e they will be ignored in LBE calculations The entire system can be set up as a porous solid by selecting the bounce back type specifying a pore fraction and clicking set pore to randomly select an appropriate number of solid lattice sites Clicking the Create button will write all the lattice space data to a lbin spa file any lattice point defined more than once will hold its latest definition 2 3 Dissipative Particle Dynamics and the DL_MESO GUI To access the DPD facilities in the DL_MESO GUI proceed as follows Click the DPD button to get the DPD panel Click the De
124. eee 3 5 Mesoscale Interaction amic e ROR A ee RE e ete gid ea eee ee de A 3 6 Summary of Lattice Boltzmann Equation 2 0 2 0000 00 00 000000000 DL_MESO_LBE Basic Definition Axl dhatticesmodels lt n hes A an eee oe Mk Bs Bd Bee a hk dn ee Be SS we E AS Ads Data structure A A A Mote ck Wad Seer dd wie Sets A SE ee he ES 4 3 The Parameters and Their Functions DL_MESO_LBE Features 5 1 Collision and Propagation Algorithms aa Ea aa e a D A a E N 5 2 Boundary conditions e neah d a A ete aot aa ee 5 3 Mesoscale interactions airada a a Bae A A eed DA E AR A A 5 4 Diffusion and heat transfer 5 5 Compressible and incompressible fluids 2 2 0 20 e 11 17 18 19 21 21 21 22 23 24 24 25 25 32 36 6 II 10 11 12 DL_MESO_LBE Package Reference GU Overviews 0 Ab ds Baths o te i hele didas Et nI a ne E 6 2 DL_MESO_LBE Subroutines and Functions 6 3 DL_MESO_LBE Data Files DL_MESO LBE Examples dal 2D PTOSSULE O tt a E E Ah ene ae eR GS PA QI Shear eb Be ek ee i oie te ka Te BE tak ol heat man ie Gi Gates fo 2D Cylinder Blow ios dos m a ae Goh PA ae Re eb Be ele eke ay 2D Warman Vortex oi o oe EAS See ke oa AA DS a1 ele tee ck we ee ee ee ee ee bd 7 6 3D_PhaseSeparation 0200 eee ee CC SD Shear 21 ee tt oe ee hae he Oe en ee ts Dissipative Particle Dynamics DPD Dissipative Particle Dynamics Basic Theory
125. eriodic boundary condition 4 back periodic boundary condition ADD OBSTACLES bounce back type x0 dx a obstacle type yo dy z0 dz j add obstacle SET AS POROUS MEDIA bounce back type pore fraction 0 8 Y Figure 2 4 Set LBE Space 2 2 2 Defining the Space Properties If this option is selected before saving the LBE system data a warning message advising that the system should be re defined will appear Figure 2 4 shows the Set LBE Space panel with the rows for data entry labelled in red numbering The following data are required 1 The top boundary condition can be selected using the pull down list from e periodic e on grid bounce back e mid grid bounce back e fixed V velocity C concentration and T temperature e fixed V and C Neumann T e fixed V and T Neumann C e fixed V Neumann C and T e fixed P pressure or density C and T e fixed P and T Neumann C e fixed P and C Neumann T e fixed P Neumann C and T 2 The bottom boundary condition can be selected using the pull down list 3 The left boundary condition can be selected using the pull down list 4 The right boundary condition can be selected using the pull down list lFor a property DL_MESO currently only calculates V 0 by using on grid bounce back on the related distribution function 2 3 DISSIPATIVE PARTICLE DYNAMICS AND THE DL_MESO GUI 11 5 The front boundary condition can be select
126. ers in the integer array bndtb1 i j for bond i with j 1 representing the first particle in the pair j 2 for the second and j 3 giving the user defined bond type The location of the first particle in each bond pair determines the processing unit which holds this data thus movement of this reference particle across processes also causes the bond list entry to be transferred with it Bond angles and bond dihedrals are stored in similar tables angtb1 i j particles j from 1 to 3 angle type at j 4 and dhdtb1 i j particles j from 1 to 4 dihedral type at j 5 respectively using the second particle in each triple or quadruple as the reference particle Prior to force calculations a list of bonded particles in each process domain and the surrounding boundary halo is constructed 1lblclst i j to allow DL MESO_DPD to find the local number for a particle j 2 from its global number j 1 using a binary search This list may include duplicates for the same global particle number the local numbers giving the shortest distance between pairs of particles are selected and used 118 CHAPTER 9 DL_MESO_DPD BASIC DEFINITION 9 2 The Parameters and Their Functions Table 9 4 lists all the globally used parameters defined in DL_MESO_DPD as given in the constants variables and numeric_container modules Because DL_MESO is an ongoing project and new parameters might be added to the package in the future it is strongly recommended that users of DL_MESO che
127. ersen 27 Peters 38 and Stoyanov Groot 49 e Constant volume NVT or constant pressure NPT simulations with Berendsen 2 or Langevin 22 barostats e User selection of interaction lengths conservative and dissipative force parameters for each species and between unlike species Bond stretching angles and dihedrals between beads in user defined molecules Density dependent conservative forces many body DPD 37 53 Electrostatic potentials between charged beads using a modified Ewald summation 11 Periodic or hard reflecting boundaries with optional short range repulsions 56 1 3 Requirements 1 3 1 Software requirements e Standard C Compiler for LBE source code DL_MESO_LBE e Standard Fortran90 Compiler for DPD source code DL_MESO_DPD e Message Passing Interface if parallel execution required e JAVA 2 Version 1 4 or higher if GUI is to be used e GNU Make included in standard Unix Linux distributions can be installed for Windows 1 3 2 System requirements DL_MESO is designed to work in both serial and parallel running it can be run on standalone machines clusters and supercomputers The code has been tested on Solaris Windows XP IBM p690 HPCx PowerPC 450 Blue Gene P and Cray XT4 XT6 HECTOR machines 1 4 The DL_MESO Directory Structure The supplied version of DL_MESO is a gzipped tar file which unpacks as directory dl_meso_2 x where x is a generation number Beneath the top level of this directory
128. ersion 1 4 or later These may be obtained from the java sun com website To build the GUI proceed as follows e Enter the DL_MESO JAVA directory e Type javac java to compile the source code e Type jar cfm GUI jar manifest mf class to create the GUI jar executable JAR file e Move to your working directory e Launch the GUI A Unix Linux script called makegui that performs the build of the GUI can be found in the JAVA subdirectory Your working directory is the directory from which you wish to work when running DL_MESO Working there will keep any files you generate separate from the DL_MESO source files Note in the current version of DL_MESO the working directory should be at the same directory level as the JAVA direction i e within the DL_MESO top directory and contain the executables of any external utilities required to set up input files and gather or process output files from simulations An example of such a working directory called WORK is present under the DL_MESO top directory this includes a makefile to compile external utilties which can be invoked by the command make f Makefile utils In your working directory you can start the GUI with the command e java jar JAVA GUI jar You may consider saving this command in a script for simple execution An example script for Unix Linux called rungui is present in the WORK subdirectory Figure 2 1 shows the DL_MESO GUI when it is started Clicking the LBE and DPD
129. es additionally dependent on local densities The free energy of an inhomogeneous system with density p r can be defined as the following in both continuous and ensemble averaged discrete forms F a 10 37 E VAD 10 38 where Y p is the free energy per particle in a homogeneous system and p 7 is a function related to the density at and near the position of particle i 7 The latter can be approximated by a function dependent on the positions of particles close to particle 6 to allow the calculation of an instantaneous free energy F vt 10 4 PARTICLE PARTICLE INTERACTIONS 131 The effective conservative force on particle i can be obtained from the spatial derivative of the free energy although only the excess part equivalent to the potential energy U is required since the kinetic motion of particles automatically accounts for the ideal contribution FE OF ric gt Ow pj 10 39 ieee Or The force can also be expressed in terms of pairwise interactions taking a form similar to the standard DPD conservative force Equation 8 5 Re E a Gy 10 40 a Op Op Tij The local density approximation can be defined as a weighted average of instantaneous densities Ba arwr i DEAR Y wlr R DaF 7 A i Y wer j 10 41 JA with we r as the weight function vanishing beyond a cutoff rq which can be equal to re or smaller and normalized so that fY 4rr w r dr 1 The most fre
130. es from a previous run can be used for restarting a previous simulation or providing a configuration for a new simulation providing the number of processing units remains the same Prior to compliation modifications may be made in the constants f90 file to the following maximum numbers e particles per process maxdim e species mxspe molecule definitions mxmoldef bond interaction definitions mxbonddef particles per molecule mxmolsize bonds per molecule mxbonds e size of statistical stack arrays mxstak The maximum sizes of transfer buffers maxbuf and potentials mxpot are automatically set by the values of maxdim and mxspe respectively and do not need to be modified by the user The value of maxdim should additionally take into account the number of particles likely to be included in boundary halos If using alternative many body DPD interactions to the vapour liquid example provided the routines manybody_force and manybody_potential in manybody_module f90 should be modified by the user as necessary the routine local_density should not be altered by the user but if an alternative weight function for calculating local densities is required weight_rho can be changed Appendix B DL MESO Utilities DL_MESO includes a number of utility programs which are not directly needed for Lattice Boltzmann or DPD simulations but are useful both for producing files required as inputs for those calculations and to process out
131. ese codes are shown below these model systems using single relaxation time collisions with mesoscopic interactions and stable boundary conditions and produce outputs in VTK format with densi ties for fluid 0 An additional version of the serial program is available which uses a boundary layer of lattice points slbecombine cpp although the code is not shown here The user may wish to modify the listings to use different outputs and formats alternative collision routines boundary conditions that move through the system etc and should refer to Appendix A for more details Serial code slbe include slbe hpp int main int argc char argv bigend fBigEndian fDefineSystem fSetSerialDomain fStartDLMESO fMemoryAllocation fInputParameters fReadSpaceParameter fGetModel O fInitializeSystem fReadInitialState fsDutputGrid Only needed for Plot3D output files timetotal fCheckTimeSerial for int i 0 i lt lbtotstep i if i lbsave 0 fsOutputVTK Substitute with required file type and or output qVersion printf hd Xt i fPrintDomainMass fPrintDomainMomentum fInteractionForceShanChen fCollisionBGK fPostCollBoundary fPropagationSwap fPostPropBoundary timetotal fCheckTimeSerial O fFreeMemory O fFinishDLMESO return 0 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS Parallel code plbe include plbe hpp
132. essure Pp recut interaction cutoff radius re rmbcut many body interaction cutoff radius rq relec short range electrostatic interaction cutoff radius re sr zcut surface repulsion cutoff distance ze rhalo size of boundary halo nrun number of calculation timesteps nseql number of equilibration timesteps tstep duration of calculation timestep At timjob maximum time available to run calculation tclose time required to shut down calculation kres calculation restart parameter nsbpo interval for printing data to OUTPUT file nsbsv interval for saving data to HISTORY file s nstk size of statistical data stack ltemp switch for temperature scaling before equilibration nsbts interval for temperature scaling lcorr switch for saving statistical data to CORREL file iscorr interval for saving statistical data to CORREL file ktype potential energy function selection itype integrator thermostat selection btype barostat selection etype electrostatic algorithm selection srftype surface boundary selection lbond switch for modelling bonds between particles langle switch for modelling bond angles 1dihed switch for modelling bond dihedrals lgbnd switch for globally storing all bond data lisoprs switch for isotropic variation of system dimensions with pressure 9 1 2 Storage of particle and bond properties The total number of particles in a system is nsyst of which nusyst particles are loose i e not bonded to other parti
133. f fluid Space is represented by a regularly distributed grid Time flow is obtained by integrating over discrete time steps Discrete velocity vectors are defined to ensure that a particle moves from one grid point to another without falling between grid points The Lattice Boltzmann method inherits the discrete concept from its ancestor Lattice Gas Cellular Automata LGCA The main difference between the Lattice Boltzmann Equation LBE method and LGCA is that LBE describes the physical state of a grid by a set of single particle distribution functions instead of a particle number This allows LBE to simulate condensed matter fluids distinct from diluted fluids in which the mean free path of the component particles is much larger than the lattice distance The Lattice Boltzmann algorithm can be summarized by the following e Fluid properties are mapped onto a discrete lattice e The physical state at each lattice point is described by a set of particle distribution functions e Macroscopic fluid variables are defined by moments of the distribution functions e The system evolution towards equilibrium is through the relaxation of the distribution function to its equilibrium form 3 2 Basic Definitions Triangular and rectangular lattices are two of the most popular grid forms The former have O6 rotation isotropy and have been widely applied in two dimensional systems e g the D2Q7 and D2Q13 models The latter have only O4 rotation isotrop
134. fGetEquilibriumIncom double feq double v double rho double rho0 e Function Calculates equilibrium distribution function for incompressible fluid e Dependencies None e Arguments feq output double pointer v input double pointer rho input double precision rho0 input double precision e Comments Equilibrium distribution function calculated here is pisu fosa which is only suitable for square lattices Further details can be found in 17 fGetMomentEquilibriumF e Header records int fGetMomentEquilibriumF double meq double p double rho e Function Calculates equilibrium distribution function in moment space for compressible fluid 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 89 e Dependencies None e Arguments meq output double pointer p input double pointer rho input double precision e Comments The equilibrium distribution function in moment space calculated here is Mea T the exact form of which is dependent on the lattice scheme given for D2Q9 by 26 and for D3Q15 and D3Q19 by 7 Parameters for calculating the square of energy e and fourth order moments Taa can be modified by the user in the 1bpMODEL module fGet MomentEquilibriumIncom e Header records int fGetMomentEquilibriumIncom double meq double p double rho double rho0 e Function Calculates equilibrium distribution function in moment space for incompressible fluid e Dependencies None e Argument
135. fine DPD System button and supply the required information The CONTROL file will be created by this step Note that currently no simulation space settings or molecular structure data can be entered using the GUI Click EXIT to finish the settings 2 3 1 Defining the System Figure 2 5 shows the Define DPD System panel with the rows for data entry labelled in red numbering The required data are as follows 1 2 3 The job header a line of text up to 80 characters long describing the simulation The system volume which can either be cubic or orthogonal cuboidal as selected from the pull down list If specifying a cubic volume the total volume should be specified while orthogonal volumes require the sizes for all three dimensions The target temperature kgT and pressure P for the system The latter is greyed out if no barostat is to be used 12 10 11 12 13 14 15 CHAPTER 2 THE DL_MESO GUI The dissipative factor y for DPD and Peters thermostats T for Lowe Andersen and Stoyanov Groot thermostats for the entire system values for individual species and between unlike species can be specified elsewhere and the size of the boundary halo for copying particle data from neighbouring subdomains or across periodic boundaries The interaction cutoff re for pairwise particle interactions and the many body cutoff rq for determining localized particle densities as used by many body
136. fluid solute and thermal parameters can be seen in Figure 2 3 with the rows labelled in red numbering multiple columns of dialogue boxes are made available for systems with multiple fluids and or solutes MBE Aud properties body force x axis o 1 body force y axis 0 z body force z axis 3 fluid densities initial 1 4 constant 1 5 top 0 6 E o ie sore properties bottom 0 Fi left 0 8 solute concentration initial 0 1 right o 9 top o 2 front 10 bottom 0 3 back 11 left o 4 relaxation time 1 12 right o 5 bulk relaxation time 1 13 front 6 interaction parameters fluid O back F fluido o 14 solute relaxation time 1 8 SAVE F CANCELF savec CANCEL c a Fluid parameters b Solute parameters DBE thermal properties heat relaxation time 1 initial T K initial dT dt wall T K dT dt K s wall T K dT dt K s 3 top bottom 4 left right 5 front back SAVE T CANCEL T c Thermal parameters Figure 2 3 Fluid solute and thermal parameter pop up windows For fluids the required data are as follows 1 Body force x axis the x component of body force on each fluid 2 Body force y axis the y component of body force on each fluid 2 2 10 11 12 13 14 L
137. forces from edges of fluid points e Dependencies fsPeriodic2D fsPeriodic3D e Comments Serial equivalent of f ForceNonBlockCommunication may be required for using combined collision propagation routines Collision Propagation in serial running 6 2 10 IbpSUB fWeakMemory e Header records inline void fWeakMemory e Function Terminates calculation if system has insufficient memory e Dependencies None e Comments If called will print error message error cannot allocate more memory 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 85 fMemory Allocation e Header records int fMemoryAllocation e Function Allocates memory for lattice Boltzmann calculations e Dependencies fWeakMemory e Comments If memory allocation is unsuccessful will print error message and stop calculation fFreeMemory e Header records int fFreeMemory e Function Frees allocated memory e Dependencies None fSetSerialDomain e Header records int fSetSerialDomain e Function Sets domain parameters for serial running e Dependencies None e Comments Default routine sets domain boundary width 1bdm bwid to zero fSetSerialDomainBuffer e Header records int fSetSerialDomainBuffer e Function Sets domain parameters for serial running including additional boundary points e Dependencies None e Comments Similar to fSetSerialDomain but does not modify the domain boundary width from its user
138. forces_ nlimit e Function Calculates all forces between particles particularly pairwise forces within cut off radius 168 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE e Dependencies mtrnd duni gaussmp manybody_force manybody_potential ewald_real ewald_reciprocal bondforcesglobal bondforceslocal hardreflect e Arguments nlimit input integer e Comments The random force is calculated using a uniform random number generator 13 by default the Mersenne twister generator mtrnd i e Gij x V12 u 0 5 which produces statistically similar results 8 and is computationally more efficient than using the Gaus sian random number subroutine gaussmp although this or the duni random number generator may be substituted if required dragforces dpdvv e Header records SUBROUTINE dragforces_dpdvv nlimit e Function Recalculates dissipative forces between particles within cut off radius e Dependencies None e Arguments nlimit input integer idnode input integer e Comments Only required for DPD Velocity Verlet DPD VV integration plcfor_ e Header records SUBROUTINE plcfor_ e Function Sets up parallel link cells structure and calculates forces including those involving particles in boundary halo e Dependencies exportdata parlnk local_density importdensitydata exportdensitydata 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 169 forces_ importdata importdata_dpdvv 11 2 14 integr
139. ftware that would have remedied or mitigated the effects of any error defect bug or deficiency in the DL_MESO Software The Licensee acknowledges that proper use of the DL_MESO Software and any Derived Work is dependent on the Licensee its employees and students exercising proper skill and care in inputting data and interpreting the output provided by the DL_MESO Software or that Derived Work STFC will not be liable for the consequences of decisions taken by the Licensee or any other person on the basis of that output STFC does not accept any responsibility for any use which may be made by the Licensee of that output nor for any reliance which may be placed on that output nor for advice or information given in connection with that output Subject to clause 3 7 STFC s liability or any breach of this Licence Agreement any negligence or arising in any other way out of the subject matter of this Licence Agreement or the use of the DL_MESO Software will not extend to any incidental or consequential damages or losses or any loss of profits loss of revenue loss of data loss of contracts or opportunity whether direct or indirect even if the Licensee has advised STFC of the possibility of those losses arising or if they were or are within STFC s contemplation Subject to clause 3 7 the aggregate liability of STFC for any and all breaches of this Licence Agreement any negligence or arising in any other way out of the subject mat
140. g z axis lbdm zdim 1 when lbsy nd 2 lbdm xs x coordinate of domain start position lbdm xe xa coordinate of domain end position lbdm xinner number of grid points along z axis in the domain lbdm xouter number of grid points along x axis in the domain including the boundary lbdm ys y coordinate of domain start position lbdm ye y coordinate of domain end position lbdm yinner number of grid points along y axis in the domain lbdm youter number of grid points along y axis in the domain including the boundary lbdm zs z coordinate of domain start position lbdm ze z coordinate of domain end position lbdm zinner number of grid points along z axis in the domain lbdm zouter number of grid points along z axis in the domain including the boundary lbdm touter total number of grid points in the domain including the boundary Table 4 6 Neighbour information parameter meaning lbnb k rank processor name of neighbour k lbnb k spos start position for sending distribution function message lbnb k rpos start position for receiving distribution function message lbnb k bspos start position for sending boundary condition message lbnb k brpos start position for receiving boundary condition message k 0 k 1 k 2 k 3 k 4 k 5 right neighbour left neighbour upper neighbour lower neighbour front neighbour back neighbour 4 3 The Parameters and Their Functions Table 4 7 lists all the parameters defined in DL_MESO_LBE The w
141. ger gt 0 NRUN number of calculation timesteps integer gt 0 NSEQL number of equilibration timesteps integer gt 0 NSBTS interval for equilibration temperature scaling integer gt 0 NSTK size of statistical data stack integer gt 0 NSBPO interval for printing data to OUTPUT file integer gt 0 NSBSV interval for saving trajectory data to HISTORY file s integer gt 0 ISCORR interval for saving statistical data to CORREL file integer gt 0 NSPE number of particle species integer gt 0 NMOLDEF number of molecule types integer gt 0 LCORR switch for saving statistical data to CORREL file Boolean false LSBSV switch for saving trajectory data to HISTORY file Boolean false LTEMP switch for temperature scaling during equilibration Boolean false LISOPRS switch for isotropic pressure tensors in barostat Boolean true LBOND switch for modelling bonds between particles Boolean false LANGLE switch for modelling bond angles Boolean false LDIHED switch for modelling bond dihedrals Boolean false LGBND switch for globally storing bond data Boolean false e name and number of individual molecules e unlike conservative potential parameters interaction lengths and drag coefficients collision frequencies MOLE name ni UNLIKE 2 j Ajj Bi Vig Vig e many body DPD parameters for each species and interacting species pair MANY 2 j Ajj Bi Cij Di Esj e the dimensions for non cubic orthorhombic systems DIM Le Ly Le
142. given for D2Q9 by 26 and for D3Q15 and D3Q19 by 7 Tuneable parameters for calculation stability can be modified by the user in the LbpMODEL module fInitializeSystem e Header records int fInitializeSystem e Function Initializes distribution function for lattice system using compressible fluids e Dependencies fGetEquilibriumF fGetEquilibriumC fGetEquilibriumT e Comments This subroutine as it stands is suitable for initializing most Lattice Boltzmann systems with compressible fluids although the user may modify it if it can otherwise be faster more stable or more suitable for a particular calculation fInitializeSystemIncom e Header records int fInitializeSystem e Function Initializes distribution function for lattice system using incompressible fluids 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 91 e Dependencies fGetEquilibriumIncom fGetEquilibriumC fGetEquilibriumT e Comments This subroutine as it stands is suitable for initializing most Lattice Boltzmann systems with incompressible fluids although the user may modify it if it can otherwise be faster more stable or more suitable for a particular calculation fSiteCollisionBGK e Header records int fSiteCollisionBGK doublex startpos double bodyforce e Function Calculates collisions at a grid point using the Bhatnagar Grook Krook BGK model for compressible fluids e Dependencies fGetOneMassSite fGetSpeedSite fGetEquilibriumF fGetEqui
143. gure 12 4 gives a comparison between the polyelectrolyte and neutral polymer at the final time step which have measured radii of gyration of 4 5 and 2 7 respectively 12 5 AMPHIPHILEMESOPHASES 181 a Polyelectrolyte b Neutral polymer Figure 12 4 Visualizations of DPD Polyelectrolyte test case red for polyelectrolyte polymer green for salt cations cyan for anions orange for counterions water omitted for clarity 12 5 AmphiphileMesophases This serial simulation example consists of four separate simulations each with 12 000 particles consisting of dimers molecules consisting of two particles one hydrophilic and the other hydrophobic with harmonic bonds of equilibrium length 0 5 between them and unbonded monomers 24 Defining the composition as the ratio of DPD particles within dimers to the total number of particles in the system simulations with dimer compositions of 30 55 70 and 90 are provided with filenames CONTROL 30 CONTROL 55 CONTROL 70 and CONTROL 90 respectively Each of these files should be renamed to CONTROL to run the simulation while the FIELD and MOLECULE files defining the amphiphile dimer molecule should be used for all four simulations These systems provide four points on a phase diagram corresponding to isotropic dimer hexagonal lamellar and isotropic monomer phases respectively The final configurations obtained for each phase can be seen in Figure 12 5 shown as isosurfaces of the hydrophobic
144. h multiple species For species the required data are as follows referring to Figure 2 6 a 1 Species name this can be up to 8 alphanumeric characters overly long names will be truncated 2 Species population the number of unbonded particles of the species in the system If molecules are to be included species that do not exist in unbonded form should still be defined here but the unbonded population should be set to zero 14 CHAPTER 2 El DPD potential properties unlike interactions 1 species i species j Aij Dij add unlike reorder unlike reset unlike Ed DPD species properties alla species 1 name ser 1 population unbonded lo 2 mass lo 3 parameter A energy 250 4 parameter B distance Lo 5 parameter C charge oo 6 parameter D dissipation oo 7 body force x axis oo 8 body force y axis loo 9 body force z axis lo o 10 SAVE SP CANCEL SP a Species parameters EJ DPD thermostat properties MAK STOYANOV GROOT 1 global coupling parameter 0 0 2 SAVE T CANCEL T c Thermostat parameters fj DPD electrostatic properties Soe EWALD SUM WITH SLATER SMEARING 1 system coupling constant Ewald real space convergence 0 0 0 0 2 savee cancere e Electrostatic parameters charge smearing 0 0 k vector x component fo 3 k vector y component o k vector z component 0 ma
145. hdtype input integer phi input real KIND dp pb input real KIND dp pc input real KIND dp a input real KIND dp b input real KIND dp c input real KIND dp d input real KIND dp force output real KIND dp potential output real KIND dp e Comments If required the user can add extra bond dihedral and improper dihedral interaction types in this sub routine as additional cases bondforceslocal e Header records SUBROUTINE bondforceslocal 156 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE e Function Calculates all bond stretching angle dihedral forces between particles in system using locally defined bond lists e Dependencies shellsort_list search_list duplicate error global_sca_and bond_force angle_force dihedral_force e Comments Assumes that bond lists contain only those bonds in the current node This is the most efficient method for parallel running but runs the risk of losing track of bonded pairs if bond lengths get longer than the subdomain halo size rhalo bondforcesglobal e Header records SUBROUTINE bondforcesglobal e Function Calculates all bond stretching angle dihedral forces between particles in system using globally defined bond lists e Dependencies shellsort_list search_list duplicate error global_sum_dble global_sca_and bond_force angle_force dihedral_force e Comments Assumes that bond lists contain all bonds in the entire system This is less efficient for paral
146. he CONTROL file using these checkboxes bond needs to be ticked if any molecules are to be included while angle and dihedral only need to be selected if those interactions are required for any of the molecules being modelled The global option makes DL_MESO_DPD store all bond data for the system in all processing units instead of by default storing bond data in each processing unit for only the bonds within its subdomain After the molecule name populations and types of interactions required are entered clicking SAVE will write this information to the existing CONTROL file 2 4 Compiling and running DL_MESO Compiling the LBE DPD code may be accomplished through the compiler panel which is activated from either of the Compile LBE Code or Compile DPD Code buttons The Compile LBE Code panel allows you to select the operating system a C compiler compiler flags and the version serial or parallel of the code you wish to build If you require a C compiler that is not included in the pull down list select other and type the command for the required compiler in the neighbouring box Clicking the COMPILE button will start the compilation and a message box will signal its completion The Compile DPD Code panel allows you to select the operating system a Fortran90 compiler compiler flags and the version serial or parallel of the code you wish to build If you require a Fortran90 compiler that is not included in the pull down list sele
147. hole range parameters are named with the prefix 1b Because DL_MESO is an ongoing project and new parameters might be added to the package in the future it is strongly suggested that users of DL_MESO would not name their own parameters with prefixes 1b dp or sp The notation column in Table 4 7 gives the restrictions applicable on the parameters A indicates an array of data followed by the number of elements in the array For example A lbsy nf means the parameter is actually an array with lbsy nf elements gt 1 means the number must be greater or equal to one while 1 or 0 means the value of the parameter can either be one or zero An asterisk in the data type for the array indicates that it is allocatable 4 3 THE PARAMETERS AND THEIR FUNCTIONS Table 4 7 DL_MESO_LBE Parameters 37 function parameter data type notation system information lbsy structure domain information 1bdm structure neighbour information lbnb structure A6 space dimension lbsy nd int number of discrete speeds lbsy nq int number of fluids lbsy nf int gt 1 number of solutes lbsy nc int temperature scalars lbsy nt int lor 0 grid points in x direction lbsy nx int grid points in y direction lbsy ny int grid points in z direction lbsy nz int domain boundary width lbsy bwid int 21 system dimension along x lbxsize double total calculation steps lbtotstep int data save interval lbsave int steering l
148. i Mats Aspn s and Jan Westerholm An efficient swap algorithm for the lattice Boltzmann method Computer Physics Communications 176 3 200 210 February 2007 BIBLIOGRAPHY 203 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 Philip M Morse Diatomic molecules according to the wave mechanics II Vibrational levels Physical Review 34 1 57 64 July 1929 David R Noble Shiyi Chen John G Georgiadis and Richard O Buckius A consistent hydrodynamic boundary condition for the lattice Boltzmann method Physics of Fluids 7 1 203 209 1995 I Pagonabarraga and D Frenkel Dissipative particle dynamics for interacting systems Journal of Chemical Physics 115 11 5015 5026 2001 E A J F Peters Elimination of time step effects in DPD EPL Europhysics Letters 66 3 311 317 May 2004 Kannan N Premnath and John Abraham Three dimensional multi relaxation time MRT lattice Boltzmann models for multiphase flow Journal of Computational Physics 224 2 539 559 2007 A G Schlijper P J Hoogerbrugge and C W Manke Computer simulation of dilute polymer solutions with the dissipative particle dynamics method Journal of Rheology 39 3 567 579 1995 Xiaowen Shan Analysis and reduction of the spurious current in a class of multiphase lattice boltzmann models Physical Review E 73 4 047701 2006 Xiaowen Shan Pressure tensor calculat
149. ian 0 for little endian 56 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE e Dependencies None fByteSwap e Header records void fByteSwap void data int len int count e Function Converts between endian types by swapping byte order of array data with byte size per entry len and count entries e Dependencies None e Arguments data input output void len input integer count input integer e Comments Primarily required for writing binary files where a specific endianness is required e g vtk files are required in big endian fCheckTimeSerial e Header records double fCheckTimeSerial e Function Outputs time in seconds since initial call e Dependencies None e Arguments fCheckTimeSerial output double e Comments Obtains calculation time based on system clock parallel calculations may obtain greater timing accuracy with fCheckTimeMPI 6 2 4 IbpGET fGetNodePosi e Header records two cases 3D and 2D inline long fGetNodePosi int xpos int ypos int zpos inline long fGetNodePosi int xpos int ypos e Function Calculates the position of the grid point in a one dimensional array from its Cartesian coordinate e Dependencies None 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 57 e Arguments xpos ypos zpos fGetNodePosi e Comments input integer input integer input integer output long integer The calculation follows the standard C data structure row major fGetCoord
150. ic boundaries do not require boundary halos these are excluded from the export of particle local densities 164 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE 11 2 11 start_module This module requires the modules constants variables error_module comms_module numeric_container and parse utils to be loaded beforehand start e Header records SUBROUTINE start e Function Sets up starting configuration for DPD calculations e Dependencies duni mtrnd revive boxfit readconfig initialvelocity error global_sca_and global_sum_dble e Comments If the user does not provide initial particle positions with optional velocities and force in a CONFIG file and the system is not restarted all unbonded particles are assigned their species randomly according to the numbers required this may be changed by the user revive e Header records SUBROUTINE revive key e Function Saves and reads restart configurations e Dependencies export files if reading restart configuration e Arguments key input integer readconfig e Header records SUBROUTINE readconfig e Function Reads initial system configuration positions velocities forces from DL_POLY style CONFIG file e Dependencies CONFIG file initialvelocity 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 165 global_sca_and global_sum_sca_int error e Comments The number of particles in CONFIG file must match up with the number described in the corresp
151. id grid bounce back Like the on grid version the mid grid bounce back condition applies a no slip condition at a boundary halfway between lattice points 50 This is applied by assigning post collisional distribution functions to the wall node based on those values at neighbouring points i e filEw t fy Bw 6 At t 5 12 This method essentially applies the actual reflection halfway between timesteps and is a spatially second order method although it is weakly non local due to its use of distribution functions from neighbouring nodes 5 2 4 Constant pressure velocity To specify either velocities or densities pressures at planar boundaries the Zou He method 60 is available in DL_MESO_LBE This is based upon applying the bounce back rule to the non equilibrium distribution functions Le ff Gay Ga 5 13 where fi fi fi with the equilibrium distribution function ff as a function of density and velocity This function can be used to determine the missing wall velocity or density along with the known distribution function values For instance for a top edge with a known velocity 7 using the D2Q9 lattice scheme V CEDF the wall density and missing distribution functions all for the boundary node at Zw are given as fot f2 fe 2 fi fr fa E 1 vw y fa fo Podes 4 8 3 fs fr H E fe f2 PV FPwUw y fs fi fo f2 H Uwe PwUw y The other form specifying the wall fluid density
152. ig At 2 Tij PR FR 8 6 in which At is the time step and w r is a switching function which imposes a finite limit on the range of the stochastic force is a random number with zero mean and unit variance The constant o is related to the temperature as is understood from the role of the stochastic force in representing a heat bath Finally the particles are subject to a drag force which depends on the relative velocity between interacting pairs of particles a a FR qijw riz Fig Uig T 8 7 ij where w rij is once again a switching function and v 0 0 The constant 7 is the drag coefficient It follows from the fluctuation dissipation theorem that for thermodynamic equilibrium to result from this method the following relations must hold of 2yijkBT 8 8 2 w rij w ri5 8 9 In practice the switching functions are defined through w rig 1 a rij lt re 8 10 8 3 EQUATION OF STATE AND DYNAMIC PROPERTIES 113 which ensures that all interactions are switched off at the range rij re In many DPD simulations the stochastic and drag coefficients are often constant for all interactions i e oij 0 and yi y although this assumption does not have to apply 8 3 Equation of state and dynamic properties The form of the conservative force determines the equation of state for a DPD fluid which can be derived using the virial theorem to express system pressure as fol
153. ile export creation interval ndump integer length of double precision MPI message Dlen integer name of processing unit idnode integer number of processing units nodes integer filename for restart files exportname character 10 switch for temperature scaling ltemp logical TorF switch for reading CONFIG file lconfig logical TorF switch for writing CORREL file lcorr logical TorF switch for writing HISTORY files lsbsv logical Tor F switch for modelling bonds lbond logical TorF switch for modelling bond angles langle logical TorF switch for modelling bond dihedrals ldihed logical TorF switch for global holding of bond information 1gbnd logical TorF switch for isotropic variations of volume with pressure lisoprs logical TorF printout selection for OUTPUT file outsel integer name of DL_MESO_DPD calculation text character 80 number of time steps for calculation nrun integer interval for writing OUTPUT file nsbpo integer gt 1 interval for writing CORREL file iscorr integer interval for writing HISTORY files nsbsv integer temperature scaling interval nsbts integer number of equilibration time steps nseql integer calculation restart parameter kres integer number of species nspe integer gt 1 lt mxspe species population of unbonded particles nspec integer A mxspe number of defined molecule types nmoldef integer lt mxmoldef size of statistical data stack nstk integer gt 1 lt mxstak total number of particles in system nsyst integer
154. ine argument is given the utility will ask the user to type in the number of processes and the CONFIG file key exportimage exportimage is a utility written in Fortran90 to produce a VTF format trajectory file export vtf from DL_MESO_DPD restart files export that can be visualized to give a snapshot of the last simulation timestep Since a limited amount of data is included in restart files the CONTROL file for the simulation is needed to provide some additional information The source code for this utility exportimage f90 can be used for export or export files created by both the serial and parallel versions of DL_MESO_DPD If the available Fortran90 compiler is invoked by the command 90 the executable export exe can be produced by typing e 90 o export exe exportimage f90 and either run at the command line or by using the GUI in Process DPD Data This utility can be run with a command line argument indicating the number of processes used to generate the restart data e g if 16 processes were used either of the following commands can be used 190 APPENDIX B DL_MESO UTILITIES e export exe 16 e export exe 16 If no command line argument is given the utility will ask the user to type in the number of processes The number of processes can be specified in the GUI before running the utility traject The Fortran90 utility traject reads in HISTORY output data files generated by DL_MESO_DPD and produces a VTF for
155. ing BGK model e Dependencies fSiteCollisionBGK e Comments This routine is fundamental to Lattice Boltzmann calculations and should not be modified 96 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE fCollisionBGKIncom e Header records int fCollisionBGKIncom e Function Collision steps for all incompressible fluids using BGK model e Dependencies fSiteCollisionBGKIncom e Comments This routine is fundamental to Lattice Boltzmann calculations and should not be modified fCollisionBGKGuo e Header records int fCollisionBGKGuo e Function Collision steps for all compressible fluids using BGK model with Guo forcing terms e Dependencies fSiteCollisionBGKGuo e Comments This routine is fundamental to Lattice Boltzmann calculations and should not be modified fCollisionBGKGuolncom e Header records int CollisionBGKGuoIncom e Function Collision steps for all incompressible fluids using BGK model with Guo forcing terms e Dependencies fSiteCollisionBGKGuolncom e Comments This routine is fundamental to Lattice Boltzmann calculations and should not be modified fCollisionMRT e Header records int fCollisionMRTO e Function Collision steps for all compressible fluids using MRT model e Dependencies fSiteCollisionMRT e Comments This routine is fundamental to Lattice Boltzmann calculations and should not be modified 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 97 fCollisionMRTIncom e H
156. initially swapped with their conjugate values f 7 then at each point the value f Z is then swapped with f 7 At These sets of swaps can be carried out either in two separate steps or in one go The use of separate swap steps requires two sweeps through the domain but the order in which the distribution functions are swapped does not matter and no boundary domain is necessary for serial calculations Simultaneous swapping cannot make use of automatic periodic boundary conditions thus a non zero domain boundary is required and requires lattice links to be additionally ordered so that the first half are all directed to lattice points that have previously gone through at least the first swap stage but only a single sweep through the domain is required By default the serial version of DL_MESO_LBE uses the propagation routine that carries out two separate swap steps fPropagationSwap while the parallel version uses the combined swap propagation routine fPropagationCombSwap If the user wishes to use the combined swap propagation routine in serial an al ternative program slbecombine cpp is available which includes the following modifications to slbe cpp e Replacement of fSetSerialDomain with fSetSerialDomainBuffer e Addition of fsBoundPeriodic and fMarkBoundArea immediately before fInitializeSystem e Addition of fsPeriodic inside main loop before calculating interaction forces it does not matter whether this is placed befo
157. ion in a class of nonideal gas lattice Boltzmann models Physical Review E 77 6 066702 June 2008 Xiaowen Shan and Hudong Chen Lattice Boltzmann model for simulating flows with multiple phases and components Physical Review E Statistical Physics Plasmas Fluids and related interdisciplinary topics 47 3 1815 1819 March 1993 Xiaowen Shan and Hudong Chen Simulation of nonideal gases and liquid gas phase transitions by the lattice Boltzmann equation Physical Review E Statistical Physics Plasmas Fluids and related interdisciplinary topics 49 4 2941 2948 April 1994 W Smith Molecular dynamics on hypercube parallel computers Computer Physics Communications 62 2 3 229 248 March 1991 W Smith A replicated data molecular dynamics strategy for the parallel Ewald sum Computer Physics Communications 67 3 392 406 January 1992 W Smith Calculating the pressure CCP3 Information Quarterly 39 14 20 October 1993 W Smith T R Forester and I T Todorov The DI POLY Classic user manual STFC STFC Daresbury Laboratory Daresbury Warrington Cheshire WA4 4AD United Kingdom version 1 0 edition December 2010 Simeon D Stoyanov and Robert D Groot From molecular dynamics to hydrodynamics A novel Galilean invariant thermostat Journal of Chemical Physics 122 11 114112 2005 Sauro Succi The Lattice Boltzmann Equation for Fluid Dynamics and Beyond Clarendon Press Oxford 2001 Michael R Swift W R Osborn an
158. irectory as the executable the number of species and their names will be read from it otherwise the user will be asked to enter this information For each molecule the user is asked for the side length for the cube bond length and number of molecule chains followed by the number of particles for each chain If the chain in question is not the first primary chain the user will also be asked for a pre existing bead number as the starting point for the chain and the neighbouring bead for determining bond angles After this the default species for the beads in the molecule will be requested the user will then be asked enter the bead numbers for each of the other species 0 can be entered to finish specifying bead numbers MOLECULE and FIELD files will be generated It should be noted that as generated by the utility e no bond angle or dihedral interaction parameters are included in the FIELD file B 2 DL_MESO_DPD 189 e the bond and bond angles are all of types 1 e no bond dihedrals are included Some editing by the user may be required before the files can be used in DPD simulations Section 11 3 1 gives the formatting of these files and the example cases involving molecules in the DEMO DPD directory should be examined by the user for guidance If the data for additional molecule types are to be included they should be appended to both the MOLECULE and FIELD files and the molecules must be in the same order in both files up to eight
159. irmed by the plot of horizontal velocity component as a function of vertical position at the last time step for the simulation 105 106 CHAPTER 7 DL_MESO LBE EXAMPLES J gt a 0 001 0 0005 4 3 a E A Bo El 2 3 E 0 0005 4 oti l i l i l 20 0 0001 0 0002 0 0003 0 0004 Vertical position y a Vector plot of system b Variation of horizontal velocity with vertical posi tion Figure 7 2 Results from LBE 2D _Shear test case 7 3 2D_CylinderFlow This is a 2D simulation of a single fluid on a 125 x 50 grid with a constant horizontal body force and a circular obstacle of radius 12 representing channel flow past an infinitely long cylinder Figure 7 3 shows this flow against a map of the system red representing the fluid and blue indicating solid objects the walls and cylinder Figure 7 3 Density and velocity vector plot from LBE 2D_CylinderFlow test case 7 4 2D KarmanVortex This is a 2D simulation of a single fluid on a 250 x 50 grid with a constant horizontal body force and a circular obstacle of radius 8 representing channel flow past an infinitely long cylinder that eventually produces a von K rm n vortex street between two solid walls Figure 7 4 shows the flow field at the final time step an AVI video file has been rendered from the calculation and can be found in the subdirectory for this case Figure 7 4 Velocity magnitude plot f
160. is Licence Agreement the website with the URL http www ccp5 ac uk DL_MESO and any website that from time to time replaces that website any virus worm time bomb time lock drop dead device trap and access code or anything else that might disrupt disable harm or 193 194 APPENDIX C DL_MESO LICENCE AGREEMENT ACADEMIC PURPOSES impede the operation of any information system or that might corrupt damage destroy or render inaccessible any software data or file on or that may allow any unauthorised person to gain access to any information system or any software data or file on it Intellectual Property patents trade marks service marks registered designs copyrights database rights design rights know how confidential information applications for any of the above and any similar right recognised from time to time in any jurisdiction together with all rights of action in relation to the infringement of any of the above the Licence Period the period beginning when the Licensee agrees to the terms and conditions of this Licence Agreement and downloads the DL_MESO Software from the DL_MESO Website and ending on the termination of this Licence Agreement under clause 5 2 2 LICENCE 2 1 STFC grants the Licensee an indefinite non exclusive non transferable royalty free licence to use copy modify and enhance the DL_MESO Software during the Licence Period on the terms and conditions of this
161. is generated and maintained by the moving boundaries in the planes normal to the y axis and a plot of the horizontal component of fluid velocity against vertical position at the last time step 0 001 T r T 0 0005 Z E 5 E 3 1 BOF F 3 x 3 E 0 0005 0 001 1 0 0 0001 0 0002 0 0003 Vertical position y a Vector plot of system b Plot of horizontal velocity with vertical position Figure 7 7 Results from LBE 3D_Shear test case Part II Dissipative Particle Dynamics DPD 109 Chapter 8 Dissipative Particle Dynamics Basic Theory 8 1 Introduction Dissipative Particle Dynamics DPD is an off lattice discrete particle method for modelling mesoscopic systems It has little in common with Lattice Boltzmann methods except in its application to systems of similar length and time scales The DPD method inherits its methodology from classical Molecular Dynamics MD particularly from Brownian Dynamics BD It differs from BD however in an important way it is Galilean invariant and for this reason conserves hydrodynamic behaviour while the BD method does not Many systems are crucially dependent on hydrodynamic interactions and it is essential to retain this feature in the model DPD is particularly useful for simulating systems on the near molecular scale such as polymers biopolymers lipids emulsions and surfactants systems in which large scale structure evolves on a time
162. ive e Arguments idnode input integer nodes input integer statis e Header records SUBROUTINE statis e Function Calculates statistical properties of system including rolling averages and fluctuations e Dependencies sclsum global_sum_dble 113 DL_MESO_DPD Data Files 11 3 1 Input files All input files for DL MESO_DPD must be in ANSI text format with keywords where necessary and numerical values separated from each other with spaces or commas tabs are currently not recognised by the parsing utilities Define system CONTROL The CONTROL file can be obtained either via use of the DL_MESO GUI or by editing existing files of that name such as those in the DEMO DPD directory With the exception of the simulation name on the first line its format consists of keywords for system parameters of which the first four characters are recognised by DLLMESO_DPD followed by their values with spaces separating them The majority of keywords are followed by a single value these are summarized in Table 11 1 While most parameters are defined by a single value certain keywords are followed by more than one parameter e essential data name number mass interaction parameters constant external forces for individual species SPEC name Ni Mi A B C yi Ti Fyri Pyg Ei 1If this is specified the system is assumed to be cubic 2 Takes the largest out of the maximum bond length interaction cutoff rcut or electrostatic cutoff rele
163. lata 3 19 X i eiaeigeiy Paguy Payup Paya plates ty 3 20 a 24 CHAPTER 3 THE LATTICE BOLTZMANN EQUATION BASIC THEORY With these definitions Equation 3 17 leads to the weakly compressible Navier Stokes equation wD Sap 300Pa D 3 pua Og puaug OgPag Og 29 puy Oo pug 0x0 0 0 3 21 where the kinetic viscosity is given by y E 3 ry 3 Similarly the convection diffusion equation governing the evolution of solute transfer and the convection conduction equation governing the evolution of thermal transfer can be obtained 3 5 Mesoscale Interaction The rate of change of momentum is proportional to the force Pua t At pug t Fy At 3 22 The Lattice Boltzmann Equation takes into account the effect of relaxation and expresses Equation 3 22 in a new format 32 Pua t At pua t Fate 3 23 F with F as the component in the a direction can be a long ranged body force or any local interactions In the Shan Chen model for multiple phases and components 43 the force on component a is defined as Fe po 2 Y gan Y wiy Z 3 24 b 7 where ga is the interaction coefficient between elements a and b 4 is related to the density of element a and can take many different forms e g y 5 1 exp p 3 25 3 6 Summary of Lattice Boltzmann Equation Lattice Boltzmann is an established numerical methodology for handling hydrodynamics in fl
164. le uwall Function Calculates the particle distribution function at a fixed density boundary for compressible fluids e Dependencies Many for different lattice schemes 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 81 e Arguments tpos input long integer prop input long integer uwall output array of doubles e Comments Planar surface calculations are based on 60 concave edges and corners assume zero speed at the boundary The array uwall is the velocity at the grid point for all fluids which is subsequently used for solute concentration and temperature boundary conditions fFixedDensityFluidIncom e Header records int fFixedDensityFluidIncom long tpos int prop double uwall e Function Calculates the particle distribution function at a fixed density boundary for incompressible fluids e Dependencies Many for different lattice schemes e Arguments tpos input long integer prop input long integer uwall output array of doubles e Comments Planar surface calculations are based on 60 concave edges and corners assume zero speed at the boundary The array uwall is the velocity at the grid point for all fluids which is subsequently used for solute concentration and temperature boundary conditions fFixedSoluteConcen e Header records int fFixedSoluteConcen long tpos int prop double uwall e Function Calculates the particle distribution function at a fixed composition boundary for compressible and incom
165. lel running than only calculating the bonds in each node but ensures longer bond lengths can be accommodated 11 2 10 domain module This module requires the modules constants variables error_module comms_module if running parallel version and bond module if running parallel version to be loaded beforehand domain e Header records SUBROUTINE domain 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 157 e Function Determines 3D domain decomposition for system number of nodes in each direction location for each node in system nearest neighbouring nodes e Dependencies parallel version only msg_receive_blocked msg_receive_sca_blocked msg_send_blocked msg_send_sca_blocked e Comments Essential for parallel version of DL_MESO_DPD the serial version of this subroutine sets values for a single processor parlnk e Header records SUBROUTINE parlnk numi num2 e Function Constructs the parallel link cells for calculations of pairwise forces between particles e Dependencies None e Arguments numi input integer num2 input integer e Comments Pairwise electrostatic forces typically act over longer lengthscales and are not considered using this sub routine a similar routine for the real space part of Ewald summation is included in ewald_real deport e Header records SUBROUTINE deport nlimit mdir mp begin final shove e Function Deports particles from boundary halo to neighbouring domain
166. libriumC fGetEquilibriumT e Arguments startpos input double pointer bodyforce input double pointer e Comments Collisions are carried out using the BGK single relaxation time 4 ta e F FE fSiteCollisionBGKIncom e Header records int fSiteCollisionBGKIncom double startpos double bodyforce e Function Calculates collisions at a grid point using the Bhatnagar Grook Krook BGK model for incompressible fluids e Dependencies fGetOneMassSite fGetSpeedIncomSite fGetEquilibriumincom fGetEquilibriumC fGetEquilibriumT e Arguments startpos input double pointer bodyforce input double pointer 92 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE e Comments Collisions are carried out using the BGK single relaxation time 4 ao Er Fi 72 fSiteCollisionBGKGuo e Header records int fSiteCollisionBGKGuo doublex startpos double bodyforce e Function Calculates collisions at a grid point using the Bhatnagar Grook Krook BGK model with the Guo forcing term for compressible fluids e Dependencies fGetOneMassSite fGetSpeedSite fGetEquilibriumF fGetEquilibriumC fGetEquilibriumT e Arguments startpos input double pointer bodyforce input double pointer e Comments Collisions are carried out using the BGK single relaxation time 4 and the Guo forcing term 14 0 tia0a fi fi fi 1 wi 3 eia Va 9 8 Fa fSiteCollisionBGKGuolncom e Header records int fSiteCollisionBGKGuoIncom d
167. lows 1 ee p pho q E n 8 11 j gt i 2T 5 R r 2 pkBT p e rA ta glr r dr 8 12 where g r is a radial distribution function for the soft sphere model 13 and p is the DPD particle density For sufficiently large densities p gt 2 g r takes the same form and the equation of state can be well approximated by p pkgT aAp 8 13 with the parameter a 0 101 0 001 This expression permits the use of fluid compressibilities to obtain conservative force parameters for bulk fluids e g for water A a Alternative equations of state may be obtained by modifying the functional form of conservative interactions to include localized densities i e many body DPD 37 53 Transport coefficients for a DPD fluid can be derived using the expressions for the drag and stochastic forces 13 25 31 The kinematic viscosity can be found to be MBkgT 27ypr x 8 14 4rypr 1575 6 while the self diffusion coefficient is given as 45kpT PEA 8 15 2rypr The ratio of these two properties the Schmidt number Sc is therefore 1 2 4y2 Sem 1 TIe 8 16 2 70875kgT and for values of the drag coefficient and density frequently used in DPD simulations this value is of the order of unity which is an appropriate magnitude for gases but three orders of magnitude too small for liquids This property of standard DPD does not rule it out for simulations of liquid phases except when hydrodynamics are important It ma
168. ltiple relaxation time MRT schemes the call to C011isionBGKx should be changed to fCollisionMRT and the tuneable parameters can be modified in the 1bpMODEL cpp file for the required lattice model D2Q9 D3Q15 or D3Q19 Additional boundary conditions should be added to the fPostCollBoundary or fPostPropBoundary routines depending on whether they are applied after collisions and before propagation or after propagation respectively lif preferred an alternative version that uses a boundary layer slbecombine cpp is available 2 may need to be replaced by on computers running Windows 183 184 APPENDIX A MANUAL COMPILIATION AND RUNNING OF DL_MESO Calls to alternative and or additional output routines can be included in the main program files provided they are called before the output file number qVersion is increased If no fluid fluid interaction forces are required the call to fInteractionForceShanChen can be commented out to increase the calculational efficiency A 2 DL MESO DPD The Fortran90 modules for DL_MESO_DPD must be compiled in a particular order to satisfy dependencies of shared variables and arrays e constants e variables e numeric_container e comms_module e error module e parse_utils e bond_module e surface_module e ewald_ module e manybody_module e domain_module e start_module e config module e field module e integrate module e statistics module e run module e dimesodpd
169. mat trajectory file traject vtf that can visualize the simulation after equilibration such that snapshots at the recorded timesteps and animations can be produced Two versions of the utility are provided traject f90 outputs every particle for all recorded timesteps and trajectselected f90 allows the user to select the number of particles and the number of timesteps to output to the trajectory file Both versions of the utility can be used with both HISTORY and HISTORY files produced by the serial and parallel versions of DL_MESO_DPD respectively If 90 is the command for the available Fortran90 compiler the executable traject exe can be produced by typing e 90 o traject exe traject f90 and either run at the command line or via the GUI in Process DPD Data A command line argument indicating the number of processes and therefore the number of HISTORY or HISTORY files can be included in a similar way to the exportimage utility The alternative version can be compiled in the same way by default the makefiles create executables named trajects exe This version does not take command line arguments and all information must be typed in by the user After the number of processes is typed in the utility displays the total number of particles and number of unbonded particles in the simulation before asking for the particle number range and number of timesteps to be included in the trajectory file This version of the utility cannot be invoked usi
170. mension nlewx integer number of link cells for electrostatics in domain cell y dimension nlewy integer number of link cells for electrostatics in domain cell z dimension nlewz integer number of link cells for electrostatics in domain cell and boundary halo x dimension nlewx2 integer number of link cells for electrostatics in domain cell and boundary halo y dimension nlewy2 integer number of link cells for electrostatics in domain cell and boundary halo z dimension nlewz2 integer electrostatic link cell length in x direction wdthewx real KIND dp electrostatic link cell length in y direction wdthewy real KIND dp electrostatic link cell length in z direction wdthewz real KIND dp species name namspe character 8 A mxspe species particle mass amass real KIND dp A mxspe potential interaction type ktype integer species interaction parameter energy A aaa real KIND dp A mxspe species interaction parameter distance B bbb real KIND dp A mxspe species interaction parameter charge C ccc real KIND dp A mxspe species interaction parameter dissipation D y or T ddd real KIND dp A mxspe Lennard Jones long range potential correction clr real KIND dp A 2 interaction parameter storage vvv real KIND dp A mxprm mxpot unlike particle A value unlikeaaa real KIND dp A mxpot unlike particle B value unlikebbb real KIND dp A mxpot unlike particle y or Ty value unlikegam real KIND dp A mxpot species number 7 for unlike
171. n this Licence Agreement are for ease of reference only they do not affect its construction or interpretation Assignment etc The Licensee may not assign or transfer this Licence Agreement as a whole or any of its rights or obligations under it without first obtaining the written consent of STFC Illegal unenforceable provisions If the whole or any part of any provision of this Licence Agreement is void or unenforceable in any jurisdiction the other provisions of this Licence Agreement and the rest of the void or unenforceable provision will continue in force in that jurisdiction and the validity and enforceability of that provision in any other jurisdiction will not be affected Waiver of rights If STFC fails to enforce or delays in enforcing an obligation of the Licensee or fails to exercise or delays in exercising a right under this Licence Agreement that failure or delay will not affect its right to enforce that obligation or constitute a waiver of that right Any waiver by STFC of any provision of this Licence Agreement will not unless expressly stated to the contrary constitute a waiver of that provision on a future occasion Entire agreement This Licence Agreement constitutes the entire agreement between the parties relating to its subject matter The Licensee acknowledges that it has not entered into this Licence Agreement on the basis of any warranty representation statement agreement or undertaking except
172. n z npz integer 2 1 list of neighbouring processes map integer A6 x position of domain cell idx integer y position of domain cell idy integer z position of domain cell idz integer position of domain cell origin in system volume x dimension delx real KIND dp position of domain cell origin in system volume y dimension dely real KIND dp position of domain cell origin in system volume z dimension delz real KIND dp domain cell length in x direction sidex real KIND dp domain cell length in y direction sidey real KIND dp domain cell length in z direction sidez real KIND dp number of particles in domain cell nbeads integer number of link cells in domain cell x dimension nlx integer number of link cells in domain cell y dimension nly integer number of link cells in domain cell z dimension nlz integer number of link cells in domain cell and boundary halo x dimension n1x2 integer 120 CHAPTER 9 DL_MESO_DPD BASIC DEFINITION Table 9 4 DL MESO_DPD Parameters continued function parameter data type notation number of link cells in domain cell and boundary halo y dimension nly2 integer number of link cells in domain cell and boundary halo z dimension nlz2 integer link cell length in x direction wdthx real KIND dp link cell length in y direction wdthy real KIND dp link cell length in z direction wdthz real KIND dp number of link cells for electrostatics in domain cell x di
173. nction Defines vector messages for system distribution functions boundary properties and interaction forces fDefineNeighbours e Header records int fDefineNeighbours e Function Calculates the names of neighbouring processes and the start points for sending and receiving messages e Comments This subroutine must not be changed fNonBlockCommunication e Header records int fNonBlockCommunication e Function Passes distribution function information for either 2D or 3D system fOutputInfo e Header records int fOutputInfo e Function Outputs number of processes and lengths of integers and floats e Comments This subroutine is necessary for gathering and rearranging the 1bout data and produces the files lbout info and lbout ext fBoundNonBlockCommunication e Header records int fBoundNonBlockCommunication e Function Passes boundary information for either 2D or 3D systems fForceNonBlockCommunication e Header records int fForceNonBlockCommunication e Function Passes interaction force information for either 2D or 3D systems 6 3 DL MESO_LBE DATA FILES 101 fCheckTimeMPI e Header records int fCheckTimeMPI O e Function Outputs time in seconds since initial call e Arguments fCheckTimeMPI output double e Comments Obtains calculation time based on MPI wall clock fPrintSystemMass e Header records int fPrintSystemMass e Function Calculates and print
174. ng Co1lisionBGKGuo 5 1 1 2 Multiple Relaxation Time MRT The MRT collision operator operates on a similar principle to the quasilinear Lattice Boltzmann Equation 18 which expresses collisions in terms of a square matrix with dimensions equal to the number of lattice links per grid point m 1bsy nq Unlike quasilinear LBE however MRT collision schemes are applied to the moments for each lattice point rather than the distribution functions 26 7 which are related to each other by M Tf 5 3 where f is a vector of all m distribution functions for the point ie fo fi seat M the vector of moments also of size m and dependent on the lattice system and T the transformation matrix that renders the moments in terms of the distribution functions Equilibrium values for the moments M 4 can be determined by transforming the standard local equilibrium distribution function into moment space by Met T fe 5 4 where fea is the vector of local equilibrium distribution functions the resulting equilibrium moments can alternatively be expressed directly as functions of fluid density and velocity Certain moments such as density and momentum must be conserved and their equilibrium values are set so that no changes are made The post collisional moments are determined by relaxation of the non equilibrium part i e M t M d t A M E t M 4 z t 5 5 where A is the collision matrix which takes the form of a diagonal ma
175. ng with HISTORY every nsbsv time steps Each binary file contains some of the system properties including an identity number for each molecule followed by the positions velocities and virial for each particle in the domain The HISTORY files can be used with the utility traject in the DPD utility directory to produce plottable VMD files with sets of bonded particles represented as residues The utility local in the same directory can calculate localized properties e g temperature composition pressure from the same files and produces VTK files with these properties as cell data Appendix B gives instructions for their use Statistical data file CORREL If the parameter 1corr is set to true DL MESO_DPD will generate an ANSI text file containing statistical data every iscorr time steps which can later be imported into a spreadsheet or used by graph plotting software The data is formatted thus t Etot Epot tot Ebpot elec Epot bond Ebpot angle Epot dihed PVT Tbond Bangle dined where t is the time E energy tot denoting total pot potential P pressure V volume T temperature rpona the mean bond length Oangie the mean bond angle and ainea the mean bond dihedral The mean bond angle and dihedral are expressed in degrees Chapter 12 DL MESO DPD Examples Test cases for Dissipative Particle Dynamics simulations using DL_MESO including the required input and sample output files can be found in the DEMO DP
176. ng the GUI local local is a utility written in Fortran90 that can read in HISTORY output data files generated by DL_MESO_DPD and produce series of VTK format files containing statistical properties number of beads density compositions per particle and molecule types pressure temperature and mean velocity in cuboidal subdivisions of the simulation volume for plotting and or visualization The source code for this utility local 90 can read both HISTORY and HISTORY files generated by the serial and parallel versions of DL_MESO_DPD respectively If 90 is the command for the available Fortran90 compiler the executable local exe can be produced by typing e 90 o local exe local f90 and either run at the command line or via the GUI in Process DPD Data after entering both the number of processes used in simulations and the number of divisions required in each dimension This utility can be run with four command line arguments the first denoting the number of processes and the others giving the number of divisions required in x y and z directions respectively If these arguments are omitted the user will be asked to enter these values The following files are produced for all the specified time steps after equilibration e scalars_ vtk with the numbers of beads density pressure and temperature for each cuboid cell B 2 DL_MESO_DPD 191 e velocity_ vtk with the mean velocity of the beads within the cells e spec
177. ng unit being allocated a subset of bonds to deal with including bonds across neighbouring units or 2 globally with all units holding all bond data and sharing bonded particle positions each carrying out all bond calculations and appropriately allocating forces to local particles The former method may require larger boundary halo sizes for the bond lengths being simulated but is more efficient for larger numbers of molecules and processing units while the latter method requires the sharing of information between all units but does not require halo information and is guaranteed to find all bonds 123 124 CHAPTER 10 DL MESO DPD FEATURES Bookkeeping arrays bndtbl angtbl and dhdtb1 list all particles involved in bonded interactions according to global index numbers and point to appropriate arrays of parameters to define the potential If the key bonded particle for a bond moves from one processing unit to another the entry in the bookkeeping array is also moved At each time step a list of bonded particles in each domain 1b1c1st is created to relate global index numbers to the local index numbers used by the processing unit in force velocity and coordinate arrays This global local index list is sorted by global index number to allow cross referencing to local index numbers by means of a binary search 10 1 2 Electrostatic interactions For systems with periodic boundary conditions DL_MESO_DPD uses the Ewald Sum to calculate Coul
178. nput array of characters e Comments The default file name is lbin sys to use different file names the argument for this routine can be changed by the user fReadSpaceParameter e Header records int fReadSpaceParameter const char filename lbin spa e Function Reads 2D or 3D space parameters from data file bin spa e Dependencies lbin spa data file fReadSpace2D fReadSpace3D e Arguments filename input array of characters e Comments The default file name is lbin spa this can be changed by the user if the space data file has a different name fReadInitialState e Header records int fReadInitialState const char filename lbin init e Function Reads initial velocities densities concentrations and temperatures from data file 1bin init and calculates initial distribution functions for compressible fluids 66 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE e Dependencies lbin init data file fReadInitialState2D fReadInitialState3D Arguments filename input array of characters Comments The default file name is lbin init this can be changed by the user if the initial state data file has a different name This routine will replace the default values of all properties at specified points and should be called after fInitializeSysten fReadInitialStateIncom Header records int fReadInitialStateIncom const char filename lbin init Function Reads initial velocities densities concentrations an
179. ny body interactions 2 species i species j Aij Bij Cij Dij Eij add mb reorder mb reset mb savep CANCELP b Potential parameters fj DPD barostat properties ej la les LANGEVIN barostat relaxation time 0 0 2 piston drag coefficient 0 0 isotropic system 3 SAVE B CANCEL B d Barostat parameters El DPD surface properties aj EI les wall directions SAVE SF CANCEL SF f Surface parameters HARD SURFACES WITH SOFT REPULSION 1 Ex By Gz 2 surface repulsion species 1 0 0 3 THE DL_MESO GUI Figure 2 6 Species potential thermostat barostat electrostatic and surface pop up windows 3 Particle mass for the species 4 Parameter A energy this parameter is equal either to the conservative interaction parameter A Aii for DPD beads or the potential well depth e ei for hard sphere systems i e Lennard Jones Weeks Chandler Anderson 5 Parameter B distance this parameter is equal either to the particle radius B for DPD or zero potential distance 0 for hard sphere models 2 3 DISSIPATIVE PARTICLE DYNAMICS AND THE DL_MESO GUI 15 6 Parameter C charge this parameter is equal to the net charge of the particle relative to the charge on an electron e This parameter does not have to be specified for systems without electrostatic interactions 7 Parameter D di
180. ode e Function Creates a Gaussian random number e Dependencies mtrnd e Arguments idnode input integer gaussmp output real KIND dp e Comments This is an implementation of the Marsaglia polar method 29 to convert linear random numbers generated by the Mersenne Twister method to Gaussian random variables with zero mean and unity variance sclsum e Header records REAL KIND dp FUNCTION sclsum n a i e Function Calculates the scalar sum of an array e Dependencies None e Arguments n input integer a input array of real KIND dp i input integer sclsum output real KIND dp 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 143 erfcdp e Header records REAL KIND dp FUNCTION erfcdp x e Function Calculates the complementary error function for x erfc x e Dependencies None e Arguments x input real KIND dp erfcdp output real KIND dp e Comments This approximation for the function is based on a Chebyshev polynomial fitting 15 11 2 3 comms module This module is essential for parallel running and does not require detailed knowledge for its use depending on the version and implementation of MPI available the user may wish to select between the lines USE MPI and INCLUDE mpif h for loading the necessary routines The serial version of the comms_module primarily consists of dummy routines to satisfy the required calls in the rest of the code initcomms e Header records SUBROUTINE ini
181. odic box undergoing time evolution as a result of the above forces In implementation it differs very little from Molecular Dynamics It should be noted that all computed interactions are pairwise which means that the principle of the conservation of momentum in the system or Galilean invariance is preserved The conservation of momentum is required for the preservation of hydrodynamic forces Chapter 9 DL MESO DPD Basic Definition 9 1 Data structure 9 1 1 Storage of running information DL_MESO_DPD contains storage for information on the system being modelled the domain and neighbour information for parallel running The parameters for these aspects of calculations can be found in Tables 9 1 9 2 and 9 3 It should be noted for the parameters in Table 9 1 that the fundamental units for the simulation are those of mass M length L and energy E the DPD unit of time is equivalent to L 4 E while temperature in the form kgT is defined as two thirds of the kinetic energy per particle 115 116 CHAPTER 9 DL_MESO_DPD BASIC DEFINITION Table 9 1 System information parameter meaning text name of calculation volm initial volume of system dimx dimy dimz initial dimensions of system nsyst total number of particles nusyst total number of unbonded particles nspe number of particle species nmoldef number of molecule types temp specified system temperature kgT prszero specified system pr
182. of bonds for molecule type nbond integer A nmoldef number of bond angles for molecule type nangle integer A nmoldef number of bond dihedrals for molecule type ndihed integer A nmoldef bond table storage for molecule insertion bdinp1 integer A nmoldef mxbonds bond table storage for molecule insertion bdinp2 integer A nmoldef mxbonds bond table storage for molecule insertion bdinp3 integer A nmoldef mxbonds angle table storage for molecule insertion anginp1 integer A nmoldef mxbonds angle table storage for molecule insertion anginp2 integer A nmoldef mxbonds angle table storage for molecule insertion anginp3 integer A nmoldef mxbonds angle table storage for molecule insertion anginp4 integer A nmoldef mxbonds dihedral table storage for molecule insertion dhdinp1 integer A nmoldef mxbonds dihedral table storage for molecule insertion dhdinp2 integer A nmoldef mxbonds dihedral table storage for molecule insertion dhdinp3 integer A nmoldef mxbonds dihedral table storage for molecule insertion dhdinp4 integer A nmoldef mxbonds dihedral table storage for molecule insertion dhdinp5 integer A nmoldef mxbonds many body DPD interaction parameter Ajj aamb real KIND dp A mxpot many body DPD interaction parameter Bij bbmb real KIND dp A mxpot many body DPD interaction parameter Cj ccmb real KIND dp A mxpot many body DPD interaction parameter D 5 ddmb real KIND dp A mxpot many body DPD interaction parameter Ej eemb real KI
183. om fGetMRTCollide fGetEquilibriumC fGetEquilibriumT e Arguments startpos input double pointer bodyforce input double pointer e Comments Collisions are carried out using multiple relaxation time MRT schemes 26 7 Os eia a fi T7 3 T7 Mis fSiteCollisionMRTGuo e Header records int fSiteCollisionMRTGuo doublex startpos double bodyforce 94 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE e Function Calculates collisions at a grid point using Multiple Relaxation Time MRT models coupled with Guo like forcing terms for compressible fluids e Dependencies fGetOneMassSite fGetSpeedSite fGetMomentEquilibriumF fGetMomentForce fGetMRTCollide fGetEquilibriumC fGetEquilibriumT e Arguments startpos input double pointer bodyforce input double pointer e Comments Collisions are carried out using multiple relaxation time MRT schemes coupled with Guo like forcing terms 39 8 ei a fi T 5 TF Mf 1 313 5 fSiteCollisionMRTGuolncom e Header records int fSiteCollisionMRTGuoIncom double startpos double bodyforce e Function Calculates collisions at a grid point using Multiple Relaxation Time MRT models coupled with Guo like forcing terms for incompressible fluids e Dependencies fGetOneMassSite fGetSpeedIncomSite fGetMomentEquilibriumincom fGetMomentForce fGetMRTCollide fGetEquilibriumC fGetEquilibriumT e Arguments startpos input double pointer bodyforce input double pointer e
184. ombic interactions see Section 10 5 Calculation of the real space component in routine ewald_real uses the same link cell algorithm as for other pairwise interactions albeit using a larger cutoff radius re and requires a larger boundary halo than for standard pairwise interactions 10 2 Thermostats and integration algorithms The integration algorithms in DL_MESO_DPD are based on the second order Velocity Verlet VV scheme 54 which yields the positions velocities and forces of particles at the same time and is generally used in molecular dynamics simulations This algorithm has two stages The first stage advances the particle velocities to time t At by integrating the forces and uses the new half step velocities to advance the position to time t At A 4 At F t Ti t At F t Ato t 4At 10 2 The positions at the end of the time step allow the forces to be recalculated before the second stage of the algorithm is applied to advance the half step velocities to the end of the time step by integrating with the new force 3 0 t At 0 t 4At a 10 3 Mi Five thermostating algorithms are currently available in DL_MESO_DPD two variants of the standard DPD thermostat and three alternative schemes which apply velocity corrections to the particles after force integration The algorithm can be selected in the CONTROL file using the keyword THERM itype is an integer specifying the thermostat integration algorithm
185. on Swaps pair of numbers a and b e Dependencies None e Arguments a input output any datatype b input output any datatype fGet Number Ordered e Header records two cases int fGetNumberOrdered int amp iox int amp ioy int amp ioz int fGetNumberOrdered int amp iox int amp ioy e Function Rearrange the integers in descending order e Dependencies None e Arguments iox input output integer reference ioy input output integer reference ioz input output integer reference 54 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE fGet Number OrderFixed e Header records two cases int fGetNumberOrderFixed int amp iox int amp ioy int amp ioz int ix int iy int iz int fGetNumberOrderFixed int amp iox int amp ioy int ix int iy e Function Rearrange a set of integers so they appear in the same order as another set of integers e Dependencies GetNumberOrdered e Arguments iox input output integer reference ioy input output integer reference ioz input output integer reference ix input integer iy input integer iz input integer fBestGrouping e Header records int fBestGrouping int totalgrid int totalgroup int indigrid int amp critigroup e Function Distribute grid points among the processes to give a maximum of one to the differences in the numbers of grid points e Dependencies None e Arguments totalgrid input integer totalgroup input integer indigrid output integer reference
186. on prior to the next collision step for incompressible fluids Dependencies Many for different lattice schemes Comments Algorithms for other boundary conditions can be added by the user although care should be taken as to when they are applied in each time step the conditions invoked in this routine apply to distribution functions after propagation fsPeriodic Header records int fsPeriodic Function Applies periodic boundary condition for serial calculations with non zero boundary domain widths by copying distribution functions from edges of fluid points Dependencies fsPeriodic2D fsPeriodic3D Comments Serial equivalent of NonBlockCommunication essential for using the combined swap propagation routine fPropagationCombSwap in serial running 84 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE fsBoundPeriodic e Header records int fsPeriodic e Function Applies periodic boundary condition for serial calculations with non zero boundary domain widths by copying boundary information from edges of fluid points e Dependencies fsPeriodic2D fsPeriodic3D e Comments Serial equivalent of BoundNonBlockCommunication may be required for using combined collision propagation routines Collision Propagation in serial running fsForcePeriodic e Header records int fsPeriodic e Function Applies periodic boundary condition for serial calculations with non zero boundary domain widths by copying interaction
187. on at the tpos th grid point e Dependencies None 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 79 e Arguments tpos input long integer e Comments This bounce back boundary condition is carried out by exchanging post collisional populations with con jugate values in neighbouring grid points fSiteBlankF e Header records int fSiteBlankF long tpos e Function Sets the fluid particle distribution function at the tpos th grid point to zero e Dependencies None e Arguments tpos input long integer e Comments This routine is used to ensure e g flows inside solid boundaries are negligible fSiteBlankC e Header records int fSiteBlankC long tpos e Function Sets the solute particle distribution function at the tpos th grid point to zero e Dependencies None e Arguments tpos input long integer e Comments This routine is used to ensure e g diffusion inside a bulk solid is negligible compared to that in a liquid fSiteBlankT e Header records int fSiteBlankT long tpos e Function Sets the temperature particle distribution function at the tpos th grid point to zero e Dependencies None e Arguments tpos input long integer e Comments This routine is used to ensure e g negligible heat transfer through an insulator 80 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE fFixedSpeedFluid e Header records int fFixedSpeedFluid long tpos int prop double uwall e Function Calculates the parti
188. on parameter d ddang real KIND dp A mxbonddef dihedral interaction parameter a aadhd real KIND dp A mxbonddef dihedral interaction parameter b bbdhd real KIND dp A mxbonddef dihedral interaction parameter c ccdhd real KIND dp A mxbonddef dihedral interaction parameter d dddhd real KIND dp A mxbonddef bond types bdtype integer A mxbonddef bond angle types angtype integer A mxbonddef bond dihedral types dhdtype integer A mxbonddef bond table bndtbl integer A numbond 3 bond angle table angtbl integer A numang 4 bond dihedral table dhdtbl integer A numdhd 5 global local particle number list lblclst integer A 4 nsyst nusyst 2 number of bonds in table nbonds integer number of bond angles in table nangles integer number of bond dihedrals in table ndiheds integer number of entries in global local particle number list nlist integer species population of bonded particles nspecmol integer A mxspe molecule name nammol character 8 A mxmoldef molecule type population nmol integer A mxmoldef bead numbers in molecule types nbdmol integer A mxmoldef species number for molecule insertion mlstrtspe integer A nmoldef mxmolsize x coordinate for molecule insertion mlstrtxxx real KIND dp A nmoldef mxmolsize y coordinate for molecule insertion mlstrtyyy real KIND dp A nmoldef mxmolsize z coordinate for molecule insertion mlstrtzzz real KIND dp A nmoldef mxmolsize cube size for molecule insertion cbsize real KIND dp A nmoldef number
189. ond Velocity Verlet step exp tiga E SAL At via t 1At St Fae Od a 1 ulAt 10 25 ga uf t At 4 0 The barostat force is then calculated using the same Gaussian random numbers for each iteration a n 1 n 1 2 n FEJD t At V Pa Po N 2 mi o puf W Cpp a 10 26 and the barostat velocity is recalculated FETO t At At 2W ula t At uga t 3At 10 27 Equations 10 25 to 10 27 are repeated until convergence in particle velocities is achieved This normally takes a few iterations per time step without requiring recalculation of particle forces or rescaling of particle coordinates 10 3 2 Berendsen barostat btype 2 The governing equation for the Berendsen barostat 2 is a simple differential equation for the pressure dP P Po 22 ne 10 28 dt Tp which can be solved to give a scaling factor for the simulation volume 7 t bAt Na t 1 Po Pa t 10 29 Tp 10 4 PARTICLE PARTICLE INTERACTIONS 129 where 8 is the isothermal compressibility of the system The exact value of this property is not critical to the algorithm since it relies on the ratio de P The barostat is implemented using a variant of the Velocity Verlet algorithm after the midstep velocities are determined the scaling factor for time t is used to modify the particle positions and resize the simulation volume Tia t At Nalt ria t Atvia t 3At 10 30 Lalt At nalt L
190. onding CONTROL file as should the information in the FIELD and MOLECULE files for any bonded particles The periodic boundary key imcon 52 is ignored since all DL_MESO systems are orthorhombic initialvelocity e Header records SUBROUTINE initialvelocity e Function Initializes particle velocities randomly in system to give required system temperature e Dependencies duni global_sum_dble boxfit e Header records SUBROUTINE boxfit side idd ds displ npls e Function Sets up lattice for unbonded particles in processor domain e Dependencies None e Arguments side input real KIND dp idd input integer ds input real KIND dp displ output real KIND dp npls output integer 11 2 12 config module This module requires the modules constants variables comms module numeric_container parse_utils domain_module and surface_module to be loaded beforehand sysdef e Header records SUBROUTINE sysdef e Function Reads in system data determines simulation properties and prints to OUTPUT file 166 CHAPTER 11 Dependencies datsys bondsys domain surfacenodes forgen elecgen error global_sca_and Arguments idnode input integer datsys Header records SUBROUTINE datsys Function Reads in system data from CONTROL file Dependencies CONTROL file getdble getint getword error global_sca_and bondsys zero Header records SUBROUTINE bondsys Function DL_MESO_DPD PACKAGE REFERENCE
191. ons The DD strategy which underpins DL_ MESO_DPD is based on the link cell algorithm 19 which requires a relatively short ranged cutoff for interparticle potentials and forces There is a need for processing units to exchange halo data i e sending the contents of link cells at the boundaries of each domain to neighbouring units so each may have all the necessary information to compute pairwise forces acting on the particles in its allotted domain Similarly the force and virial contributions from particles in boundary halos need to be returned to their original processing units for summation The link cell algorithm is also applied in serial by duplicating system data to create the boundary halo across periodic boundaries The size of the boundary halo rhalo which can be specified in the CONTROL file should not be greater than the minimum system dimension per domain for good parallel performance it is recommended that the halo size should be no larger than one third of the smallest subdomain dimension The value of maxdim in the constants module giving the maximum sizes of force velocity and position arrays should be large enough to hold all particles in each domain plus any particles in boundary halos including duplicates when running in serial or using smaller numbers of processing units 10 1 1 Intramolecular interactions Intramolecular interactions may be handled in two different ways either 1 locally with each processi
192. onships between particles for each molecule type follow which are preceded by the line indicating the numbers of bond pairs angle triples and dihedral quadruples to follow if no angles or dihedrals are to be listed these values can be omitted MOLECULE Npond Nangle Ndihedral Bond pairs then follow in the format m1 1 11 2 by m2a 1 2 2 ba Mhn 1 Mn 2 bn where m j is the particle number of the molecule for the ith bond pair and the jth particle in the pair and b is an integer indicating the bond as previously listed e g 2 uses the second set of listed bond properties for pair 7 Similarly bond angle triples are formatted thus M 71 2 Mig Q1 M21 M22 M23 Q2 Mn 1 Mn 2 Mn 3 An while bond dihedral quadruples are formatted m11 M12 M13 M1 4 d M21 M22 M23 M2 4 d2 Mn 1 Mn 2 Mn 3 Mn 4 dn with a and d indicating the angle and dihedral properties to use for the triple quadruple Both of these files can be created using the supplied C program molecule generate in the directory DPD utility see Appendix B for more details although some editing is currently required to specify the bond angle dihedral types and properties Example files in the DEMO DPD directory can be examined for this purpose Initial state CONFIG An optional CONFIG file can be included to define the initial state of the system which can include the positions velocities and forces for each particle Before this information is inclu
193. oublex startpos double bodyforce e Function Calculates collisions at a grid point using the Bhatnagar Grook Krook BGK model with the Guo forcing term for incompressible fluids e Dependencies fGetOneMassSite fGetSpeedIncomSite fGetEquilibriumincom fGetEquilibriumC fGetEquilibriumT e Arguments startpos input double pointer bodyforce input double pointer e Comments Collisions are carried out using the BGK single relaxation time 4 and the Guo forcing term 14 8i eia a fi E fi F 1 2 wi lB ein 0 9G D Fa Tf 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 93 fSiteCollisionMRT e Header records int fSiteCollisionMRT double startpos double bodyforce e Function Calculates collisions at a grid point using Multiple Relaxation Time MRT models for compressible fluids e Dependencies fGetOneMassSite fGetSpeedSite fGetMomentEquilibriumF fGetMRTCollide fGetEquilibriumC fGetEquilibriumT e Arguments startpos input double pointer bodyforce input double pointer e Comments Collisions are carried out using multiple relaxation time MRT schemes 26 7 Os eia a fi TO 3 T7 M fSiteCollisionMRTIncom e Header records int fSiteCollisionMRTIncom double startpos double bodyforce e Function Calculates collisions at a grid point using Multiple Relaxation Time MRT models for incompressible fluids e Dependencies fGetOneMassSite fGetSpeedIncomSite fGetMomentEquilibriuminc
194. ous run using export files a new run which takes a starting state particle positions and velocities for a new simulation from export files and rescaled does the same as a new run but additionally rescales the particle velocities to give the specified system temperature The number of species involved in the simulation a maximum of 6 species can be specified in DL_MESO_DPD and defined using the GUI The properties for the species must then be set by clicking on set species see below for more details The parameters for pair potentials thermostats barostats electrostatics and surfaces can only be set once those for the species are defined The potential key gives the definition of potential energy functions used in calculations the standard DPD model by Groot and Warren 13 is the default but many body density dependent DPD 37 53 Lennard Jones 23 and Weeks Chandler Anderson hard sphere models can also be selected Note that the Lennard Jones and Weeks Chandler Anderson 57 models are non stochastic and not DPD models although the DPD thermostat can be used with them to maintain system temperature Parameters for unlike species interactions if more than one species is specified and for many body DPD interactions can be set by clicking on set potentials see below for more details The system thermostat the dissipative and random forces as defined for DPD with the standard molec ular dynamics form of the Velocity Verlet
195. particles to highlight the distinctions between the phases Cx Y Y il 4 L a 6 0 30 b 0 55 c 0 70 d 0 90 zi Figure 12 5 Visualizations of DPD AmphiphileMesophases test case at final time step isosurfaces of hydropho bic particles 12 6 Mixture_Par This parallel simulation example consists of 512 000 particles with 2 species Species 1 has a population of 256 000 and species 1 has a population of 256 000 particles The particle types have identical sizes and masses but different energy parameters using the default values for unlike particle parameters 12 7 VesicleFormation This parallel simulation example consists of 37 440 unbonded water particles and 1008 molecules each consisting of one hydrophilic head particle and three hydrophobic tail particles bonded together with stiff harmonic bonds 182 CHAPTER 12 DL MESO DPD EXAMPLES of equilibrium length 1 0 between them 58 The molecules represent amphiphiles and during the course of the simulation self assemble into a vesicle and encapsulate a number of water particles Figure 12 6 shows the self assembled vesicle both in three dimensions and in a cross section to show the encapsulated water An AVI video file of the simulation can be found in the subdirectory a Vesicle water omitted for clarity b Cross section orthogonal to x axis Figure 12 6 Visualizations of DPD VesicleFormation test case at t 43000 red for hydrophile green fo
196. perator C normally in the form of a matrix ABO fi Z t Ci and movement of these distribution functions to neighbouring lattice nodes Ji Z iAt t At fi z t which combine to give the governing equation for calculations 5 1 1 Collision algorithms The forms of collision currently available in DL MESO LBE are the Bhatnagar Gross Krook BGK single relaxation time 4 and Multiple Relaxation Time MRT schemes 5 1 1 1 BGK single relaxation time The BGK collision operator is defined by Ci fi t F 5 1 E where Ty is the relaxation time related to the kinetic shear viscosity of fluid by E 1 Az A 173 with the kinematic bulk viscosity of the fluid v set equal to 3v 6 This reduces the equation for post collisional distribution functions to fi zt f Et Fi T t A At Tf the same as Equation 3 12 Interfacial and external forces are dealt with either by adding uE to the velocity of the fluid 32 when calculating the equilibrium distribution function fj or by adding a forcing term to the collisional distribution function 14 1 Bie coma F 1 Wi A eij F 5 2 42 CHAPTER 5 DL MESO_LBE FEATURES where 7 is defined as equal to t ap and used in both the expression above and as the velocity for calculating equilibrium distribution functions The former is used by default in CollisionBGK while the latter method by Guo et al can be invoked usi
197. prop at each lattice point for system in Plot3D format e Dependencies fOutputQcA2D fOutputQCA3D e Arguments filename input array of characters e Comments The default output file name is 1bout q and the mass fraction of fluid 0 is output by default These can be changed by specifying an output file name and fluid number up to 1bsy nf 1 when calling the routine fOutputQCAIncom e Header records int fOutputQCAIncom const char filename l1bout int iprop 0 e Function Outputs mass fraction and speeds along x y and z directions of incompressible fluid iprop at each lattice point for system in Plot3D format e Dependencies fOutputQCA2DIncom fOutputQCA3DIncom e Arguments filename input array of characters e Comments The default output file name is 1bout q and the mass fraction of fluid 0 is output by default These can be changed by specifying an output file name and fluid number up to 1bsy nf 1 when calling the routine 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 69 fOutputQCB e Header records int fOutputQCB const char filename lbout int iprop 0 e Function Outputs concentration of solute iprop and compressible fluid velocity at each lattice point for system in Plot3D format e Dependencies fOutputQCB2D fOutputQCB3D e Arguments filename input array of characters e Comments The default output file name is 1bout q and the concentration of solute 0 is output by default These can
198. pter 7 DL MESO LBE Examples Test cases for Lattice Boltzmann Equation simulations using DL_MESO including the required input and sample output files can be found in the DEMO LBE subdirectory They can be run using either the serial or parallel versions of DL_MESO_LBE They consist of the following cases 7 1 2D Pressure This is a 2D simulation of a single fluid on a 42 x 42 grid with fixed pressure density boundary conditions on the left and right boundaries and solid walls represented by bounce back at the top and bottom A visualization with vector glyphs and a plot of fluid speed against vertical position can be seen in Figure 7 1 which show the boundary conditions result in a laminar flow with a parabolic velocity profile 0 005 3435 0 004 4 ae z t E 0 003 4 E 3 amp 3 2 E E E 0 002 4 T 0 001 4 Mop Ssh ef Se ee Slee fi i i 0 0 0001 0 0002 0 0003 0 0004 Vertical position y a Vector plot of system b Variation of fluid speed with vertical position Figure 7 1 Results from LBE 2D_Pressure test case 7 2 2D Shear This is a 2D simulation of a single fluid on a 42 x 42 grid with a shear boundary condition The vector plot in Figure 7 2 demonstrates the ability of the applied boundary conditions to generate a linear shearing Couette flow throughout the grid which is conf
199. put files for visualization These may be found in the LBE utility and DPD utility directories Compilation can either be carried out individually or collectively using makefiles each utility directory includes a makefile to compile all the utilities therein and the working directory WORK includes one to compile both sets for use with the GUI The latter can be invoked using the command e make f Makefile utils Some further details on these utilities can be found in the README files in the source directories B 1 DL_MESO LBE lbeplot3dgather lbeplot3dgather is a utility written in C to gather Plot3D output files produced by the parallel version of DL_MESO_LBE and produce a single structure file lbtout xyz or 1btout xy and a single set of solution files Ibtout q for visualization of the entire system If c is the command for the available C compiler the executable plot3d exe can be produced by typing e c o plot3d exe lbeplot3dgather cpp and either run at the command line plot3d exe or plot3d exe or via the GUI under Gather LBE Data All lbout xyz and lbout q files should be copied to the directory including the executable if necessary before running as well as the 1bout info file to give information on the sizes of integers and floating point numbers No user input is required although the utility will stop with an error message if no 1bout info file is available No other error messages are produced so care
200. quently used form for the weight function is 15 eN P p LA Ls w rij E 1 i rij lt ra 10 42 which reduces to standard DPD when the excess free energy per particle is set to 4 p 30 AP Multiple component many body DPD is also possible by defining partial local densities e g for component a p Y wr 10 43 jeaj i and generalizing Equation 10 40 pe a ae meee eta c i Tij 10 44 i Pes Pei rij with c i as the component to which particle belongs e g if i a c i a and 5 as the set of local densities of different components at the position of particle 7 The manybody_module includes a routine local_density to calculate local densities for each species using Equation 10 43 the overall density for each particle can be obtained by a simple sum over all species The user can modify the routines manybody_force and manybody_potential to apply their own choices for many body forces and potentials respectively Up to five many body interaction parameters per species pair can be specified in the CONTROL file The many body DPD example provided with DL_MESO_DPD produces a van der Waals like equation of state and can be used to model vapour liquid interfaces 55 The potential excess free energy per particle is given by 4 T TY ex 5 Aj p dBi 6 10 45 PN B 3 Aue Gy PaP 10 45 where P is equivalent to p but with the cutoff set to re instead of rg The associated pairwis
201. quid interactions used for many body forces 55 although only the density dependent potential energies are calculated using this routine the rest are calculated in manybody_force Users may wish to modify this routine to use their own many body interaction models 11 2 9 bond_module This module requires the modules constants variables error_module and comms_module to be loaded be forehand shellsort_list e Header records SUBROUTINE shellsort_list e Function Reorders the global local particle number list in terms of global particle number using a Shell sort e Dependencies None search_list e Header records INTEGER FUNCTION search_list aim 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 153 e Function Determines the index for the global local particle number list for a specified global particle number aim using a binary search e Dependencies None e Arguments aim input integer search_list output integer e Comments This function returns a negative value if the global particle number cannot be found in the list If it is not the only entry for the global particle number the function returns the index plus the value of nlist number of list items to flag up duplicate entries duplicate e Header records SUBROUTINE duplicate global1 global2 indi ind2 e Function Determines the indices out of duplicate entries from the global local particle number list for a pair of specified global par
202. r e Comments This bounce back boundary condition is carried out using f lbopv i f 1bv i fBounceBackT e Header records int fBounceBackC long tpos e Function Performs an on grid bounce back for the temperature distribution function at the tpos th grid point 78 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE e Dependencies None e Arguments tpos input long integer e Comments This bounce back boundary condition is carried out using f lbopv i 1bv i fMidBounceBackF e Header records int fMidBounceBackF long tpos e Function Performs a mid link bounce back for the fluid distribution function at the tpos th grid point e Dependencies None e Arguments tpos input long integer e Comments This bounce back boundary condition is carried out by exchanging post collisional populations with con jugate values in neighbouring grid points fMidBounceBackC e Header records int fMidBounceBackC long tpos e Function Performs a mid link bounce back for the solute distribution function at the tpos th grid point e Dependencies None e Arguments tpos input long integer e Comments This bounce back boundary condition is carried out by exchanging post collisional populations with con jugate values in neighbouring grid points fMidBounceBackT e Header records int fMidBounceBackC long tpos e Function Performs a mid link bounce back for the temperature distribution functi
203. r hydrophobe blue for water Appendix A Manual compiliation and running of DL_MESO A 1 DL MESO LBE DL_MESO LBE has been written in C in a modular form where the main codes for serial and parallel running slbe cpp and plbe cpp automatically include all the required modules To compile the code and produce an executable in the working directory at a command line either type e c LBE slbe cpp o lbe exe to compile the serial version or e c LBE plbe cpp o lbe exe to compile the parallel version assuming that c is the name of the available C compiler and lbe exe is the name of the executable required After creating or copying the required input files lbin sys and 1bin spa into the working directory DL_MESO_LBE can then either be run at the command line lbe exe or 1be exe depending on the operating system or launched using a run script By default DL MESO_LBE models compressible fluids using the Bhatnagar Gross Krook single relaxation time scheme and applies a correction of uE to the velocity for interfacial and body forces To model incompressible fluids calls to the following routines in slbe cpp or plbe cpp should have Incom appended to the end fInitializeSystem fReadInitialState e fOutput or fsOutput e fCollision fPostCollBoundary fPostPropBoundary The Guo forcing scheme can be applied by appending Guo to the call for the collision scheme e g Co11isionBGKGuo To use mu
204. r than the first but more readable fGetSpeedSite e Header records three cases int fGetSpeedSite double speed double startpos int fGetSpeedSite double speed long tpos int fGetSpeedSite double speed int xpos int ypos int zpos e Function Calculates the macroscopic speed of all compressible fluids at the grid point e Dependencies None e Arguments startpos input double pointer tpos input long integer xpos input integer ypos input integer zpos input integer speed output double precision array e Comments The calculation is based on va tes ista The second and third cases are more readable than the first but also slightly slower fGetSpeedIncomSite e Header records three cases int fGetSpeedIncomSite double speed double startpos int fGetSpeedIncomSite double speed long tpos int fGetSpeedIncomSite double speed int xpos int ypos int zpos e Function Calculates the macroscopic speed of all incompressible fluids at the grid point e Dependencies None 62 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE e Arguments startpos input double pointer tpos input long integer xpos input integer ypos input integer Zpos input integer speed output double precision array e Comments pee Filio ie The second and third cases are more readable than the first The calculation is based on vg but also slightly slower fGetOneMomentSite e Header records three cases int fGetOneMomentSite
205. raction distance between components i and j 10 5 LONG RANGED ELECTROSTATIC COULOMBIC POTENTIALS 133 of each other forming infinitely strong ion pairs The charges are therefore spread out over a finite volume using a smearing charge distribution f r The current approach uses a Slater type distribution i e f r a exp 10 50 where A is the decay length of the charge This gives a potential energy between charged particles and j of Pag U rij 4 1 1 Brij exp 28r55 10 51 TYi5 and the corresponding electrostatic force gt TPg g Pi Fj rig 3 1 exp 28 rig 1 28r 1 Briz 10 52 Arr F Tij Te where re is the electrostatic short range real space cutoff and 6 For large particle separations these expressions reduce down to the standard expressions for point charges and thus the reciprocal space part of the Ewald sum can be calculated without modification The real space terms for the Ewald sum are modifications of the above expressions evaluated for particle pairs when rj lt re The short range potential energy between particles 7 and j is given as sr 495 Uj ie 1 1 Brij exp 26ri 10 53 and the pairwise force Be ENG a PO extend US expt Sabre ere Bre 10 54 a Anri VT 2 J J J A Tij where a is a convergence parameter that controls the real space contribution chosen to give negligible contribu tions beyond the real space cutoff If
206. rameters per grid point lbsitelength int number of grid points in yz plane lbyz int grid spacing 1bdx double time step 1bdt double speed of sound 1bsoundv double Reynolds number lbreynolds double processor name 1bdm rank int total number of processors lbdm size int x coordinate for domain lbdm xcor int y coordinate for domain lbdm ycor int z coordinate for domain lbdm zcor int number of processors along x axis lbdm xdim int number of processors along y axis lbdm ydim int number of processors along z axis lbdm zdim int x coordinate of domain start position 1bdm xs int x coordinate of domain end position 1bdm xe int y coordinate of domain start position lbdm ys int y coordinate of domain end position lbdm ye int z coordinate of domain start position lbdm zs int z coordinate of domain end position 1bdm ze int inner grid points along x lbdm xinner int outer grid points along x lbdm xouter int inner grid points along y lbdm yinner int outer grid points along y lbdm youter int inner grid points along z lbdm zinner int outer grid points along z lbdm zouter int total grid points lbdm touter int name of neighbouring processor lbnb rank int position for sending distribution message position for receiving distribution message lbnb spos lbnb rpos unsigned long unsigned long 2Excluding boundary points 3Including boundary points 4 3 THE PARAMETERS AND THEIR FUNCTIONS Table 4 7
207. rature 298K with Nm molecules per DPD particle T 20 00Nm Alternatively the total electrostatic potential energy can be expressed as a sum of Coulombic energies which also include periodic images i e T qidj 10 48 An De IRAE razl i J gt 1 One method to determine electrostatic interactions is currently included in DL_MESO_DPD This can be spec ified in the CONTROL file using the keyword ELEC etype is an integer specifying the electrostatic calculation method and a b c and d are parameters whose meanings are given in Table 10 3 Table 10 3 Electrostatic coefficients used in DL MESO_DPD etype a b cld 1 Standard Ewald sum exponential smearing DP a 8 10 5 1 Standard Ewald sum with exponential charge smearing etype 1 The method currently used in DL_MESO_DPD to determine the electrostatic potential is an Ewald summation 11 U aw U u NT 10 49 where U is a short range potential energy term that sums quickly in real space U is a long range term that sums quickly in Fourier or reciprocal space U is a self energy correction term and U is a correction for systems with a net charge While the original form of the short range electrostatic potential uses point charges this cannot be used un modified for DPD simulations soft beads used in combination with unlike point charges would collapse on top 3 Bij in this case should not be confused with the inte
208. re If the GUI is used this can be selected using the pulldown list in the Gather LBE Data panel If the argument is omitted the utility will ask the user to enter the required property No other user input is required but error messages will be produced if either of the files lbout info and lbout ext are missing No other error messages are produced so care should be taken to ensure no VTK files for the pieces are missing particularly since these files are required for plotting as the linking files do not include the data B 2 DLMESO DPD molecule generate molecule generate is a utility written in C to generate the input files required for modelling particles in DPD simulations that are bonded together i e molecules A random flight generation system is used to generate the coordinates of bonded beads which can form branched molecule chains a constant distance apart within a cube of a size specified by the user which will be used by DL_MESO_DPD to insert the molecule into the system This utility can be compiled to produce the executable molecule exe with the command e c o molecule exe molecule generate cpp if c is the command for the available C compiler This utility can be run from the command line or via the GUI in Set DPD Molecules which runs the utility in a new command line shell window When running the utility the user will be asked for the number of molecules required If a CONTROL file exists in the same d
209. re or after the call for creating output files e Addition of fsForcePeriodic inside main loop after calculating interaction forces if these are not to be calculated this call can be omitted 5 2 Boundary conditions To apply boundary conditions to a Lattice Boltzmann Equation simulation the distribution functions f at boundary lattice points have to be modified or replaced during each time step to give the required fluid velocity or pressure concentration temperature This may take place either between the collision and propagation stages or at the end of each time step the subroutines fPostCollBoundary and fPostPropBoundary respectively are used to invoke the boundary conditions The space property lbphi is used to define the boundary conditions for each lattice node in the system 5 2 1 Periodic Periodic boundaries are used to simulate bulk fluids sufficiently far away from the actual boundaries of a real physical system so that surface effects can be neglected As the fluid moves out of one face of the system volume it reappears on the opposite face with the same velocity density etc DL_MESO_LBE applies periodic boundary conditions in two different ways depending on the size of the bound ary domain 1bdm bwid If there is no boundary domain the default for serial running periodic boundary conditions are automatically applied during the propagation step by using the function fCppMod to adjust the destination of distribution function
210. red LBE model can be selected from the pull down list the D2Q9 D3Q15 D3Q19 and D3Q27 square lattice schemes are available 2 The number of grid points sets the size of the system For 2D systems the number of grid points in the z direction must equal 1 selecting a two dimensional lattice model greys out this box 3 The total steps for the simulation and the save span number of timesteps between system outputs are given in this row 2 2 LATTICE BOLTZMANN AND THE DL_MESO GUI 7 10 11 12 13 14 L MESO LBE DPD SPH Manual Help Exit Lattice Boltzmann LBE madal p2a9 y 1 number of grid points x 40 y 30 2 2 total steps 10 save span 5 3 boundary width A 4 sound speed 540 5 kinetic viscosity 0 001 6 top boundary speed vx 0 vy 0 Vz T Change LBE Code bottom boundary speed vx 0 Vy 0 Vz 8 left boundary speed vx 0 Vy 0 Vz 9 Compile LBE Code right boundary speed vx 0 vy 0 vz 10 front boundary speed Vx Vy Vz 11 Run LBE Program back boundary speed Vx Vy vz 12 noise 0 1 o phase field Y no phase fidd Gather LBE Data number of fluids phases 1 gt gt set fluid parameters 14 number of solutes o gt gt set solute parameters15 Piot LBE Results using temperature scalar O yes gt gt set thermal parameters16 Figure 2 2 Define LBE System
211. requires the calculation of the wall velocity which can be simplified by setting non orthogonal velocity components to zero For the analogous example at the top wall for D2Q9 P CEDF the same equations for f4 f3 and fs can be used together with PwUVw x 0 PwUw y fo fo fe 2 fi fr fs Pw denotes any valid letters for the solute composition and temperature 1 9 gt 46 CHAPTER 5 DL MESO_LBE FEATURES One complication for three dimensional lattices is the requirement to apply the non equilibrium bounce back to all unknown distribution functions which ordinarily overspecifies the system but can be counteracted using transverse momentum corrections for directions other than orthogonal to the boundary which are non zero for e g shearing flows It should be noted that if the wall velocity is set to zero the boundary condition reduces to on grid bounce back DL_MESO_LBE includes implementations of the Zou He boundary condition for all four lattice schemes D2Q9 D3Q15 D3Q19 and D3Q27 Only planar surfaces can be dealt with using this method concave edges in three dimensions and concave corners instead use the equilibrium distribution function for the given density velocity and either zero velocity or the density at the nearest fluid grid point 5 2 5 Constant solute concentration temperature To specify constant solute concentrations or temperatures at planar boundaries the Inamuro method 21 20 is used in DL_MESO
212. rgy stkbe real KIND dp A mxstak data stack for angle potential energy stkae real KIND dp A mxstak data stack for dihedral potential energy stkde real KIND dp A mxstak data stack for electrostatic potential energy stkee real KIND dp A mxstak duni random number generator state uni integer A 102 mtrnd random number generator state mt integer A 0 624 Chapter 10 DL MESO DPD Features 10 1 Domain decomposition and linked list cell calculations The Domain Decomposition DD strategy is one of several ways to parallelize particle based simulations 45 Its basis is the division of the simulated system into equal sized spatial blocks or domains each of which is allocated to a specific processing unit of a parallel computer The arrays defining the coordinates velocities and forces for all N particles in the system are divided into sub arrays of size x on each of the P processing units with the particles allocated geometrically among them In order for the strategy to work efficiently the simulated system should possess a reasonably uniform density so that each processing unit is allocated as equal a portion of particle data as possible The computation of forces and integration of the equations of motion are shared more or less equally between the processing units and to a large extent can be computed independently on each unit While tricky to program this method is conceptually simple and particularly suited to large scale simulati
213. ro the lines can be omitted Each particle is represented by a block record with at least two lines of information e The species name 8 characters or number integer and the global particle number integer optional e The x y and z coordinates for the particle real The x y and z components of particle velocity real if levcfg gt 0 e The x y and z components of force on the particle real if levcfg gt 1 If global particle numbers are not included these are generated automatically for the particles in the order specified by the CONFIG file Care should be taken that any particles belonging to molecules are numbered correctly since the bond information is assigned in an identical fashion to unspecified systems i e numbering after all loose particles in the relative order specified by the MOLECULE and FIELD files CONFIG files can be created from restart files of previous simulations using the supplied Fortran90 program exportconfig in the directory DPD utility see Appendix B for more details 11 3 2 Output files General output file OUTPUT This ANSI text file is generated by all DPD calculations and contains e The system and bond angle dihedral properties used for calculations e Domain decomposition details Parallel version only e The starting positions and velocities of a particle sample e The calculation time current values and rolling averages for the total energy potential energy total elec tro
214. rom LBE 2D_KarmanVortex test case scale blue to red 7 5 2D_LidCavity This is a 2D simulation of a single fluid on a 128 x 128 grid with a shear boundary condition at the top and solid walls surrounding the other edges of the system resulting in lid driven cavity flow Figure 7 5 shows the fully developed velocity field for a Reynolds number of 100 at the final time step 7 6 83D PHASESEPARATION 107 velocity X 0 16 0 1 0 08 Figure 7 5 Magnitude plot of x component velocity and velocity vector plot from LBE 2D_Lidcavity test case 7 6 3D _PhaseSeparation This is a 3D simulation of a binary fluid on a 100 x 100 x 100 grid with a fixed boundary condition and non zero mesoscopic interaction parameters producing phase separation of the two fluids Figure 7 6 shows the phase separation process in a number of snapshots two AVI video files have been rendered from an example calculation one giving a 3D view of the system the other showing a plane normal to the z axis which can be found in the directory for this case L a a t 400 b t 1200 c t 2000 Figure 7 6 Progressive density plots in plane normal to z axis from LBE 3D_PhaseSeparation test case red for fluid 0 blue for fluid 1 108 CHAPTER 7 DL_MESO LBE EXAMPLES 7 7 3D_Shear This is a 3D simulation of a single fluid on a 40 x 30 x 25 grid with a shear boundary condition Figure 7 7 shows a vector plot for this system demonstrating that linear shear
215. rs can be set by clicking on set barostat and the target system pressure can be specified 17 The electrostatics scheme for the simulation an Ewald sum method with Slater type exponential charge smearing 11 is available in DL _MESO_DPD If selected the short range electrostatic cutoff can be edited and the parameters for the Ewald sum can be specified by clicking on set electrostatics 18 The surfaces to be applied to the system by default periodic boundary conditions are used but hard walls can be applied at the boundaries with soft repulsions and specular reflection 56 Selecting this method allows the surface cutoff to be edited and the wall parameters can be specified by clicking on set surfaces Note that the DPD code uses reduced units in which the unit of length is the particle size the unit of mass is the particle mass and the unit of energy is the primary energy parameter of the potential energy function From these the time unit may be derived The temperature is defined to be 2 of the system kinetic energy The CONTROL file for input into DL_MESO_DPD is created by clicking the SAVE button 2 3 1 1 Species potential thermostat barostat electrostatic and surface parameters Examples of the pop up windows for species potential thermostat barostat electrostatic and surface parame ters can be seen in Figure 2 6 with the rows labelled in red numbering multiple columns of dialogue boxes are made available for systems wit
216. rtpos double fGetFracSite int fpos long tpos double fGetFracSite int fpos int xpos int ypos int zpos e Function Calculates the mass fraction of fluid fpos in the site e Dependencies None e Arguments startpos input double pointer fpos input integer tpos input long integer xpos input integer ypos input integer zpos input integer fGetFracSite output double precision e Comments The calculation is based on z 4 This function operates in a similar way to fGetOneMassSite with i the second and third cases slightly slower than the first but more readable fGetOneConcSite e Header records two cases double fGetOneConcSite int cpos long tpos double fGetOneConcSite int cpos int xpos int ypos int zpos e Function Calculates the concentration of solute cpos at the grid point e Dependencies fGetOneMassSite e Arguments tpos input long integer cpos input integer xpos input integer ypos input integer Zpos input integer fGetFracSite output double precision fGetTemperatureSite e Header records two cases double fGetTemperatureSite long tpos double fGetTemperatureSite long xpos long ypos long zpos e Function Calculates the scalar temperature at the grid point 60 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE e Dependencies fGetOneMassSite e Arguments tpos input long integer xpos input long integer ypos input long integer Zpos input long integer fGetFracSite output do
217. s filename input array of characters Comments The default output file name is lbout vts and the mass fraction of fluid 0 is output by default These can be changed by specifying an output file name and fluid number up to 1bsy nf 1 when calling the routine 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 75 fOutputVTKCB e Header records int fOutputVTKCB const char filename l1bout int iprop 0 e Function Outputs concentration of solute iprop and compressible fluid velocity at each lattice point for system in Structured Grid XML VTK format e Dependencies fOutputVTKCB2D fOutputVTKCB3D e Arguments filename input array of characters e Comments The default output file name is 1bout vts and the concentration of solute 0 is output by default These can be changed by specifying an output file name and solute number up to 1bsy ns 1 when calling the routine fOutputVTKCBIncom e Header records int fOutputVTKCBIncom const char filename l1bout int iprop 0 e Function Outputs concentration of solute iprop and incompressible fluid velocity at each lattice point for system in Structured Grid XML VTK format e Dependencies fOutputVTKCB2DIncom fOutputVTKCB3DIncom e Arguments filename input array of characters e Comments The default output file name is 1bout vts and the concentration of solute 0 is output by default These can be changed by specifying an output file name and solute number up to 1bsy ns
218. s meq output double pointer p input double pointer rho input double precision rho0 input double precision e Comments The equilibrium distribution function in moment space calculated here is Mea Tf the exact form of which is dependent on the lattice scheme given for D2Q9 by 26 and for D3Q15 and D3Q19 by 7 Parameters for calculating the square of energy e and fourth order moments Taa can be modified by the user in the 1bpMODEL module fGetMomentForce e Header records int fGetMomentForce double source double v double force e Function Calculates Guo like forcing terms in moment space e Dependencies None e Arguments meq output double pointer v input double pointer force input double pointer 90 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE e Comments The forcing terms in moment space 39 calculated here are T 3 T 9 0 F the exact form of which is dependent on the lattice scheme given for D2Q9 D3Q15 and D3Q19 fGetMRTCollide e Header records int fGetMRTCollide double collide double omegashear double omegabulk e Function Calculates collision vector for Multiple Relaxation Time MRT scheme with specified fluid relaxation frequencies e Dependencies None e Arguments collide output double pointer omegashear input double precision omegabulk input double precision e Comments The exact form of the collision vector 5 is dependent on the lattice scheme
219. s and positions of the particles are calculated by integration of slightly modified differential equa tions dVi a Esa 1 dt 7 Ug aVi a Np e 2 Uga 10 19 d 1 0 a Via Ug aise 10 20 where the force on particle 2 F includes any thermostating forces the time integral of these forces can be determined for all thermostat types The implementation of this barostat is carried out using the Velocity Verlet scheme to integrate the equations of motion for both the particles and the barostat The first Velocity Verlet stage integrates the forces on the particles At Fa t Vi a t At Vio t jeak Atug aUi a 10 21 which is followed by a similar integration for the barostat velocity At Fy t Ug a t At uga t gt wW 10 22 before the positions and simulation box dimensions are updated Tia t At exp uga t 4At At ria t Atvi a t 5At 10 23 La t At exp uga t 4At At La t 10 24 At this point the forces at the end of the time step are calculated including thermostating forces along with the system pressure Since the the barostat force requires correct velocities for both the particles and barostat an iterative procedure is required which begins by calculating an initial guess for the barostat velocity at the end of the time step 2F a t At g u t At uga t At 4 g 0 Each iteration starts by calculating the particle velocities in a slightly modified sec
220. s between the two particles The velocities of particles and 7 thus become Ui Vi Bis es Tiz vg ij 10 5 Mi gt gt Hij gt gt A j j E ig ig Vy i 10 6 l The probability of a particle pair being thermostatted is equal to rAt where I is defined as the collision frequency with a maximum value of 25 and the velocity corrections to particle pairs are applied in a random order to prevent biasing The viscosity and self diffusion generated by this thermostat for a single species are norr 10 7 75m ane kgT p BBE PD 10 8 m where Tp is the decay time for velocity correlations and inversely proportional to the collision frequency The Schmidt number is therefore proportional to Er and can thus reach values up to O 10 This thermostat is suited to systems with higher viscosities and low diffusitivies while giving the correct system temperature for a wide range of time step sizes within numerical errors due to Velocity Verlet force integration Its implementation in parallel running uses a replicated data strategy to carry out the velocity corrections this requires additional memory on each processing unit for the velocities of all particles and the data required to modify the velocities of particle pairs The efficiency of the Lowe Andersen thermostat thus decreases with increasing numbers of particles in the entire system and I 126 CHAPTER 10 DL MESO_DPD FEATURES 10 2 4
221. s leaving the system so they are placed at the opposite side No periodic boundary using this implementation needs to be defined by the space property which can be left equal to zero as for the bulk fluid 5 2 BOUNDARY CONDITIONS 45 For systems that include a non zero boundary domain size the distribution functions at the edges of the actual system are copied into the buffer at the opposing sides before collisions and propagation take place This requires the use of the space property lbphi i 10 to determine the location of the domain buffer which can be set up using the routine fMarkBoundArea and either fsPeriodic and similar routines for serial running or fNonBlockCommunication and similar routines for running in parallel to copy the distribution functions into the buffers 5 2 2 On grid bounce back The on grid bounce back condition applies a no slip condition i e zero fluid velocity at a boundary that lies halfway between grid points This is applied after the propagation stage by reversing the distribution functions sitting on each wall node Zw i e Fi w t fj Zw t 5 11 where j is the opposite lattice link to i i e The reflection of distribution functions occurs on grid On grid bounce back is a first order approximation of the boundary condition 36 i e the error is proportional to the lattice spacing Az but it is completely local i e only uses distribution functions at the wall node 5 2 3 M
222. s total and individual fluid masses in entire system e Dependencies fGetTotMassDomain fGetOneMassDomain fGlobalValue fPrintSystemMomentum e Header records int fPrintSystemMomentum e Function Calculates and prints the total fluid momentum in entire system e Dependencies fGetTotMomentDomain fGlobalValue 6 3 DL MESO LBE Data Files 6 3 1 Input files Define system lbin sys The use of the DL_MESO GUI is recommended for producing lbin sys although existing files of that name can also be edited Superfluous parameters can be deleted and new parameters introduced although the latter requires an addition to the parameter recognition loop in the fInputParameters subroutine Define space lbin spa The GUI is recommended for creating the 1bin spa file which stores the data in the following format 102 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE X y Z grid property An empty 1bin spa file represents all boundaries as periodic Define initial condition 1bin init This optional file cannot currently be created by the GUI the user must create this file if it is required The following format is required for each lattice point whose default velocity fluid densities solute concentrations or temperature needs replacing X y Z Ux gt Uy Uz fPo Plbsy nf 1 CO Clbsy nc 1 ST Note that solute concentrations and temperature can be omitted if these are not required in LBE simulations 6 3 2 Output files
223. scale that is too long to be modelled effectively by MD DPD may also be used when such systems experience shear and flow gradients The DPD algorithm can be summarized by the following A condensed phase system may be modelled as a system of free particles interacting directly through soft forces e The system is coupled to a heat bath via stochastic forces which act on the particles in a pairwise manner e The particles also experience a damping or drag force which also acts in a pairwise manner e Thermodynamic equilibrium is maintained through the balance of the stochastic and drag forces i e the method satisfies the fluctuation dissipation theorem e At equilibrium or steady state the properties of the system are calculated as averages over the individual particles as in Molecular Dynamics 8 2 Outline of Method In DPD the system is modelled as a system of free particles which are spherical and interact over a range that is of the same order as their diameters The particles can be thought of as assemblies or aggregates of molecules such as solvent molecules or polymers or more simply as carriers of momentum lThe outline of the DPD method supplied here is based on 13 111 112 CHAPTER 8 DISSIPATIVE PARTICLE DYNAMICS BASIC THEORY The equations governing the time evolution in a DPD simulation resemble those of ordinary MD di F So es 1 dr ey 8 2 T v 8 2 in which r v and F are the po
224. should be taken to ensure no solution files are missing Since the data is copied into the combined structure and solution files the original output files can be deleted after this utility is run lbevtkgather lbevtkgather is a utility written in C to gather Structured Grid XML formatted VTK output files produced by the parallel version of DL_MESO_LBE lbout vts and produce a set of linking files lbtout pvts for lif using the GUI and the utilties are to be compiled manually or in their source directories copies of the executables are required in the directory from which the GUl is to be launched e g WORK 187 188 APPENDIX B DL_MESO UTILITIES visualization of the entire system If c is the command for the available C compiler the executable vtk exe can be produced by typing e c o vtk exe lbevtkgather cpp and either run at the command line or via the GUI under Gather LBE Data All 1bout vts files should be copied to the directory including the executable if necessary before running as well as the lbout info and lbout ext files to give information about the number of processors used for the simulation and the extents of each piece This utility can be run with a command line argument to give the scalar property required e g for Windows and Unix Linux computers respectively e vtk exe 1 e vtk exe 1 where 0 is used for fluid density 1 for mass fraction 2 for solute concentration and 3 for temperatu
225. sition velocity and force of the ith particle which has mass m The force on the particle is a sum of pair forces N R PC BD BR K ES a t FR 8 3 IA in which ES represents the force exerted on particle due to the presence of particle 7 Additional pairwise forces may be BD and ES are the conservative drag and random or stochastic pair forces respectively Each included for more complicated systems such as those involving chains of particles bonded together 40 The conservative interactions are usually soft i e weakly interacting so that the particles can pass by each other or even through each other relatively easily so that equilibrium is achieved quickly A common form of interaction potential is an inverse parabola 2 1 Tij where ri r jl re is a cutoff radius and Aj is the interaction strength A may be the same for all particle pairs or may be different for different particle types Equation 8 4 gives rise to a repulsive force of the form RO Cip Tij Tij Tig Fg Aig rij Aij 1 man 8 5 ij c ij This is the deterministic or conservative force Fe exerted on particle i by particle j Note the switching function w rij and the force are zero when rj gt re and thus the particles have an effective diameter of 1 in units of the cutoff radius re The stochastic forces experienced by the particles is again pairwise in nature and takes the form R 1 ij OijWw rij C
226. ssipation this parameter is equal either to the dissipative factor y as used by the DPD and Peters thermostats or to the collision frequency T for the Lowe Andersen and Stoyanov Groot thermostats If this value is not specified or set to zero the value for the entire system will apply 8 Body force x axis x component of constant external body force on each particle 9 Body force y axis y component of constant external body force on each particle 10 Body force z axis z component of constant external body force on each particle Potentials require the following data see Figure 2 6 b 1 The energy Aij distance Bij and dissipation Dij parameters between species and j If any of these values are set to zero DL_MESO_DPD will use default values i e Ajj y A Aj Bij 3 Bi Bj Vij y Clicking add unlike will store the entered parameters to memory and add a record to the text box which cannot be edited by the user reorder unlike will reorder the specified interactions and reset unlike will clear all of them 2 The many body DPD parameters Aij Bij Cij Dij and Eij between species i and j Parameters for all possible species pairs must be specified by the user Clicking add mb stores the entered parameters to memory and add a record to the text box reorder mb reorders the specified interactions and reset mb will clear all many body DPD parameters The thermostat pop up window is formatted as in Figure 2 6 c
227. static and from bond stretching angles and dihedrals virial kinetic energy pressure and temperature every NSBPO time steps e Final averages and fluctuations standard deviations over all time steps after equilibration e The final positions and velocities of a particle sample e Elapsed and average times for the calculation 11 3 DL_MESO_DPD DATA FILES 177 Restart file s export Each processing unit produces a restart file with a name beginning with export every ndump time steps This binary file contains the following information for the time step e Name of DPD calculation e Numbers of particles bonds angles and dihedrals in the processing unit e Specified temperature number of time steps system volume e Current state of random number generators e Particle global identity numbers species and molecule numbers e Bond angle and dihedral tables e Particle Cartesian coordinates velocities forces and virials e Current time step and statistical properties current values cumulative sums and fluctuations rolling averages and stacks In combination the export files can restore a stopped DPD calculation They can also be used to generate a plot at the given time step using the utility exportimage in the DPD utility directory Appendix B gives instructions for its use Trajectory file s HISTORY If the parameter 1sbsv is set to true each processing unit will generate a trajectory file with a name beginni
228. system viscosities and diffusivities with good temperature control and hydrodynamics using the collision frequency to obtain the required Schmidt number The replicated 10 3 BAROSTATS 127 data strategy is again used for the Lowe Andersen part which requires memory in each processing unit to store the velocities of all particles and the data for particle pair modification using the Lowe Andersen scheme the Nos Hoover scheme calculates the thermostating forces locally and does not require additional memory 10 3 Barostats In addition to a thermostat a barostat may be included in simulations to obtain a desired average pressure Po by adjusting the size and shape of the simulation cell DL_MESO_DPD includes two such algorithms a Langevin type barostat 22 and the Berendsen barostat 2 both of which have been coupled to all five available thermostats The instantaneous isotropic pressure in a system is calculated using the virial theorem 1 Pt ap m0 DAO 0 10 15 while for anisotropic orthorhombic systems the pressure in dimension a related to the instantaneous stress tensor component Caa t is defined as 1 Pot VO X Miti alt D Falat 10 16 All barostat definitions are expressed for the more general anisotropic case these can be applied for isotropic systems by setting P t P t and P t all equal to P t The flag lisoprs can be used to determine whether the system pressure is isotropic
229. t xyz for 3D systems and lbout xy for 2D systems fOutputQ e Header records int fOutputQ const char filename l1bout int iprop 0 e Function Outputs macroscopic mass density and velocity speeds along z y and z directions at each lattice point of compressible fluid iprop for system in Plot3D format e Dependencies fOutputQ2D fOutputQ3D e Arguments filename input array of characters e Comments The default output file name is 1bout q and the density of fluid 0 is output by default This can be changed by specifying an output file name and fluid number up to 1bsy nf 1 when calling the routine fOutputQIncom e Header records int fOutputQIncom const char filename lbout int iprop 0 e Function Outputs macroscopic mass density and velocity speeds along x y and z directions at each lattice point of incompressible fluid iprop for system in Plot3D format 68 CHAPTER 6 DL_MESO_LBE PACKAGE REFERENCE e Dependencies fOutputQ2DIncom fOutputQ3DIncom e Arguments filename input array of characters e Comments The default output file name is 1bout q and the density of fluid 0 is output by default This can be changed by specifying an output file name and fluid number up to 1bsy nf 1 when calling the routine fOutputQCA e Header records int fOutputQCA const char filename lbout int iprop 0 e Function Outputs mass fraction and speeds along x y and z directions of compressible fluid i
230. t output integer KIND li 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 149 getdble e Header records REAL KIND dp FUNCTION getdble txt n e Function Reads the nth word of the string txt to obtain a double precision number e Dependencies parsedble e Arguments txt input character 160 n input integer getdble output real KIND dp 11 2 6 surface module This module requires the variables module to be loaded beforehand surfacenodes e Header records SUBROUTINE surfacenodes e Function Identifies nodes containing surfaces or other boundary conditions e Dependencies None hardreflect e Header records SUBROUTINE hardreflect k e Function Applies boundary condition for hard reflecting walls calculates short range forces on particles when k 1 applies bounce back condition when k 2 for particles about to pass through boundary e Dependencies surfacebounce e Arguments k input integer e Comments The boundary condition is given by 56 the applied wall potential on each particle is given as 1 2 Uwaiti Z zwali 2 for z lt Ze 150 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE surfacebounce e Header records SUBROUTINE surfacebounce ddd vdd sided e Function Applies boundary condition on a leaving particle by means of specular reflection move particle back into system and invert velocity component normal to wall e Dependencies None e Arguments ddd input o
231. t vector is p y e T M 1 p 61500 Jy ly Pees Pey with p as the density e the energy e the square of energy j momentum q energy flux Pax the diagonal stress tensor component and pz the off diagonal stress tensor component The transformation matrix is 1 1 1 1 1 1 T chy 2 4 2 1 2 1 2 1 2 I 4 1 2 1 2 1 2 1 2 0 1 1 1 0 1 1 1 0 T 0 1 2 1 0 1 2 1 0 0 1 0 1 1 1 0 1 1 0 1 0 1 2 1 0 1 2 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1l 0 1 O The equilibrium moment vector is de Jaz gt A 5 ig E Ht tee p p 3p 3p with we and wej as adjustable parameters for convergence to the single relaxation time BGK scheme these are set equal to 1 and 3 respectively Of the 9 collision parameters available sy s3 and s5 have no effect except when applying external forces when they should be set equal to one as the associated moments are preserved and sa s4 and sg are tuneable parameters for calculational stability with the condition that s4 sg The remaining parameters represent the viscosities of the fluid i e 1 1 SEN AE ANA 1 1 Ar Tr Sep 2 De als 2 AE aL 2 A PE 2 Ax 1 Esmas 5 Ax 6 s 2 At 6 2 At Le Te x and Tf bulk A If the moment transformed source terms are to be used the vector S for this lattice scheme is defined as 0 6 uz Fs uy Fy 6 uz Fr uy Fy Fo F Fy f 2 Ua Fy UyFy Ug Fy Uy Py Y II Similar schemes are
232. tatic potential energy aveee real KIND dp system potential energy fluctuation flcpe real KIND dp system virial fluctuation flevir real KIND dp system kinetic energy fluctuation flcke real KIND dp system total energy fluctuation flcte real KIND dp system pressure fluctuation flcprs real KIND dp system volume fluctuation flcvlm real KIND dp system temperature fluctuation flcttp real KIND dp system bond potential energy fluctuation flcbe real KIND dp system angle potential energy fluctuation flcae real KIND dp system dihedral potential energy fluctuation flcde real KIND dp system electrostatic potential energy fluctuation flcee real KIND dp system potential energy accumulator zumpe real KIND dp system virial accumulator zumvir real KIND dp system kinetic energy accumulator zumtke real KIND dp system volume accumulator zumvlm real KIND dp system bond potential energy accumulator zumbe real KIND dp system angle potential energy accumulator zumae real KIND dp system dihedral potential energy accumulator zumde real KIND dp system electrostatic potential energy accumulator zumee real KIND dp system potential energy at current step stppe real KIND dp system virial at current step stpvir real KIND dp system kinetic energy at current step stptke real KIND dp system total energy at current step stptke real KIND dp system pressure at current step stpprs real KIND dp system volume at current step stpvlm real KIND dp
233. tcomms e Function Starts Message Passing Interface MPI e Dependencies None exitcomms e Header records SUBROUTINE exitcomms e Function Closes Message Passing Interface MPI in a controlled manner e Dependencies None abortcomms e Header records SUBROUTINE abortcomms e Function Terminates Message Passing Interface MPI e Dependencies None 144 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE gsync e Header records SUBROUTINE gsync e Function Pauses running until all processes are synchronized e Dependencies None global_and e Header records SUBROUTINE global_and iii nnn e Function Finds global logical AND from a Boolean array iii0 of size nnn e Dependencies None e Arguments iii input output array of logical nnn input integer e Comments For Boolean scalars the alternative SUBROUTINE global_sca_and iii is available global_sum_dble e Header records SUBROUTINE global_sum_dble aaa nnn e Function Globally sums double precision array aaa of size nnn e Dependencies None e Arguments aaa input output array of real KIND dp non input integer e Comments For double precision scalars the alternative SUBROUTINE global_sca_sum_dble aaa is available global_sum_int e Header records SUBROUTINE global_sum_int iii nnn e Function Globally sums integer array iii of size nnn e Dependencies None 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTION
234. te is given as At gi 7 e At t At g 1 t 9 t 9 2 5 18 Ts where g is the distribution function for the solute and 7 the solute relaxation time which is related to its ad 1 Az D 5 7 5 At diffusivity and the Schmidt number can be determined by Vy 277 1 a a Taking the concentration of the solute as C J gi the equilibrium distribution function for the solute is given by a simpler form of Equation 3 11 g 1 wiCs i Al 5 19 c where the velocity used is that of the bulk fluid Heat transfers can be coupled to the system in a similar manner using a thermal distribution function h and a thermal relaxation time 7 which gives the governing equation At hi Z At t At hi Z t hi Z t hj 5 20 Tt The temperature at each lattice point relative to a mean value can be determined as the sum of the distribution functions T hi which can be used to determine the equilibrium distribution function e s u again using the bulk fluid velocity The thermal relaxation time is related to the thermal diffusivity vl 1 Az e T M with the Prandtl number for the system determined by a ratio of relaxation times i e P ve 2rf l1 CT 27 1 It should be noted that if multiple relaxation time MRT schemes are to be used these only apply to fluids all diffusion and heat transfer processes are calculated using the single relaxation
235. ter of this Licence Agreement or the use of the DL_MESO Software will not exceed in total 5000 Nothing in this Licence Agreement limits or excludes STFC s liability for death or personal injury caused by its negligence or for any fraud or for any sort of liability that by law cannot be limited or excluded The express undertakings and given by STFC in this Licence Agreement and the terms of this Licence Agreement are in lieu of all warranties conditions terms undertakings and obligations on the part of STFC whether express or implied by statute 197 common law custom trade usage course of dealing or in any other way All of these are excluded to the fullest extent permitted by law 4 INTELLECTUAL PROPERTY RIGHTS AND ACKNOWLEDGEMENTS 4 1 4 2 5 2 5 3 5 4 Nothing in this Licence Agreement assigns or transfers any Intellectual Property Rights in any of the DL_MESO Software Those rights are reserved to STFC The Licensee will ensure that if any of its employees or students publishes any article or other material resulting from or relating to a project or work undertaken with the assistance of any part of the DL_MESO Software that publication will contain the following acknowledgement DL_MESO is a mesoscale simulation package written by R Qin W Smith and M A Seaton and has been obtained from STFC s Daresbury Laboratory via the website http www ccp5 ac uk DL_MESO TERMINATION This
236. this frequently produces a system temperature higher than that specified by the user and requires a small time step At to reduce the offset to tolerable levels 10 2 2 DPD thermostat with DPD Velocity Verlet integration DPD VV itype 1 As with the MD VV scheme this algorithm uses the drag and random forces as the system thermostat which are combined with all other forces before being integrated using the Velocity Verlet scheme The drag force is subsequently recalculated after the second stage using the velocities at the end of the time step 10 The recalculation of drag forces after force integration helps to alleviate the temperature offset produced by the MD VV and hence larger time steps may be used for reasonable temperature control It does require the re use of the linked list cells and inter processor communications to recalculate the drag forces which can significantly increase the time required per time step compared to the MD VV scheme 10 2 3 Lowe Andersen thermostat itype 2 The Lowe Andersen thermostat 27 is an alternative to the use of drag and random forces in the DPD thermostat which uses a variant of the Andersen thermostat 1 After all other forces conservative bonding etc are integrated using the Velocity Verlet scheme a random sample of particle pairs have their relative velocity replaced by a value from a Maxwellian distribution i e kgT Us Cig 10 4 j ij where Hij eee is the reduced mas
237. ticle numbers global1 and global2 that produce the shortest distance between the particles e Dependencies None e Arguments globali input integer global2 input integer ind1 input output integer ind2 input output integer contract _bndtbl e Header records SUBROUTINE contract_bndtbl e Function Strips out all bond pairs from bond table that have been reassigned to neighbouring processors e Dependencies None e Comments Only called for parallel version of DL_MESO_DPD when bond tables include only local bonds in each processor contract _angtbl e Header records SUBROUTINE contract_angtbl 154 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE e Function Strips out all bond angle triples from angle table that have been reassigned to neighbouring processors e Dependencies None e Comments Only called for parallel version of DL_MESO_DPD when angle tables include only local angles in each processor contract_dhdtbl e Header records SUBROUTINE contract_dhdtbl e Function Strips out all bond dihedral quadruples from dihedral table that have been reassigned to neighbouring processors e Dependencies None e Comments Only called for parallel version of DL_ MESO_DPD when dihedral tables include only local dihedrals in each processor bond_force e Header records SUBROUTINE bond_force bondtype r a b c d force potential e Function Determines the stretching force and potential energy be
238. time schemes described in this section 48 CHAPTER 5 DL_MESO_LBE FEATURES 5 5 Compressible and incompressible fluids The standard Lattice Boltzmann Equation scheme is capable of modelling compressible fluids Incompressible fluids can be modelled by making a simple modification to the local equilibrium distribution function 60 B e ir 2 0 2 fi wi f po 2 2 ds e 5 22 2c4 2c where po is the fixed density of the incompressible fluid and the density p becomes an analogue to pressure c2p While Equation 3 1 is still applicable to calculate p the fluid velocity is now calculated by q Poua ys Filia 5 23 i 0 The initialization collision and file output routines therefore need to be modified for incompressible fluid systems DL MESO LBE supplies additional routines ending in Incom to allow the user to model incompress ible fluids e g the subroutine call for the default BGK collision scheme fCollisionBGK can be changed to fCollisionBGKIncom All of these routines use a specified constant density for each fluid in the system 1bincp as the value of p Chapter 6 DL MESO LBE Package Reference 6 1 Overview DL_MESO LBE consists of seven packages e 1bpMODEL Contains subroutines to assign lbw 1bv and 1bopv for D2Q9 D3Q15 D3Q19 and D3Q27 lattice models e 1bpBASIC Contains general purpose subroutines which can be used elsewhere e g a random number generator producing values between 1 and 1
239. total number of unbonded particles in system nusyst integer total number of molecules in system nummol integer total number of bonds in system numbond integer total number of bond angles in system numang integer total number of bond dihedrals in system numdhd integer number of defined bond types nbonddef integer lt mxbonddef number of defined bond angle types nangdef integer lt mxbonddef number of defined bond dihedral types ndhddef integer lt mxbonddef time step number nstep integer force calculation time accumulator timfrc real KIND dp step time accumulator timstp real KIND dp specified system temperature kgT temp real KIND dp size of time step At tstep real KIND dp halo boundary size rhalo real KIND dp interaction cutoff radius re rcut real KIND dp gt 0 square of interaction cutoff radius rct2 real KIND dp gt 0 many body DPD cutoff radius rg rmbcut real KIND dp square of many body DPD cutoff radius rmbct2 real KIND dp short range electrostatic cutoff re relec real KIND dp square of electrostatic cutoff rel2 real KIND dp maximum calculation time timjob real KIND dp calculation close time tclose real KIND dp system volume volm real KIND dp size of system in x dimension dimx real KIND dp size of system in y dimension dimy real KIND dp size of system in z dimension dimz real KIND dp number of domain cells in x npx integer gt 1 number of domain cells in y npy integer AT number of domain cells i
240. trix of m collision parameters represented by 5 A diag s 5 6 Some of the collision parameters can be specified by the user to set both kinematic and bulk viscosities a few others can be tuned to improve simulation stability and the remainder i e those for density and momentum are fixed to conserve macroscopic hydrodynamics If all values of s are set to a the scheme reduces to BGK single relaxation time collisions Since the collisional matrix is diagonal equation 5 5 can be rewritten in terms of each moment i e Mi amp t M E t s M E t ME E 0 5 7 2 Multiplying M Z t by the inverse of the transformation matrix T gives the post collisional distribution functions Interfacial and external forces can be applied either by the addition of rF to the fluid momentum 32 or by the use of a moment transformed source term S whose terms are defined by 39 s u EAF 5 8 and applied by the following to correct the post collisional moments AM I A SAt 5 9 where I as an identity matrix The above equation can be re expressed as AM 1 tsi S At 5 10 The above equations after inverse transformation reduce to equation 5 2 when the collision parameters are set to a All collision parameters for conserved moments should be set to unity when applying external forces 5 1 COLLISION AND PROPAGATION ALGORITHMS 43 An example can be given for the D2Q9 lattice system 26 the momen
241. ture to calculate local densities for each component pi gt we r A omitting self contributions for each particle 53 weight_rho e Header records REAL KIND dp FUNCTION weight_rho rrr e Function Calculates normalized weight function for local densities e Dependencies None e Comments The default weight function 55 is which may be changed by the user manybody_force e Header records SUBROUTINE manybody_force i j k rrr mbforce e Function Calculates many body DPD interaction force and non density dependent potential energies between par ticles i and j using parameter set k and inter particle distance rrr e Dependencies None 152 CHAPTER 11 DL MESO DPD PACKAGE REFERENCE e Arguments i input integer j input integer k input integer rrr input real KIND dp mbforce output real KIND dp e Comments The default form for the many body DPD interaction force is a two term style suitable for modelling vapour liquid mixtures 55 mC _ 4 4_ TH peep raj Ps FE jAy 1 2 Biz pi ps 1 3 ne This subroutine may be changed by users who wish to use different many body interaction functional forms manybody_potential e Header records SUBROUTINE manybody_potential e Function Calculates the self energy for every particle resulting from density dependent many body interactions Dependencies None Comments The default form is based on the two term vapour li
242. tween a pair of bonded particles e Dependencies None e Arguments bondtype input integer r input real KIND dp a input real KIND dp b input real KIND dp c input real KIND dp d input real KIND dp force output real KIND dp potential output real KIND dp e Comments If required the user can add extra stretching bond interaction types in this subroutine as additional cases angle_force e Header records SUBROUTINE angle_force angtype theta rab rcb a b c d force potential virial dfab df cb 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 155 e Function Determines the bond angle force potential energy and virial across a triple of bonded particles e Dependencies None e Arguments angtype input integer theta input real KIND dp rab input real KIND dp rcb input real KIND dp a input real KIND dp b input real KIND dp c input real KIND dp d input real KIND dp force output real KIND dp potential output real KIND dp virial output real KIND dp dfab output real KIND dp dfcb output real KIND dp e Comments If required the user can add extra bond angle interaction types in this subroutine as additional cases dihedral_force e Header records SUBROUTINE dihedral_force dhdtype phi pb pc a b c d force potential e Function Determines the bond dihedral force and potential energy across a quadruple of bonded particles e Dependencies None e Arguments d
243. uble precision fGetOneSpeedSite e Header records three cases int fGetOneSpeedSite double speed double startpos int fGetOneSpeedSite double speed int fpos long tpos int fGetOneSpeedSite double speed int fpos int xpos int ypos int zpos e Function Calculates the macroscopic speed of compressible fluid fpos at the grid point e Dependencies None e Arguments startpos input double pointer fpos input integer tpos input long integer xpos input integer ypos input integer Zpos input integer speed output double precision array e Comments The second and third cases are slightly slower than the first The calculation is based on vg auis but more readable fGetOneSpeedIncomSite e Header records three cases int fGetOneSpeedIncomSite double speed double startpos double rho0 int fGetOneSpeedIncomSite double speed int fpos long tpos int fGetOneSpeedIncomSite double speed int fpos int xpos int ypos int zpos e Function Calculates the macroscopic speed of incompressible fluid fpos at the grid point e Dependencies None 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 61 e Arguments startpos input double pointer rho0 input double precision fpos input integer tpos input long integer xpos input integer ypos input integer Zpos input integer speed output double precision array e Comments The calculation is based on va Pa fits The second and third cases are slightly slowe
244. uid It is particu larly suited to simulating systems with complex boundary conditions and is also suitable for systems with phase transitions since mesoscale interactions can be merged into the method easily Chapter 4 DL MESO LBE Basic Definition 4 1 Lattice models DL_MESO_LBE utilizes a right handed Cartesian coordinate system with the x axis from left to right in the horizontal direction the y axis from low to high in the vertical direction and z axis from back to front D2Q9 D3Q15 D3Q19 and D3Q27 lattice models have been included The speed vectors and weight factors are arranged to allow the use of swap algorithms for propagation 34 The transformation matrices T for Multiple Relaxation Time MRT schemes are included for each lattice except D3Q27 along with the equilibrium moments M expressed for the incompressible case for the compressible case pg is substituted with p gt collision operators 5 definition of bulk viscosity 1 and forcing source terms 5 D2Q9 Weight factor 2 Wi 0 4 2 4 6 8 i 153 957 z Speed vector 1 Cia Ciy 0 0 0 1 1 1 2 1 0 3 1 1 4 0 1 5 1 1 6 1 0 7 1 1 8 0 1 25 26 CHAPTER 4 DL_MESO_LBE BASIC DEFINITION 1 1 1 1 1 1 i ame A l 4 2 1 2 1 2 1 2 I 4 1 2 12 1 2 1 2 0 1 1 1 0 1 1 1 0 T 0 1 2 1 0 1 2 1 0 0 1 0 1 1 1 0 1 1 0 1 0 1 2 1 0 1 2 0 0 1 0 1 0 1 0
245. ulation The utility should be compiled beforehand to give the executable molecule exe refer to Appendix B for more details Each input file for molecules FIELD and MOLECULE can be examined and or edited using DL_MESO s built in editor which is launched by the button edit file In addition to identifying the particles with a molecule that are bonded together and where necessary make up any angles or dihedrals the FIELD file also requires the parameters for each bond interaction The pull down list allows the user to select bond angle or dihedral interactions The interaction type is selected using this pull down list whose contents are dependent on whether a bond angle or dihedral is required 2 4 COMPILING AND RUNNING DL_MESO 17 The parameters are entered in these boxes corresponding to Tables 10 4 10 5 and 10 6 for bonds angles and dihedrals respectively The button add to FIELD will add the interaction and its parameters to the existing FIELD file The molecule name of up to 8 characters in length and molecule population for each molecule described in the input files is required in the CONTROL file clicking add molecule stores the entered information in memory in the order it will appear in the CONTROL file while reset molecules removes them all from memory By default a maximum of 8 molecule types can be included in DL_MESO_DPD The types of bonded interactions can be selected and specified in t
246. unction Reads calculation parameters lattice scheme numbers of fluids solutes temperature scalars and phase field order parameters grid size from input system file 1bin sys e Dependencies lbin sys data file e Arguments filename input array of characters e Comments The default file name is 1bin sys to use different file names the argument for this routine can be changed by the user fPrintDomainMass e Header records int fPrintDomainMass e Function Calculates and prints total and individual fluid masses in domain e Dependencies double fGetTotMassDomain double fGetOneMassDomain e Comments This routine only produces the masses for the entire system if running in serial to obtain system masses in parallel running fPrintSystemMass would be needed fPrintDomainMomentum e Header records int PrintDomainMomentum 6 2 DLMESO_LBE SUBROUTINES AND FUNCTIONS 65 e Function Calculates and prints the total fluid momentum in domain e Dependencies int fGetTotMomentDomain double momentum e Comments This routine only produces the entire system momentum if running in serial to obtain system momentum in parallel running fPrintSystemMomentum would be needed fInput Parameters e Header records int fInputParameters const char filename lbin sys e Function Reads system parameters in from 1bin sys data file e Dependencies lbin sys data file e Arguments filename i
247. ure and phase field the distribution functions are in the order of ida tattoo tatoo Les Ji f To Ti To To Ta To To Tr Ts pf Therefore the number of particle distribution functions at each grid point is lbsitelength lbsy nf lbsy nc lbsy nt x 1bsy nq 1bsy ph where 1bsy nf lbsy nc 1bsy nt 1bsy nq and 1bsy ph are respectively number of fluids number of solutes number of temperature scalars number of discrete speeds and number of the phase field order parameters lbsy nt and lbsy ph can only take the values of 1 or 0 representing systems with or without temperature scalar and phase field Also if 1bsy nc 4 0 lbsy nf cannot be set larger than 1 4 2 2 Storage of space property The space property is represented by an integer value in DL MESO_LBE For example 1bphi 100 0 represents the 100th grid point as a liquid site and 1bphi 101 12 shows the 101st grid point as an on grid bounce back boundary Table 4 1 lists the categories of space properties The orientation of a solid liquid boundary is also represented by the value of an integer For example a planar surface with normal vector along the y axis is denoted by 21 while a concave corner face at the top right front corner is denoted by 31 It must be pointed out that only those space positions located in the surface of a face centered cube have been included and translated in DL_MESO_LBE Points with random orientations e g 47 plane have not been included
248. utput real KIND dp vdd input output real KIND dp sided input real KIND dp e Comments Momentum tangential to the boundary is preserved i e this applies a free slip boundary condition 11 2 7 ewald module This module requires the constant variables numeric_container and comms_module modules to be loaded beforehand ewald_real e Header records SUBROUTINE ewald_real nlimit e Function Calculates real space terms for Ewald summation with Slater type exponential decay charge distribu tions e Dependencies erfcdp e Arguments nlimit input integer e Comments Calculates short range Coulombic forces and potential energies for Ewald summation using the following smeared non point charge distribution 11 to amp exp 2 Te ewald_reciprocal e Header records SUBROUTINE ewald_reciprocal e Function Calculates reciprocal space terms for Ewald summation e Dependencies None 11 2 DL_MESO_DPD SUBROUTINES AND FUNCTIONS 151 e Comments Calculates long range Coulombic forces and potential energies using standard Ewald summation including self energy corrections for charged systems 11 2 8 manybody_module This module requires pre loading of the constants and variables modules local_density e Header records SUBROUTINE local_density e Function Calculates local densities for many body DPD interactions e Dependencies weight_rho e Comments Uses the parallel link cell struc
249. y also be argued that the self diffusion of DPD particles might not correspond to that of individual molecules and thus a Schmidt number of the order 10 is unnecessary for modelling liquids 38 Alternative thermostats are available which can model systems with higher Schmidt numbers 27 49 8 4 Derivation of Equilibrium The derivation of the DPD algorithm is based on the Fokker Planck equation Op C 8 17 a LP 8 17 where p is the equilibrium distribution function and is the evolution operator which may be split into conservative and dissipative parts Estel 8 18 114 CHAPTER 8 DISSIPATIVE PARTICLE DYNAMICS BASIC THEORY with INS 3 N AS eee Lo pee pas Oe 8 19 i 1 tj A a af p 7 2 o 0 Lo Yeu ae om i ig gt wh ras isle ae 8 20 e Tis where jj 2 ij When a y 0 then Equation 8 17 becomes ane Lop 8 21 for which the equilibrium solution is evidently N je ix pt E exp EBT 2 om 7 26 rs 8 22 which is of course the Boltzmann distribution function for an equilibrium system Thus it is apparent that for the simulation based on Equation 8 17 to maintain the same distribution function the terms in the operator LP of Equation 8 20 must sum to zero It follows that the conditions given in Equations 8 8 and 8 9 must apply 8 5 Summary of Dissipative Particle Dynamics DPD is a simple method All that is required is a system of spherical particles enclosed in a peri
250. y but more easily handle the simulation of three dimensional systems with complex boundary conditions The equilibrium for D2Q9 and D3Q27 lattice models can be derived a priori from the Maxwell equilibrium distribution D3Q15 and D3Q19 models appear to be more popular than D3Q27 because the latter is much more expensive in terms of computing cost It is required that the equilibrium state should be able to reproduce elementary macroscopic fluid variables p Yi 3 1 22 CHAPTER 3 THE LATTICE BOLTZMANN EQUATION BASIC THEORY q Pua 5 filia 3 2 i 0 where p is the density f the ith particle distribution function the ith speed vector and ua the macroscopic velocity along the a axis In the Lattice Boltzmann method is not the thermal velocity of a particle and therefore EX Y hea 33 Equation 3 3 implies that the temperature cannot be derived from the lattice particle distribution function However temperature can be defined at each grid point 3 3 Derivation of Equilibrium There are two methods by which the Lattice Boltzmann equilibria can be constructed The bottom up method obtains the equilibrium from the Maxwell Boltzmann equilibrium distribution The top down method constructs the equilibrium so that the required macroscopic properties can be reproduced Only the bottom up method is shown here the top down method can be found in 5 The Maxwell Boltzmann single particle equilibrium distribution function is

DL MESO USER MANUAL

Contents

Download Pdf Manuals

Related Search

Related Contents