D2.6 - OBMD library for LAMMPS

OBMD library for LAMMPS MultiXscale Deliverable 2.6 Deliverable Type: Report Delivered in June, 2025 MultiXscale EuroHPC Centre of Excellence for Multiscale Modelling Acknowledgement Funded by the European Union. This work has received funding from the European High Performance Computing Joint Undertaking (JU) under grant agreement No 101093169. Disclaimer Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European High Performance Computing Joint Undertaking (JU). Neither the European Union nor the granting authority can be held responsible for them. MultiXscale Deliverable 2.6 Page ii Project and Deliverable Information Project Title MultiXscale: EuroHPC Centre of Excellence for Multiscale Modelling Project Ref. Grant Agreement 101093169 Project Website https://www.multixscale.eu EuroHPC Project Officer Matteo Mascagni Deliverable ID D2.6 Deliverable Nature Report Dissemination Level Public Contractual Date of Delivery Project Month 30 (30th June, 2025) Actual Date of Delivery 27th June, 2025 Description of Deliverable OBMD library with plugin for LAMMPS (including user guide, programmer’s manual, architectural document) Document Control Information Document Title: OBMD library for LAMMPS ID: D2.6 Version: As of June, 2025 Status: Accepted by Steering Committee Available at: https://www.multixscale.eu/deliverables Document history: Internal Project Management Link Review Review Status: Reviewed Authorship Written by: Petra Papež (NIC) Contributors: Tilen Potisk (NIC), Nikolaos Ntarakas (NIC), Matej Praprotnik (NIC) Reviewed by: Tilen Potisk (NIC), Matej Praprotnik (NIC) Approved by: Matej Praprotnik (NIC) Document Keywords Keywords: MultiXscale, OBMD , equilibrium , shear flow , acoustic wave 27th June, 2025 Disclaimer: This deliverable has been prepared by the responsible Work Package of the Project in accordance with the Consortium Agreement and the Grant Agreement. It solely reflects the opinion of the parties to such agreements on a collective basis in the context of the Project and to the extent foreseen in such agreements. Copyright notices: This deliverable was co-ordinated by Petra Papež1(NIC) on behalf of the MultiXscale consortium with contributions from Tilen Potisk (NIC), Nikolaos Ntarakas (NIC), Matej Praprotnik (NIC). This work is licensed under the Creative Commons Attribution 4.0 International License. To view a copy of this license, visit: http://creativecommons.org/licenses/by/4.0 cb 1[email protected] MultiXscale Deliverable 2.6 Page iii Contents Executive Summary 1 1 Introduction 2 1.1 Scope of the deliverable ................................................ 2 1.2 Target audience . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Report outline ...................................................... 2 2 Overview of the OBMD method 3 3 Implementation of OBMD method in LAMMPS 5 4OBMD simulations of liquid water employing DPD and SPC water models 8 4.1 DPD water model .................................................... 8 4.2 SPC water model ..................................................... 10 5 Scaling study 12 6 Conclusion and perspectives 13 A Appendices 14 A.1 OBMD User Guide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 A.2 OBMD Programmers Manual . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 A.3 OBMD Architectural Document . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 Acknowledgements 53 References 53 List of Figures 1 Schematic representation of the OBMD simulation setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Schematic representation of the weighting functions g∥and g⊥that are used to distribute normal and tangential forces, respectively. ............................................ 4 3 Directory tree showing the contents of the OBMD package and the newly added pair style of the KOKKOS package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 4 Schematic illustration of the timestep and functions called in the pre_exchange() method. Blue colored rectangle contains the newly implemented function (blue text) and variables defined (black text) in the obmd extension. To increase readability, comments are also provided (gray text). . . . . . . . . . . . . 5 5 Schematic illustration of the timestep and functions called in the post_force() method. Same color coding is used as in Fig. 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 6 Number density profile along the open direction of the simulation box. The desired density of the DPD water is depicted using the black dotted line and is equal to ρ=3.0 U−3 L. The two orange dotted lines indicate boundaries buffer-ROI.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 7RDF for the DPD water. Black line and dots correspond to the RDF computed from NVT simulation, while blue line and dots depict the RDF computed from the OPEN-BOUNDARY MOLECULAR DYNAMICS (OBMD) simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 8 Computed density variation through the ROI for a simulated acoustic wave with frequency of 1.84 U−1 T and ∆p=0.5pxx at times t=0 (left panel), t≈t0/2 (middle panel), and t=t0(right panel). The analytical solution is shown using the blue solid line, while black dots represent simulation data. . . . . . . . . . . . 9 9 Calculated velocity profile vyalong the open direction of the simulation box at pxy =0.8 UEU−3 L. The ROI is sandwiched between the two black dotted lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 10 Mass density profile along the open direction of the simulation box. The desired density of the SPC water is depicted using the black dotted line and is equal to ρ=1.2 UMU−3 L. The two orange dotted lines indicate the boundaries buffer-ROI.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 11 RDFs for the SPC water. Black line and dots correspond to RDF calculated from the NVT simulation, while blue line and dots depict RDF computed from OBMD simulations. In both cases, the rf method was used to describe electrostatic interactions beyond the cutoff. Black dotted lines denote distances and corresponding values of g(r) reported in the literature for the SPC water. . . . . . . . . . . . . . . . . . . . 11 12 Computed density variation through the ROI for simulated acoustic wave with frequency of 0.21 U−1 T and ∆p=2.0pxx at times t=1/6t0(left panel), t≈1/3t0(middle panel), and t=t0(right panel). The analytical solution is shown using the blue solid line, while black dots represent simulation data. . . . . 11 MultiXscale Deliverable 2.6 Page iv 13 Calculated velocity profile vyalong the open direction of the simulation box at pxy =2.0 UEU−3 L. The ROI is located between the two black dotted lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 14 Strong scalability study performed using the DPD fluid. The dashed lines indicate results obtained when the PBCs are applied, while solid lines show data collected using the obmd extension. . . . . . . . . . . . . 12 MultiXscale Deliverable 2.6 Page 1 Executive Summary This deliverable provides a comprehensive overview of the open-boundary molecular dynamics (OBMD) method and its implementation within the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) simulation package. Typically, molecular dynamics (MD) simulations are conducted in the canonical ensemble, which maintains a constant number of particles, volume, and temperature. However, many real systems are open and can exchange mass, momentum, and energy with their surroundings. This necessitates the use of the grand-canonical ensemble to perform simulations. The OBMD method makes this possible by allowing particle exchange through the opening of the system’s boundaries. To this end, the OBMD method enables imposition of the external boundary conditions defined by the normal component of the momentum flux, where the governing equations of motion in the bulk remain intact. The implementation of the OBMD method introduces a new fix style called obmd, which manages the deletion and insertion of particles, measures the outflow of linear momentum of deleted particles, and applies external forces to particles in the buffer regions, thereby enabling the imposition of the user-defined boundary conditions. The execution of the newly implemented obmd extension is invoked by specific methods within the LAMMPS code. To demonstrate the applicability of the newly implemented obmd extension, we performed both equilibrium and nonequilibrium simulations of liquid water using the mesoscopic dissipative particle dynamics (DPD) and simple point charge (SPC) water models. Our findings clearly show that open boundaries do not affect the structural properties of the liquid. Furthermore, we showcased the application of the obmd extension as a virtual ultrasound machine and virtual rheometer, emphasizing the significant potential of the OBMD method for biomedical applications. In addition, we investigated the strong scalability of the DPD fluid in equilibrium using the LAMMPS code with the obmd extension enabled. We compared these results to simulations where periodic boundary conditions (PBCs) were utilized. The results demonstrated that the introduction of the obmd extension does not significantly affect the strong scaling behavior of the LAMMPS code. We observed excellent scalability up to at least 80 A100 GPUs, which represents approximately one-third of the EuroHPC Vega GPU partition. The observed behavior is crucial for large-scale simulations, ensuring that our methods remain efficient when applied on advanced supercomputing infrastructures. The appendices to this deliverable include essential resources such as the user guide, programmer’s manual, and architectural document, all designed to further support users and developers of the obmd extension. The user guide offers step-by-step instructions for compiling LAMMPS with the obmd extension and conducting OBMD simulations. The programmer’s manual provides insight into the structure and functions of the new header and implementation files, while the architectural document provides details about the communication patterns associated with the newly added obmd fix style. MultiXscale Deliverable 2.6 Page 2 1 Introduction 1.1 Scope of the deliverable This deliverable outlines the OBMD method, which opens the system boundaries and allows imposition of the external boundary conditions, while leaving the governing equations of motion in the bulk intact. It includes a detailed description of the implementation of the OBMD method as a newly integrated obmd extension within the LAMMPS code (version: 7Feb2024). To demonstrate the versatility of applications enabled by the OBMD method, we present and discuss the results obtained by performing equilibrium and nonequilibrium OBMDsimulations employing DPD and SPC water model. Additionally, the deliverable features a user guide and a programmer’s manual. The user guide provides comprehensive instructions on compiling LAMMPS with the obmd extension and conducting OBMD simulations with the provided Python script that generates an input script for LAMMPS. The programmer’s manual provides information on the structure and functions of the added header and implementation files, which are utilized to execute OBMD simulations. The architectural document addresses the communication patterns and key components of the new obmd fix style. 1.2 Target audience This deliverable is aimed (but not limited) to: • people employing MD simulations, • people using LAMMPS MD simulation code, • people who want to exploit OBMD method to perform simulations, • people who want to further extend the OBMD method, • people who want to test scaling of the system. 1.3 Report outline This deliverable is structured as follows. Sec. 2provides a detailed description of the OBMD method and gives an insight into its implementation within the LAMMPS code. The primary focus of Sec. 3is to illustrate the execution of a simulation timestep and to identify the stages at which the functions of the obmd extension are invoked. In Sec. 4, we define the external boundary conditions required for the imposition of a constant normal load, inducing shear flow, or exciting acoustic waves. In the same section, we present and discuss the results obtained from the OBMD simulations utilizing the specified boundary conditions. Finally, Sec. 5examines the impact of the obmd extension on the scaling behavior of the LAMMPS code. Additionally, the user guide (Appendix A.1), programmer’s manual (Appendix A.2), and architectural document (Appendix A.3) are included at the end of this deliverable, providing further resources for users and developers. MultiXscale Deliverable 2.6 Page 3 2 Overview of the OBMD method The MD simulations are usually performed in the canonical (NVT) ensemble, mimicking the experimental conditions of constant temperature (T), volume (V), and number of particles (N). However, many real systems are open; therefore, they exchange the mass, momentum, and energy with the external world. To characterize such systems, the grand-canonical (µVT) ensemble is used. The latter is specified by the chemical potential (µ), volume (V), and temperature (T). To this end, grand-canonical MD simulations with a fluctuating number of particles are enabled by the OBMD method [1–3], which by opening of the system’s boundaries allows for particle exchange. In OBMD, the simulation box is opened in one direction, while PBCs are applied in the remaining ones. The box is divided into three regions, where the central region of interest (ROI) is enclosed by two buffer regions (as depicted in Fig. 1). Buffers Figure 1: Schematic representation of the OBMD simulation setup. serve as particle reservoirs, where particles are removed from and inserted into the system. The number of particles in the buffers is maintained by the following feedback algorithm: ∆Nb=δt τb (〈Nb〉− Nb), (1) where 〈Nb〉and Nbare the desired and current numbers of particles inside the buffer, respectively, while τbrepresents the relaxation time of the buffer. When there are more particles in the buffer than 〈Nb〉, that is, when ∆Nb<0, the particles need to be deleted from the system. In practice, particles are not deleted immediately, but rather left to diffuse across the outer boundary of the buffer and then deleted. In contrast, new particles are inserted when ∆Nb>0, implying that there are fewer particles inside the buffer than desired. The insertion of new particles is carried out using the iterative USHER algorithm [4]. The algorithm begins by selecting a random position within the buffer for the particle to be inserted. Then, the potential energy of the particle is calculated relative to all other particles in the buffer. If the calculated potential energy is below a predefined threshold (U0), the position is accepted and the particle is inserted. On the other hand, if the energy is higher, the position is iteratively corrected according to the following update rule: rn+1=rn+fn |fn|δsn. (2) Here, rnis the position of the particle at n-th iteration and fnis the total force on a given particle in the same iteration step. The displacement δsndepends on the local potential energy Unand force magnitude |fn|. If the computed energy Unis larger than Uovlp , which is chosen to be a very large energy, representing an overlap position, δsnis computed using δsn=∆sovlp =rσ−(4ε/Un)1/2. Here, the position of the particle will be translated for a distance rσ from the center of mass of the overlapped particle. When the calculated energy is lower than Uovl p ,δsis determined using δsn=min(︁∆s,(Un−U0)/|fn|)︁. The iterative procedure is performed until the computed energy is below the predefined threshold. If this does not succeed within the specified number of steps, the position is rejected, and the insertion is repeated by selecting another random position. The buffers are also utilized to introduce external boundary conditions, while keeping the equations of motion unchanged. To impose boundary conditions, OBMD employs external forces fext iacting only on particles within buffers, i.e., fext i=0 everywhere else [3]. In addition, in OBMD the total linear momentum is conserved, which directly stems from the Navier-Stokes equation, coming on top of the linear momentum conservation law ∂(ρv) ∂t= −∇·Jp. (3) Here, ρ,v, and Jpare the fluid density, fluid velocity, and momentum flux tensor, respectively. The latter is expressed as Jp=ρv⊗v+Π, (4) MultiXscale Deliverable 2.6 Page 4 where Πis the stress tensor, given by Π=(P+π)I+Π. (5) The first term on the right-hand side of Eq. 5is diagonal, where Pis the thermodynamic pressure, and π= −ζ∇·v, i.e., the isotropic stress. The second term is a traceless symmetric tensor, whose components are Παβ = −η(∂αvβ+∂βvα−2∂γvγδαβD−1). (6) Here, Dstands for spatial dimension, while symbols ηand ζdenote the shear and bulk viscosity coefficient, respectively. In OBMD, boundary conditions are defined by the normal component of the momentum flux, i.e., Jp·n, where Jp stands for the already defined momentum flux tensor and nis the unit vector normal to the interface buffer-ROI (pointing towards the center of the ROI). To determine external forces used to impose the external boundary conditions, the amount of momentum created by such forces over one timestep (δt) should be evaluated, and the result equated to the desired amount of momentum that needs to be added to or extracted from the system. By reformulating the momentum flux balance equation for the OBMD boundary of the surface A, Jp·nAδt=∑︂ i∈B fext iδt+∑︂ i′ ∆(mi′vi′), (7) we arrive at the expression for the total external force Fext : Fext =∑︂ i∈B fext i=A⎛ ⎜ ⎝Jp·n− ∑︂ i′∈B ∆(mi′vi′) Aδt ⎞ ⎟ ⎠. (8) Indices iand i′run over all the particles in the buffers and over all particles that have been inserted into or deleted from the system in the last timestep, respectively. Accordingly, the momentum change is given by ∆(mi′vi′)= ±mi′vi′ if the particle i′is inserted (+) or deleted (-), where mi′and vi′represent the mass and velocity of the particle i′, respectively. Ultimately, the total external force Fext is distributed among the particles in the buffer using fext i=G(ri)Fext , (9) where G(ri) is a weighting function, generally in tensorial form, allowing distribution of normal and tangential forces (see Fig. 2) [1,3]. Figure 2: Schematic representation of the weighting functions g∥and g⊥that are used to distribute normal and tangential forces, respectively. MultiXscale Deliverable 2.6 Page 5 3 Implementation of OBMD method in LAMMPS In this section, we describe how the OBMD method (presented in Sec. 2) is implemented in the LAMMPS code by closely inspecting: • the newly added OBMD package (see Fig. 3), which includes the new obmd extension for LAMMPS simulation package, • deletion and insertion of particles, • measurement of the outflow of the linear momentum of particles that have diffused across the outer boundary of the buffers, • imposition of the external boundary conditions using external force. Figure 3: Directory tree showing the contents of the OBMD package and the newly added pair style of the KOKKOS package. The execution of the obmd extension within a timestep is triggered by the pre_exchange() and post_force() methods (see Fig. 4and Fig. 5). In the pre_exchange() method, particles are inserted or deleted before rebuilding the neighbor list. In the pre_exchanged() method, the try_deleting() member function is called first. Here, particles that have left the outer boundary of the buffers are marked for deletion. In the process of deletion, a loop is made over Figure 4: Schematic illustration of the timestep and functions called in the pre_exchange() method. Blue colored rectangle contains the newly implemented function (blue text) and variables defined (black text) in the obmd extension. To increase readability, comments are also provided (gray text). MultiXscale Deliverable 2.6 Page 12 5 Scaling study In this section, we investigate whether the obmd extension affects the scaling behavior of the LAMMPS code. 48 16 32 64 80 No. of GPUs 10−3 10−2 10−1 100 101 Seconds per timestep 255k 525k 1M 2M 4M 8M 33M 1 Figure 14: Strong scalability study performed using the DPD fluid. The dashed lines indicate results obtained when the PBCs are applied, while solid lines show data collected using the obmd extension. To this end, we have conducted a strong scalability study of the DPD fluid in equilibrium using the obmd extension and compared it with the NVT simulations, where PBCs were applied. The results obtained demonstrate that the strong scalability behavior of the LAMMPS code is not greatly affected by the introduction of the obmd extension. In addition, we observed an excellent strong scalability behavior up to at least 80 A100 GPUs, which is about one-third of the EuroHPC Vega GPU partition. To put this into perspective, each A100 GPU delivers approximately 10 TFlops in double precision and 20 TFlops in single precision. This means that the 80 GPU workload we tested was running on hardware with raw compute capabilities between 0.75 and 1.5 Petaflops, highlighting the immense computational power used. Such scalability is crucial for large-scale simulations, as it ensures that our methods remain efficient even when applied on advanced supercomputing infrastructures. MultiXscale Deliverable 2.6 Page 13 6 Conclusion and perspectives In conclusion, this deliverable presented a comprehensive overview of the OBMD method and its integration into the LAMMPS code (version 7 Feb 2024). The newly added obmd extension will allow performing OBMD simulations in and out-of-equilibrium, as demonstrated in this deliverable. At the same time, we showed that scalability behavior of the LAMMPS code is not greatly affected by the introduction of the presented obmd extension. We aim to integrate the presented obmd extension into the official LAMMPS release, upon successful completion of Tasks 3.2 and 4.4. The latter are specifically designed to enhance the functionality of the OBMD method. Prior to this integration, we will deploy our builds and example scripts to dev.eessi.io, i.e., the development repository of EESSI. This will make it easier to test development versions of our LAMMPS code extensions on different system architectures, where dev.eessi.io CernVM-FS is available. This document serves as a guide for users and developers working with the o b md extension in LAMMPS. Therefore, general knowledge about LAMMPS and molecular dynamics simulations is assumed. Additionally, basic familiarity with compiling and running C++ code is expected. To perform open-boundary molecular dynamics (OBMD) simulations using the o b md extension, we provide a Python script named input . py , which generates in . simul a tion input script for LAMMPS. LAMMPS is then started with this input as a parameter, for example, The script input . py is located here. INITIALIZATION SETTINGS A code snippet of the function in the input . py Python script, which writes the input script for the LAMMPS OBMD simulation of a DPD fluid, is given and explained below. First, the generated output file name, which serves as the input for LAMMPS ( in _ file ), is set. The default value is in . simul a tion . Next, the write _ in ( in _ file ) function is defined, taking the name of the input file to be created as an argument. Using an f-string, the block of code that forms the input script is stored in c ontent _ 4 in . Finally, the content is written to a file using file _ to _ write . write ( c ontent _ 4 in ) . OBMD simulations should be performed using reduced units ( units lj ) to avoid errors in the calculation of the external force and the style of boundaries ( b ound a ry ) for the global simulation box in each dimension should be set to f p p . The latter applies the same style to the lower and upper sides of the simulation box, i.e., in the case of OBMD simulations, the non-periodic style f is applied in the x -direction, while the box is periodic in y - and z - directions by choosing the periodic style p . Following is the explanation of the inputs. lmp _ mpi - in in . simul a tion in _ file = " in . simul a tion " # -------------------------------------------------------- # def write _ in ( in _ file ): file _ to _ write = open ( in _ file , ' w ' ) c ontent _ 4 in = f "" " # ----------------- Va ri ab les S e c tion ----------------- units lj b ound a ry f p p a tom _ style a tomi c c omm _ modify vel yes newton on region left B b lo c k { xlo } { b uffer _ size } { ylo } { yhi } { zlo } { zhi } region right B b lo c k { xhi - b uffer _ size } { xhi } { ylo } { yhi } { zlo } { zhi } region leftshe a r b lo c k { xlo } { she a r _ size } { ylo } { yhi } { zlo } { zhi } region rightshe a r b lo c k { xhi - she a r _ size } { xhi } { ylo } { yhi } { zlo } { zhi } region left B in b lo c k { offset } { b uffer _ size - offset } { offset } { yhi - offset } { offset } { zhi - offset } region right B in b lo c k { xhi - b uffer _ size + offset } { xhi - offset } { offset } { yhi - offset } { offset }{ zhi - offset } region roi b lo c k { b uffer _ size } { xhi - b uffer _ size } { ylo } { yhi } { zlo } { zhi } a tom _ style v a lue The style of atoms in the system. The o b md extension was tested using a tomi c , full , and mole c ul a r atom A Appendices The following appendices specific versions of the documents that are extracted from the OBMD GitHub repository (the latest versions can be found there). A.1 OBMD User Guide When performing OBMD simulations, a momentum-conserving thermostat should be applied. Usually, the DPD thermostat is used, where velocities should be communicated between neighboring processors and stored as properties of ghost atoms using c omm _ modify vel yes . Importantly, Newton's third law for pairwise and bonded interactions should be enabled via newton on . One should also define geometric regions. Here, 6 regions (namely left B , right B , leftshe a r , rightshe a r , left B in , and right B in ) must be specified. All regions must be defined using the following format: xlo xhi ylo yhi zlo zhi (for instructions on b lo c k style, readers are referred to the LAMMPS documentation). SIMULATION SETTINGS AND SYSTEM DEFINITION styles. left B b lo c k a rgs Region defining the left buffer, where R egion - ID = left B . b uffer _ size is the length of the buffer region, xlo , ylo , and zlo are the lower bounds of the simulation box in the x -, y -, and z -direction, respectively, while xhi , yhi , and zhi are the upper bounds. right B b lo c k a rgs Region defining the right buffer, where R egion - ID = right B . leftshe a r b lo c k a rgs Region defining part of the left buffer where shear forces are applied. Here, R egion - ID = leftshe a r . she a r _ size is the length of the region where shear flow is imposed. If only normal pressure forces are applied, e.g., in the case of equilibrium simulations, she a r _ size , xlo , xhi , ylo , yhi , zlo , and zhi should be set to 0.0 . rightshe a r b lo c k a rgs Region defining part of the right buffer where shear forces are applied. Here, R egion - ID = rightshe a r . If only normal pressure forces are applied, e.g., in the case of equilibrium simulations, she a r _ size , xlo , xhi , ylo , yhi , zlo , and zhi should be set to 0.0 . left B in b lo c k a rgs Region defining part of the left buffer where new particles are inserted. Here, R egion - ID = left B in . offset is the value used to narrow the insertion region. right B in b lo c k a rgs Region defining part of the right buffer where new particles are inserted. Here, R egion - ID = right B in . # ----------------- I nter ac tion S e c tion ----------------- p a ir _ style dpd { temp } { r c } { seed _ dpd } # ----------------- A tom D efinition S e c tion ----------------- re a d _ d a t a dpd _ 8 m a p _ o b md . d a t a # ----------------- S ettings S e c tion ----------------- p a ir _ c oeff * * { a ij } { g a mm a _ dpd } { r c } neigh b or { skin } b in neigh _ modify del a y 0 every 1 timestep { dt } fix 1 a ll nve p a ir _ style dpd is used to perform simulations employing DPD force field. According to the LAMMPS documentation it takes temperature ( temp ), cutoff of the interaction ( r c ), and seed ( seed _ dpd ) as input parameters. Alongside, the following coefficients ( p a ir _ c oeff ) must be given: interaction parameter between atoms i and j ( a ij ) friction parameter ( g a mm a _ dpd ) cutoff of the interaction ( r c ) MultiXscale Deliverable 2.6 Page 15 The above fix invokes the newly added o b md extension in LAMMPS. The associated parameters are given as follows: re a d _ d a t a is used to read a pre-equilibrated configuration of the DPD fluid given in the dpd _ 8 m a p _ o b md . d a t a file. The latter contains information about the positions of the atoms in reduced units. It is important that the position of the atoms do not exceed xlo or xhi of the simulation box; if so, they should be manually "deleted" before conducting the OBMD simulations. fix 2 a ll o b md { ntype } { nfreq } { seed } { pxx } { pxy } { pxz } { dpxx } { freq } { a lph a } { t a u } { n b uf } region 1 left B region 2 right B region 3 leftshe a r region 4 rightshe a r region 5 left B in region 6 right B in b uffersize { b uffer _ size } gf ac { gf ac } stepp a r a llel { step _ p a r a llel } stepperp { step _ perp } m a x a ttempt { m a x a ttempt } usher { usher _ fl a g } { et a rget } { ds 0 } { dthet a 0 } { uovlp } { dsolvp } { eps } { n a ttempt } c h a rged { c h a rge _ fl a g } FIX_ID Unique user-assigned name for the fix. Here, FIX _ ID = 2 . GROUP _ ID ID of the group of atoms to which fix o b md is applied. The group used here must already exist. Here, GROUP _ ID = a ll . FIX _ STYLE Name of a fix style to be used. In the case of OBMD simulations, one must choose o b md . ntype Type of inserted particles. Commonly, the ntype is chosen to correspond to the type of particles that are part of the GROUP _ ID group, while for molecules, it should be set to 0 . nfreq Specifies the number of timesteps after which the fix o b md is executed. Typically it set to 1 , meaning that the o b md extension is invoked at every timestep. seed A random number used to generate the initial position of newly inserted particles. pxx Pressure applied onto the buffer regions. It represents the equilibrium pressure of a fluid and can be derived from the equation of state, or obtained from an NVT simulation by monitoring the pressure to which the system relaxes. Can be provided as a real constant number or as a variable using v _ pxx . pxy Shear stress with forces along the y -axis. It can be given as a real constant number or as a variable using v _ pxy . pxz Shear stress with forces along the z -direction. It can be given as a real constant number or as a variable using v _ pxz . dpxx v a lue Pressure amplitude added to the variable pxx and applied only onto the left buffer. It can be utilized to simulate a mechanical pressure wave with a frequency freq and can be provided as a real constant number or as a variable using v _ dpxx . freq Frequency of the mechanical pressure wave (if excited, otherwise should be set to 0.0 ). It can be given as a real constant number or as a variable using v _ freq . a lph a Parameter to (further) reduce the desired number of particles in the buffers. t a u Characteristic relaxation time of the buffers. n b uf Desired number of particles in the buffer. It can be computed using nbuf = N ∗ lb / lx , where N is the total number of particles in the simulation box, lb is the buffer length, and lx is the length of the simulation box. MultiXscale Deliverable 2.6 Page 16 SIMULATION EXECUTION region 1 v a lue Region ID of the left buffer, where normal forces are applied. Here, region ID is v a lue = left B . region 2 v a lue Region ID of the right buffer, where normal forces are applied. Here, region ID is v a lue = right B . region 3 v a lue Region ID of the left buffer, where tangential shear forces are applied. Here, region ID is v a lue = leftshe a r . region 4 v a lue R egion ID of the right b uffer , where t a ngenti a l she a r for c es a re a pplied . H ere , value = rightshear`. region 5 v a lue Region ID of the left buffer, where new particles will be inserted. Here, v a lue = left B in . region 6 v a lue Region ID of the right buffer, where new particles will be inserted. Here, v a lue = right B in . b uffersize v a lue Length of the buffer region. Here, v a lue = b uffer _ size . gf ac v a lue Smoothing length. Its value must be between 0.0 and 1.0 . Here, v a lue = gf ac . stepp a r a llel v a lue Weighting function used to distribute the total external force acting in normal direction. For now only "smooth" weighting function is implemented and v a lue = 1 must be used. Here, v a lue = step _ p a r a llel = 1 . stepperp v a lue Weighting function used to distribute the total external force acting in tangential direction. For now only Heaviside step weighting function is implemented and v a lue = 0 must be used. Here, v a lue = step _ perp = 0 . m a x a ttempt v a lue The maximum number of attempts to perform the insertion by calling the try _ inserting () member function. Here, v a lue = m a x _ a ttempt . If m a x a ttempt v a lue is not given, the default value is set to 1 . usher options I nsertion a lgorithm , where the options` are: usher _ fl a g usher _ fl a g must be set to 1 if one aims to employ the USHER insertion algorithm. et a rget Target energy of the newly inserted particle. ds 0 Displacement step used by the USHER algorithm, which should be small but not too small compared to the interaction range. dthet a 0 Angular steps used when inserting a molecule. uovlp Very large energy characterizing the overlap position, which is commonly set to 10000 . dsolvp Parameter set to the value of the first RDF maximum. eps Parameter used in case of overlaps. n a ttempt Maximum number of iterations performed by the USHER algorithm, which is typically set to 40 . # ----------------- C ompute S e c tion ----------------- # # ----------------- O utput S e c tion ----------------- run { steps } MultiXscale Deliverable 2.6 Page 17 The run option is used to execute the simulation, where steps represent the number of simulation steps to be performed. Finally, the in . simul a tion is written using file _ to _ write . write ( c ontent _ 4 in ) . COMBINED INPUT SCRIPT TO RUN OBMD SIMULATIONS USING LAMMPS Below are example values to generate a file for the OBMD simulation. "" " file _ to _ write . write ( c ontent _ 4 in ) file _ to _ write . c lose () # -------------------------------------------------------- # # write input s c ript if __ n a me __ = "__ m a in __" : write _ in ( in _ file ) print ( f " file { in _ file } is written " ) # ----------------- Va ri ab les S e c tion ----------------- units lj b ound a ry f p p a tom _ style a tomi c c omm _ modify vel yes newton on region left B b lo c k 0 . 0 5 . 0391 0 . 0 11 . 198 0 . 0 11 . 198 region right B b lo c k 28 . 5549 33 . 594 0 . 0 11 . 198 0 . 0 11 . 198 region leftshe a r b lo c k 0 . 0 0 . 0 0 . 0 0 . 0 0 . 0 0 . 0 region rightshe a r b lo c k 0 . 0 0 . 0 0 . 0 0 . 0 0 . 0 0 . 0 region left B in b lo c k 0 . 0 5 . 0391 0 . 0 11 . 198 0 . 0 11 . 198 region right B in b lo c k 28 . 5549 33 . 594 0 . 0 11 . 198 0 . 0 11 . 198 region roi b lo c k 5 . 0391 28 . 5549 0 . 0 11 . 198 0 . 0 11 . 198 # ----------------- I nter ac tion S e c tion ----------------- p a ir _ style dpd 1 . 0 1 . 0 2616 # ----------------- A tom D efinition S e c tion ----------------- re a d _ d a t a dpd _ 8 m a p _ o b md . d a t a # ----------------- S ettings S e c tion ----------------- p a ir _ c oeff * * 209 . 6 4 . 5 1 . 0 neigh b or 0 . 4 b in neigh _ modify del a y 0 every 1 timestep 0 . 001464 fix 1 a ll nve fix 2 a ll o b md 1 1 6111 188 . 0 0 . 0 0 . 0 0 . 0 0 . 0 0 . 7 0 . 005 1327 & region 1 left B region 2 right B region 3 leftshe a r & region 4 rightshe a r region 5 left B in region 6 right B in & b uffersize 5 . 0391 gf ac 0 . 25 stepp a r a llel 0 stepperp 1 & m a x a ttempt 1 usher 1 31 . 03 1 . 0 0 . 02 10000 1 . 5 1 . 0 40 c h a rged 0 # ----------------- C ompute S e c tion ----------------- # # ----------------- O utput S e c tion ----------------- run 2000000 MultiXscale Deliverable 2.6 Page 18 If you are using fix o b md , please cite the following articles: R. Delgado-Buscalioni, J. Sablić, and M. Praprotnik, Eur. Phys. J. Special Topics 224, 2331-2349 (2015). R. Delgado-Buscalioni, K. Kremer, and M. Praprotnik, J. Chem. Phys. 128, 114110 (2008). J. Sablić, M. Praprotnik, and R. Delgado-Buscalioni, Soft Matter 12, 2416-2439 (2016). L. Delle Site and M. Praprotnik, Phys. Rep. 693, 1-56 (2017). P. Papež and M. Praprotnik, J. Chem. Theory Comput. 18, 1227-1240 (2022). R. Delgado-Buscalioni and P. V. Coveney, J. Chem. Phys. 119, 978-987 (2003). MultiXscale Deliverable 2.6 Page 19 In the Programmer's manual, instructions on how to compile LAMMPS (version 7 Feb 2024) with the o b md extension are provided, along with explanation of the structure, classes, and function of the added lj / c ut / rf pair style. The code introduced below can be downloaded here. OBMD PACKAGE The OBMD package is added to the LAMMPS code as a subdirectory with the package name in capital letters. OBMD package represents the extension that allows grand-canonical simulations in and out-of-equilibrium using particles. To this end, it implements routines that take care of the deletion and insertion of particles, as well as computation and distribution of external forces, used to impose user-defined external boundary conditions (e.g., constant normal load, shear flow, mechanical pressure wave). Build LAMMPS using CMake with OBMD extension enabled To perform OBMD simulations, LAMMPS code needs to be compiled with the OBMD package enabled. Instructions on how to compile LAMMPS can be found here, while the OBMD package can be included using - D PKG _ OBMD = on when compiling. Additions of the OBMD package FIX STYLE obmd In this section, we explain the o b md extension, which is implemented as a new fix style within the F ix Ob md M erged class and is used for conducting OBMD simulations both in and out-of-equilibrium. This fix manages the deletion and insertion of particles while concurrently monitoring the outgoing and incoming linear momentum of both deleted and newly inserted particles, as implemented in the fix _ o b md _ merged . c pp and its corresponding fix _ o b md _ merged . h header file. Accurately measuring the outgoing and ingoing linear momentum is essential for imposing the external boundary conditions through external forces, which is also performed by the o b md fix style. The LAMMPS input syntax is as follows: The FIX _ ID serves as the unique identifier for the fix, while GROUP _ ID specifies the ID of the group to which the fix will be applied. The name of the o b md fix style is followed by 11 essential parameters that are required for its execution. Namely: #!/ b in / sh mkdir - p b uild _ mpi c d b uild _ mpi c m a ke .. / c m a ke \ - D PKG _ OBMD = on m a ke - j 12 sr c / OBMD fix o b md fix FIX _ ID GROUP _ ID o b md ntype nfreq seed pxx pxy pxz dpxx freq a lph a t a u n b uf keyword v a lue ntype Atom type to assign to the newly inserted particles (offset when molecule is to be inserted). A.2 OBMD Programmers Manual Additionally, there are several (optional) parameters ( keyword v a lue ) that can be included to optimize the performance, execution, and various factors within the simulation. Keyword: nfreq Specifies the number of timesteps after which the fix is executed. seed A random number (positive integer). pxx Pressure (real constant number or equal style variable). pxy Shear stress with forces acting along y -axis (real constant number or equal style variable). pxz Shear stress with forces acting along z -axis (real constant number or equal style variable). dpxx Pressure amplitude (real constant number or equal style variable). freq Frequency (real constant number or equal style variable). a lph a Parameter to (further) reduce the desired number of particles in the buffers (real constant number or equal style variable). t a u Relaxation time of the buffers (real constant number or equal style variable). n b uf Desired number of particles in the buffer (real constant number or equal style variable). region 1 v a lue Region ID of the left buffer, where normal forces are applied. region 2 v a lue Region ID of the right buffer, where normal forces are applied. region 3 v a lue Region ID of the left buffer, where tangential shear forces are applied. region 4 v a lue Region ID of the right buffer, where tangential shear forces are applied. region 5 v a lue Region ID of the left buffer, where new particles will be inserted. region 6 v a lue Region ID of the right buffer, where new particles will be inserted. b uffersize v a lue Length of the buffer region. gf ac v a lue Smoothing length (value between 0.0 and 1.0 ). stepp a r a llel v a lue Weighting function used to distribute the total external force acting in normal direction. v a lue = 1 represents the "smooth" weighting function. stepperp v a lue Weighting function used to distribute the total external force acting in tangential direction. v a lue = 0 represents the Heaviside step weighting function. m a x a ttempt v a lue The maximum number of attempts to perform the insertion by calling the try _ inserting () member function. usher v a lue To select usher as the insertion algorithm, the v a lue should be set to 1 , and followed by the parameters below. MultiXscale Deliverable 2.6 Page 21 They can be set to EQUAL or ATOM , depending whether the variable is of equ a l or a tom style, and are used in the pre _ ex c h a nge () method. If " v _" is not given, the arguments are treated as numeric constant values. In the init () method, checks are performed to ensure that the regions required for conducting OBMD simulations are defined, and error is raised if they are not. During the OBMD simulation, particles will be deleted from and inserted into the system. Ensuring that the IDs of particles remain unique, the m a xt a g _ a ll (representing current maximal atom ID) or/and m a xmol _ a ll (representing current maximal molecule ID), is/are also determined (and distributed across all processes). t a uv a r = input - > v a ri ab le - > find ( t a ustr ); if ( t a uv a r < 0 ) error - > a ll ( FLERR , " Va ri ab le {} for fix o b md does not exist ", t a ustr ); if ( input - > v a ri ab le - > equ a lstyle ( t a uv a r )) t a ustyle = EQUAL ; else if ( input - > v a ri ab le - > a tomstyle ( t a uv a r )) t a ustyle = ATOM ; else error - > a ll ( FLERR , " Va ri ab le {} for fix o b md is inv a lid ", t a ustr ); } if ( n b ufstr ) { n b ufv a r = input - > v a ri ab le - > find ( n b ufstr ); if ( n b ufv a r < 0 ) error - > a ll ( FLERR , " Va ri ab le {} for fix o b md does not exist ", n b ufstr ); if ( input - > v a ri ab le - > equ a lstyle ( n b ufv a r )) n b ufstyle = EQUAL ; else if ( input - > v a ri ab le - > a tomstyle ( n b ufv a r )) n b ufstyle = ATOM ; else error - > a ll ( FLERR , " Va ri ab le {} for fix o b md is inv a lid ", n b ufstr ); } if ( input - > v a ri ab le - > equ a lstyle ( pxxv a r )) pxxstyle = EQUAL ; else if ( input - > v a ri ab le - > a tomstyle ( pxxv a r )) pxxstyle = ATOM ; pxx = utils :: numeri c ( FLERR , a rg [ 6 ], f a lse , lmp ); pxxstyle = CONSTANT ; // get regions iregion = dom a in - > get _ region _ b y _ id ( idregion ); // LEFT BUFFER if (! iregion ) error - > a ll ( FLERR , " R egion ID for fix o b md does not exist "); iregion 2 = dom a in - > get _ region _ b y _ id ( idregion 2 ); // RIGHT BUFFER if (! iregion 2 ) error - > a ll ( FLERR , " R egion ID for fix o b md does not exist "); iregion 3 = dom a in - > get _ region _ b y _ id ( idregion 3 ); // SHEAR FLOW LEFT if (! iregion 3 ) error - > a ll ( FLERR , " R egion ID for fix o b md does not exist "); iregion 4 = dom a in - > get _ region _ b y _ id ( idregion 4 ); // SHEAR FLOW RIGHT if (! iregion 4 ) error - > a ll ( FLERR , " R egion ID for fix o b md does not exist "); iregion 5 = dom a in - > get _ region _ b y _ id ( idregion 5 ); // INSERTION LEFT if (! iregion 5 ) error - > a ll ( FLERR , " R egion ID for fix o b md does not exist "); iregion 6 = dom a in - > get _ region _ b y _ id ( idregion 6 ); // INSERTION RIGHT if (! iregion 6 ) error - > a ll ( FLERR , " R egion ID for fix o b md does not exist "); MultiXscale Deliverable 2.6 Page 28 The allocated memory is freed in the destructor that follows the constructor. The execution of the o b md fix during the timestep is controlled by the setm a sk () method. The latter invokes methods of the F ix Ob md M erged class at the pre _ ex c h a nge () and post _ for c e () stages of the timestep. The deletion and insertion of particles should be performed before the neighbor list is rebuilt, therefore, these routines are implemented in the pre _ ex c h a nge () method. On the other hand, external forces, which impose the external boundary conditions, should be applied after evaluating the bonded and non-bonded forces. Therefore, these forces are applied to the particles in the buffers in the post _ for c e () method. // find c urrent m a x a tom a nd mole c ule ID s if ne c ess a ry if ( idnext ) find _ m a xid (); /* ---------------------------------------------------------------------- */ void F ix Ob md M erged :: find _ m a xid () { t a gint * t a g = a tom - > t a g ; t a gint * mole c ule = a tom - > mole c ule ; int nlo ca l = a tom - > nlo ca l ; t a gint m a x = 0 ; for ( int i = 0 ; i < nlo ca l ; i ++ ) m a x = MAX ( m a x , t a g [ i ]); MPI _ A llredu c e ( & m a x , & m a xt a g _ a ll , 1 , MPI _ LMP _ TAGINT , MPI _ MAX , world ); if ( mode == MOLECULE && mole c ule ) { m a x = 0 ; for ( int i = 0 ; i < nlo ca l ; i ++ ) m a x = MAX ( m a x , mole c ule [ i ]); MPI _ A llredu c e ( & m a x , & m a xmol _ a ll , 1 , MPI _ LMP _ TAGINT , MPI _ MAX , world ); } } /* ---------------------------------------------------------------------- */ F ix Ob md M erged :: ~ F ix Ob md M erged () { delete r a ndom ; delete [] molfr ac ; delete [] idrigid ; delete [] idsh a ke ; delete [] idregion ; delete [] idregion 2 ; delete [] idregion 3 ; delete [] idregion 4 ; delete [] idregion 5 ; delete [] idregion 6 ; memory - > destroy ( c oords ); memory - > destroy ( im a gefl a gs ); memory - > destroy ( list ); memory - > destroy ( m a rk ); } /* ---------------------------------------------------------------------- */ int F ix Ob md M erged :: setm a sk () { int m a sk = 0 ; m a sk | = PRE _ EXCHANGE ; m a sk | = POST _ FORCE ; return m a sk ; } MultiXscale Deliverable 2.6 Page 29 Invoking the pre _ ex c h a nge () method, call to the try _ deleting () member function is made first. This function takes R egion _ id , vnewl , and vnewr as arguments. Here, R egion - ID is set to iregion when considering the left buffer, and to iregion 2 when considering the right buffer. In the try _ deleting () member function, a loop over all nlo ca l particles is made. If the particle's position in the open direction (i.e., x [][ 0 ] , corresponding to the position in the x -direction) is less than the lower limit ( b oxl ) or greater than the upper limit ( b oxh ) of the simulation box in the x -direction, the particle's local index ( i ) is added to the list of particles that need to be deleted ( list [ n c ount ++ ] = i ). Each process stores the number of particles to be deleted in variable n c ount . The total number of particles to be deleted is stored in value n a ll , representing the summed value of n c ount from all processes. The sum of n c ount from all processes with lower rank (i.e., the processes before the current one) is stored in n b efore . After, the deletion is performed. If mode == ATOM , the index of a particle to be deleted is randomly chosen across processors and marked for deletion ( m a rk [] = 1 ). Variables n c ount , ndel , and n a ll are updated accordingly. Using the post-decrement operator, the values of n c ount and n a ll are decreased by 1 , while using post-increment operator, the value of variable del is increased by 1 . The latter stores the total number of deletions performed. /* ---------------------------------------------------------------------- perform p a rti c le / mole c ule insertions / deletions ------------------------------------------------------------------------- */ void F ix Ob md M erged :: pre _ ex c h a nge () { // delete p a rti c les th a t c ross the open b ound a ries a fter first h a lf of velo c ity - V erlet a lgorithm try _ deleting ( iregion , vnewl , vnewr ); try _ deleting ( iregion 2 , vnewl , vnewr ); int n c ount = 0 ; for ( i = 0 ; i < nlo ca l ; i ++ ) { if ( x [ i ][ 0 ] < b oxl || x [ i ][ 0 ] > b oxh ) { std :: c out << " D eleting x [ i ][ 0 ] = " << x [ i ][ 0 ] << " i = " << i << " type [ i ] = " << type [ i ] << " glo ba l id = " << a tom - > t a g [ i ] << std :: endl ; list [ n c ount ++ ] = i ; } } int n a ll , n b efore ; MPI _ A llredu c e ( & n c ount , & n a ll , 1 , MPI _ INT , MPI _ SUM , world ); MPI _ Sca n ( & n c ount , & n b efore , 1 , MPI _ INT , MPI _ SUM , world ); n b efore - = n c ount ; // a tomi c deletions // c hoose a toms r a ndomly ac ross a ll pro c s a nd m a rk them for deletion // shrink eligi b le list a s my a toms get m a rked // keep ndel , n c ount , n a ll , n b efore c urrent a fter e ac h a tom deletion if ( mode == ATOM ) { while ( n a ll ) { iwhi c hglo ba l = st a ti c _ ca st < int > ( n a ll * r a ndom - > uniform ()); if ( iwhi c hglo ba l < n b efore ) n b efore --; else if ( iwhi c hglo ba l < n b efore + n c ount ) { iwhi c hlo ca l = iwhi c hglo ba l - n b efore ; m a rk [ list [ iwhi c hlo ca l ]] = 1 ; list [ iwhi c hlo ca l ] = list [ n c ount - 1 ]; n c ount --; } ndel ++ ; n a ll --; MultiXscale Deliverable 2.6 Page 30 If a molecule is to be deleted, a similar routine is executed, followed by the deletion of bonds, angles, proper dihedral angles, and improper dihedral angles (if any). The vital part of the implemented OBMD method is monitoring the outgoing linear momentum. The positions ( x [][ 0 ] ) of all nlo ca l particles marked for deletion are extracted and checked to determine if they correspond to the crossing of the outer boundary of the left or right buffer region. If the particle has exited the left buffer, the variable vnewl is incremented by its linear momentum, while vnewr is incremented if the particle has exited the right buffer. The deletion of particles in the pre _ ex c h a nge () method is followed by the computation of the number of particles to be inserted. The number of particles in the left and right buffer is stored in variables c nt _ left and c nt _ right , respectively, which are computed utilizing the group - > c ount ( GROUP _ ID , R egion - ID ) implemented in the LAMMPS code. Equation Δ Nb = δt τb ( ⟨ αNb ⟩ − Nb ) gives the number of particles to be inserted into the left and right buffer, which are stored in variables ninsert _ left and ninsert _ right , respectively. By calling try _ inserting () member function (for each buffer region separately), the insertion procedure starts. The try _ inserting () method takes R egion - ID , ninsert _ left / right , vnewl , and vnewr as arguments. The try _ inserting () member function supports both ne a r and usher insertion algorithms, however, here we focus on the latter. Readers interested in the former are referred to the LAMMPS documentation. } } for ( i = nlo ca l - 1 ; i >= 0 ; i --) { if ( m a rk [ i ]) { if ( x [ i ][ 0 ] < 0 . 5 * ( b oxh + b oxl )) { vnewl [ 0 ] += m a ss [ type [ i ]] * vel [ i ][ 0 ]; vnewl [ 1 ] += m a ss [ type [ i ]] * vel [ i ][ 1 ]; vnewl [ 2 ] += m a ss [ type [ i ]] * vel [ i ][ 2 ]; } else { vnewr [ 0 ] += m a ss [ type [ i ]] * vel [ i ][ 0 ]; vnewr [ 1 ] += m a ss [ type [ i ]] * vel [ i ][ 1 ]; vnewr [ 2 ] += m a ss [ type [ i ]] * vel [ i ][ 2 ]; } a ve c - > c opy ( a tom - > nlo ca l - 1 , i , 1 ); a tom - > nlo ca l --; } } // c ount num b er of p a rti c les in left / right b uffer c nt _ left = group - > c ount ( igroup , iregion ); c nt _ right = group - > c ount ( igroup , iregion 2 ); // ca l c ul a te num b er of p a rti c les needed for insertion ninsert _ left = - st a ti c _ ca st < int > (( st a ti c _ ca st < dou b le > ( c nt _ left ) / mol _ len - a lph a * n b uf ) * upd a te - > dt / t a u ); ninsert _ right = - st a ti c _ ca st < int > (( st a ti c _ ca st < dou b le > ( c nt _ right ) / mol _ len - a lph a * n b uf ) * upd a te - > dt / t a u ); // tries inserting ninsert _ left / right p a rti c les into left / right b uffer try _ inserting ( iregion 5 , ninsert _ left , vnewl , vnewr ); try _ inserting ( iregion 6 , ninsert _ right , vnewl , vnewr ); MultiXscale Deliverable 2.6 Page 31 After selecting a random position within the buffer region for the particle to be inserted, its (potential) energy is computed by calling the usher () member function and stored in the variable entmp . In the usher () method, call to energy _ a tomisti c _ o b md () or energy () member function is made. The former is used when simulating charged particles and takes the following parameters: When simulating non-charged particles, the energy () member function is called and it takes as arguments. In the energy _ a tomisti c _ o b md () member function, the downcast conversion of the pointer to the pair potential class Pa ir LJC ut RF is made, so that newly added single _ a tomisti c _ o b md () method implemented in the Pa ir LJC ut RF class can be accessed. else if ( usherfl a g ) { me = - 1 ; entmp = usher ( iregion _ v a r , c oords , et a rget , n a tom , imol , iter ); if ( entmp < et a rget + EPSILON ) { if ( c omm - > me == 0 ) std :: c out << " USHER acc epts a t E = " << entmp << " in a ttempt N o . " << a ttempt << " with " << iter << " iter a tions " << std :: endl ; } else { if ( c omm - > me == 0 ) std :: c out << " USHER denies a t E = " << entmp << " a t a ttempt N o . " << a ttempt << std :: endl ; fl a g = 1 ; } } iregion _ v a r Region ID where particle to be inserted is located. onemols []- > q [] Charge of the particle to be inserted. type _ temp Type of the particle to be inserted. c oords [] Randomly selected position of the particle to be inserted (vector of length three). fusher Force acting on the particle to be inserted. index (Not relevant when simulating DPD particles.) type _ temp c oords [] fusher if ( c h a rgefl a g ) { entmp += energy _ a tomisti c _ o b md ( iregion _ v a r , onemols [ imol ]- > q [ m ], type _ temp , c oords [ m ], fusher ); } else { entmp += energy ( 1 , type _ temp , c oords [ m ], fusher ); // i - does not m a tter , itype inserted , c oords [ m ], fusher } // potenti a l ac ting b etween a tomisti c p a rti c les a uto p a ir = dyn a mi c _ ca st < Pa ir LJC ut RF * > ( for c e - > p a ir _ m a t c h (" lj / c ut / rf ", 1 )); MultiXscale Deliverable 2.6 Page 32 Looping over all nlo ca l particles, the distance between the particle to be inserted and the j -th particle of nlo ca l particles is calculated. If the distance ( rsq ) is less than a cutoff ( c utsq ), the energy and the force acting on the particle to be inserted are computed. The energy _ a tomisti c _ o b md () member function returns the total energy of the particle to be inserted, which is stored in variable tot a l _ energy . Similar is performed when a call to the energy () function is made. The computed values of the total energy and force for each process are stored in variables entmp and fusher , respectively. The summed values of energy and force over all processes are stored in variables entmp _ a ll and fusher _ a ll , respectively. The computed total energy ( entmp _ a ll ) is first compared to the predefined threshold, and if it is lower, the position is accepted and the usher () member function is exited. Otherwise, the position is iteratively corrected according to the following update rule: r n +1 = r n + f n | f n | δsn , where displacement δsn depends on the local potential energy (stored in variable entmp _ a ll ) and force magnitude (stored in variable f ab s ). If the computed energy is larger than uovlp , which is chosen to be a very large energy representing the overlap position, δsn is computed using δsn = Δ sovlp = rσ − ( 4 ε Un ) 1 2, where particle position will be translated for a distance rσ away from the center of the overlapped particle. for ( int j = 0 ; j < nlo ca l ; j ++ ) { delx = c oord [ 0 ] - x [ j ][ 0 ]; dely = c oord [ 1 ] - x [ j ][ 1 ]; delz = c oord [ 2 ] - x [ j ][ 2 ]; dom a in - > minimum _ im a ge ( delx , dely , delz ); rsq = delx * delx + dely * dely + delz * delz ; jtype = type [ j ]; if ( rsq < c utsq [ itype ][ jtype ]) { tot a l _ energy += p a ir - > single _ a tomisti c _ o b md ( qi , j , itype , jtype , rsq , f ac tor _ c oul , f ac tor _ lj , fp a ir ); fusher [ 0 ] += fp a ir * delx ; fusher [ 1 ] += fp a ir * dely ; fusher [ 2 ] += fp a ir * delz ; } } MPI _ A llredu c e ( fusher , fusher _ a ll , 3 , MPI _ DOUBLE , MPI _ SUM , world ); MPI _ A llredu c e ( & entmp , & entmp _ a ll , 1 , MPI _ DOUBLE , MPI _ SUM , world ); if ( entmp _ a ll < et a rget + EPSILON ) b re a k ; else if ( entmp _ a ll > uovlp ) { f ab s = sqrt ( fusher _ a ll [ 0 ] * fusher _ a ll [ 0 ] + fusher _ a ll [ 1 ] * fusher _ a ll [ 1 ] + fusher _ a ll [ 2 ] * fusher _ a ll [ 2 ]); if ( f ab s < EPSILON ) c ontinue ; ds = dsovlp - pow ( 4 * eps / entmp _ a ll , 1 . 0 / 12 . 0 ); // loop over a toms of mole c ule MultiXscale Deliverable 2.6 Page 33 When the computed energy is lower than uovlp , δsn is determined using δsn = min ( Δ s , Un − U 0 | f n | ) . Here, U 0 represents the target energy, which is stored in variable et a rget , while Δ s stands for the displacement step, which is stored in variable ds 0 . In case of molecule insertion, rotation of the molecule about its center of mass is also performed. The position obtained by the iterative procedure is verified to remain within the buffer region by calling the c he c k _ mol _ region () member function. If the iterative correction of the position leads to particle leaving the buffer region, the usher () method is exited, and the insertion procedure is terminated. On the other hand, if the iterative procedure is successful, a new particle/molecule is inserted. In case of molecule insertion, bonds, angles, proper dihedrals, and improper dihedrals (if any) are also created. As the newly created particles are inserted with zero velocity, the incoming linear momentum is also zero. Therefore, the computed change in the linear momentum comes only from the particles that have exited the outer boundary of the left and right buffer region. Values for the left and right buffer are for each process stored in variables vnewl and vnewr , respectively, while the corresponding summed values from all processes are stored in variables vnewl _ a ll and vnewr _ a ll , respectively. for ( m = 0 ; m < n a tom ; m ++ ) { c oords [ m ][ 0 ] += fusher _ a ll [ 0 ] / f ab s * ds ; c oords [ m ][ 1 ] += fusher _ a ll [ 1 ] / f ab s * ds ; c oords [ m ][ 2 ] += fusher _ a ll [ 2 ] / f ab s * ds ; } int c he c k = c he c k _ mol _ region ( iregion _ v a r , c oords , n a tom ); if ( c he c k == 1 ) b re a k ; // exit - p a rti c le is now within the b uffer region } else { f ab s = sqrt ( fusher _ a ll [ 0 ] * fusher _ a ll [ 0 ] + fusher _ a ll [ 1 ] * fusher _ a ll [ 1 ] + fusher _ a ll [ 2 ] * fusher _ a ll [ 2 ]); if ( f ab s < EPSILON ) c ontinue ; ds = std :: min (( entmp _ a ll - et a rget ) / f ab s , ds 0 ); if ( mode == MOLECULE ) { torq ab s = sqrt ( torq _ a ll [ 0 ] * torq _ a ll [ 0 ] + torq _ a ll [ 1 ] * torq _ a ll [ 1 ] + torq _ a ll [ 2 ] * torq _ a ll [ 2 ]); dthet a = std :: min (( entmp _ a ll - et a rget ) / torq ab s , dthet a 0 ); Ma th E xtr a :: norm 3 ( torq _ a ll ); Ma th E xtr a :: a xis a ngle _ to _ qu a t ( torq _ a ll , dthet a , qu a t ); Ma th E xtr a :: qu a t _ to _ m a t ( qu a t , rotm a t ); } for ( m = 0 ; m < n a tom ; m ++ ) { c oords [ m ][ 0 ] += fusher _ a ll [ 0 ] / f ab s * ds ; c oords [ m ][ 1 ] += fusher _ a ll [ 1 ] / f ab s * ds ; c oords [ m ][ 2 ] += fusher _ a ll [ 2 ] / f ab s * ds ; if ( mode == MOLECULE ) { Ma th E xtr a :: m a tve c ( rotm a t , c oords [ m ], c oords _ tmp ); Ma th E xtr a :: c opy 3 ( c oords _ tmp , c oords [ m ]); } } int c he c k = c he c k _ mol _ region ( iregion _ v a r , c oords , n a tom ); if ( c he c k == 1 ) b re a k ; // exit } MPI _ A llredu c e ( vnewl , vnewl _ a ll , 3 , MPI _ DOUBLE , MPI _ SUM , world ); MPI _ A llredu c e ( vnewr , vnewr _ a ll , 3 , MPI _ DOUBLE , MPI _ SUM , world ); MultiXscale Deliverable 2.6 Page 34 The insertion procedure is followed by the computation of the external forces: F ext = ∑ i ∈ B f ext i = A J p ⋅ n − ∑ i ′∈ B Δ( mi ′ vi ′) δt , which are utilized in the post _ for c e () method. The first term on the right-hand side of the above equation reads as: A J p ⋅ n = pxxA n . when a system under a constant normal load applied from both sides of the simulation domain is simulated. To introduce a sinusoidal disturbance, it should be rewritten into A J p ⋅ n = pxxA n + Δ p sin(2 πνt ) A n while it is expressed as: A J p ⋅ n = pxxA n + pxyA t when a shear flow is introduced. Variables momentum F or c e []_ left / right and she a r F or c e _ left / right [] represent the array of length three and carry information about the external boundary condition of the constant normal load (or mechanical pressure wave) and shear stress to be applied in the x -, y -, and z -directions, respectively. The o b md fix is also executed at the post _ for c e () stage of the timestep, as defined in the setm a sk () method. At this stage, external forces computed in the pre _ ex c h a nge () method are applied to the particles in the buffers. Forces acting in normal direction are distributed first by calling the reg _ for c e () member function, which takes R egion - ID , momentum F or c e _ left / right , and step _ p a r a llel as arguments. Here, R egion - ID denotes iregion corresponding to the left buffer region, where the total external force stored in variable momentum F or c e _ left will be applied. The same applies to the imposition of the total external force ( momentum F or c e _ right ) to the right buffer. The step _ p a r a llel parameter defines the weighting function used to distribute the total external force among particles within the buffers. // a re a of the b uffer - ROI interf ac e dou b le a re a = ly * lz ; simul a tion _ time += upd a te - > dt ; dou b le f ac tor = pxx + dpxx * sin ( 2 . 0 * MY _ PI * freq * simul a tion _ time ); // ca l c ul a tes momentum for c es on the left b uffer Ma th E xtr a :: zero 3 ( momentum F or c e _ left ); momentum F or c e _ left [ 0 ] = vnewl _ a ll [ 0 ] / ( upd a te - > dt ) + f ac tor * a re a ; momentum F or c e _ left [ 1 ] = vnewl _ a ll [ 1 ] / ( upd a te - > dt ); momentum F or c e _ left [ 2 ] = vnewl _ a ll [ 2 ] / ( upd a te - > dt ); she a r F or c e _ left [ 0 ] = 0 . 0 ; she a r F or c e _ left [ 1 ] = pxy * a re a ; she a r F or c e _ left [ 2 ] = pxz * a re a ; // ca l c ul a tes momentum for c es on the right b uffer Ma th E xtr a :: zero 3 ( momentum F or c e _ right ); momentum F or c e _ right [ 0 ] = vnewr _ a ll [ 0 ] / ( upd a te - > dt ) - pxx * a re a ; // c onst a nt norm a l lo a d momentum F or c e _ right [ 1 ] = vnewr _ a ll [ 1 ] / ( upd a te - > dt ); momentum F or c e _ right [ 2 ] = vnewr _ a ll [ 2 ] / ( upd a te - > dt ); she a r F or c e _ right [ 0 ] = 0 . 0 ; she a r F or c e _ right [ 1 ] = - pxy * a re a ; she a r F or c e _ right [ 2 ] = - pxz * a re a ; /* ---------------------------------------------------------------------- */ void F ix Ob md M erged :: post _ for c e ( int vfl a g ) { if ( upd a te - > ntimestep % nevery ) { return ; } MultiXscale Deliverable 2.6 Page 35 In the reg _ for c e () method the total mass of particles in the corresponding buffer is computed by calling the g _ p a r _ glo ba l _ c h a rged () member function. The computed value is stored in variable gtmp . For now, only smooth weighting function is enabled (i.e., step _ p a r a llel == 0 ), where parameter gf ac is used to define the region in which the position dependent sigmoidal function is multiplied with the mass of a particle. // p a r a llel for c es reg _ for c e ( vfl a g , iregion , momentum F or c e _ left , step _ p a r a llel ); reg _ for c e ( vfl a g , iregion 2 , momentum F or c e _ right , step _ p a r a llel ); /* ---------------------------------------------------------------------- */ void F ix Ob md M erged :: reg _ for c e ( int vfl a g , R egion * region , dou b le * momentum F or c e , int step ) { if ( region ) region - > prem a t c h (); v _ init ( vfl a g ); dou b le ** x = a tom - > x ; dou b le ** f = a tom - > f ; dou b le * m a ss = a tom - > m a ss ; int * m a sk = a tom - > m a sk ; im a geint * im a ge = a tom - > im a ge ; dou b le v [ 6 ]; int nlo ca l = a tom - > nlo ca l ; dou b le m a ss _ tmp ; int * type = a tom - > type ; dou b le gtmp = g _ p a r _ glo ba l _ c h a rged ( region , step ); dou b le glo c tmp ; dou b le unwr a p [ 3 ]; dou b le xv a l , yv a l , zv a l ; /* ---------------------------------------------------------------------- */ dou b le F ix Ob md M erged :: g _ p a r _ glo ba l _ c h a rged ( R egion * region , int step ) { if ( region ) region - > prem a t c h (); int nlo ca l = a tom - > nlo ca l ; t a gint * mole c ule = a tom - > mole c ule ; int * rep _ a tom = a tom - > rep _ a tom ; dou b le ** c ms = a tom - > c ms _ mol ; dou b le * m a ss = a tom - > m a ss ; int * type = a tom - > type ; dou b le ** x = a tom - > x ; int imol ; int * t a g = a tom - > t a g ; dou b le ** f = a tom - > f ; dou b le upper _ x = dom a in - > b oxhi [ 0 ]; // upper dou b le lower _ x = dom a in - > b oxlo [ 0 ]; // lower // init for distri b ution fun c tion dou b le g _ p a r _ a ll = 0 . 0 ; dou b le g _ p a r _ a ll _ tmp = 0 . 0 ; // for now only smooth tr a nsition - > ADD for step = 1 dou b le ca rg ; if (! step ) { for ( int i = 0 ; i < nlo ca l ; i ++ ) { // m a ss of a p a rti c le MultiXscale Deliverable 2.6 Page 36 The total mass computed by each process is stored in variable g _ p a r _ a ll , while the value returned by the g _ p a r _ glo ba l _ c h a rged () method is stored in variable g _ p a r _ a ll _ tmp , which represents the summed value of g _ p a r _ a ll from all processes. In the reg _ for c e () member function, the mass (or weighted mass) of a particle in the buffer region (stored in variable glo c tmp ) is divided by the total mass of the buffer region ( gtmp ), which is computed by calling the g _ p a r _ glo ba l _ c h a rged () method. Multiplying this ratio by the momentum F or c e _ left / right yields the force acting on the particle within the corresponding buffer region. The latter is applied to each particle within the buffer region. dou b le m a ss _ temp = m a ss [ type [ i ]]; // ( looping over nlo ca ls ) c he c k if in pres c ri b ed region , c ontinue if not if ( region && ! region - > m a t c h ( x [ i ][ 0 ], x [ i ][ 1 ], x [ i ][ 2 ])) c ontinue ; // if here ... mole c ule in b uffer if ( x [ i ][ 0 ] < lower _ x + b uffer _ size ) { // LEFT BUFFER if ( x [ i ][ 0 ] < ( lower _ x + ( 1 . 0 - g _ f ac ) * b uffer _ size )) { // only m a ss g _ p a r _ a ll += m a ss _ temp ; } else { // sigmoid a l m a ss ca rg = 1 . 0 / g _ f ac * MY _ PI * ( x [ i ][ 0 ] - b uffer _ size - lower _ x ) / (- b uffer _ size ) - MY _ PI ; g _ p a r _ a ll += 0 . 5 * ( 1 . 0 + c os ( ca rg )) * m a ss _ temp ; } } if ( x [ i ][ 0 ] > upper _ x - b uffer _ size ) { // RIGHT BUFFER if ( x [ i ][ 0 ] > ( upper _ x - ( 1 . 0 - g _ f ac ) * b uffer _ size )) { g _ p a r _ a ll += m a ss _ temp ; } else { ca rg = 1 . 0 / g _ f ac * MY _ PI * ( x [ i ][ 0 ] - upper _ x + b uffer _ size ) / ( b uffer _ size ) - MY _ PI ; g _ p a r _ a ll += 0 . 5 * ( 1 . 0 + c os ( ca rg )) * m a ss _ temp ; } } // OTHER IS NOT OF MY INTEREST ... ROI } } else { error - > a ll ( FLERR , " F or now implemented only for ! step "); } // sum a nd distri b ute MPI _ A llredu c e ( & g _ p a r _ a ll , & g _ p a r _ a ll _ tmp , 1 , MPI _ DOUBLE , MPI _ SUM , world ); return g _ p a r _ a ll _ tmp ; } for ( int i = 0 ; i < nlo ca l ; i ++ ) { if ( region && ! region - > m a t c h ( x [ i ][ 0 ], x [ i ][ 1 ], x [ i ][ 2 ])) c ontinue ; dom a in - > unm a p ( x [ i ], im a ge [ i ], unwr a p ); m a ss _ tmp = m a ss [ type [ i ]]; if ( mode == MOLECULE ) m a ss _ tmp = mtot ; glo c tmp = g _ p a r _ lo ca l _ c h a rged ( m a ss [ type [ i ]], region , x [ i ][ 0 ], step ); xv a l = momentum F or c e [ 0 ] * glo c tmp / gtmp ; yv a l = momentum F or c e [ 1 ] * glo c tmp / gtmp ; zv a l = momentum F or c e [ 2 ] * glo c tmp / gtmp ; f [ i ][ 0 ] += xv a l ; f [ i ][ 1 ] += yv a l ; MultiXscale Deliverable 2.6 Page 37 The a llo ca te () member function handles memory allocation. While the arguments of the p a ir _ style lj / c ut / rf command are processes in the settings () method, the arguments provided to the p a ir _ c oeff command are handled in the c oeff () member function. A total of 5 or 7 arguments should be provided to the p a ir _ c oeff command, where global cutoffs are set as default, i.e., c ut _ lj _ one = c ut _ lj _ glo ba l and c ut _ c oul _ one = c ut _ c oul _ glo ba l , unless the sixth and seventh parameter are specified. Therefore, if the r c _ lj and r c _ rf arguments are not provided, the global LJ and rf cutoffs specified in the p a ir _ style command are utilized. If only one cutoff is given, it is applied for both interactions. If both cutoffs are provided, they are used as the LJ and rf cutoffs for this type pair. Ultimately, the 2D arrays, namely epsilon , sigm a , c ut _ lj , c ut _ c oul , and epsilon _ rf , are created for each atom pair. /* ---------------------------------------------------------------------- a llo ca te a ll a rr a ys ------------------------------------------------------------------------- */ void Pa ir LJC ut RF :: a llo ca te () { a llo ca ted = 1 ; int np 1 = a tom - > ntypes + 1 ; memory - > c re a te ( setfl a g , np 1 , np 1 , " p a ir : setfl a g "); for ( int i = 1 ; i < np 1 ; i ++ ) for ( int j = i ; j < np 1 ; j ++ ) setfl a g [ i ][ j ] = 0 ; memory - > c re a te ( c utsq , np 1 , np 1 , " p a ir : c utsq "); memory - > c re a te ( c ut _ lj , np 1 , np 1 , " p a ir : c ut _ lj "); memory - > c re a te ( c ut _ ljsq , np 1 , np 1 , " p a ir : c ut _ ljsq "); memory - > c re a te ( c ut _ c oul , np 1 , np 1 , " p a ir : c ut _ c oul "); memory - > c re a te ( c ut _ c oulsq , np 1 , np 1 , " p a ir : c ut _ c oulsq "); memory - > c re a te ( epsilon , np 1 , np 1 , " p a ir : epsilon "); memory - > c re a te ( sigm a , np 1 , np 1 , " p a ir : sigm a "); memory - > c re a te ( lj 1 , np 1 , np 1 , " p a ir : lj 1 "); memory - > c re a te ( lj 2 , np 1 , np 1 , " p a ir : lj 2 "); memory - > c re a te ( lj 3 , np 1 , np 1 , " p a ir : lj 3 "); memory - > c re a te ( lj 4 , np 1 , np 1 , " p a ir : lj 4 "); memory - > c re a te ( offset , np 1 , np 1 , " p a ir : offset "); memory - > c re a te ( epsilon _ rf , np 1 , np 1 , " p a ir : epsilon _ rf "); } /* ---------------------------------------------------------------------- set c oeffs for one or more type p a irs ------------------------------------------------------------------------- */ void Pa ir LJC ut RF :: c oeff ( int n a rg , c h a r ** a rg ) { if ( n a rg < 5 || n a rg > 7 ) error - > a ll ( FLERR , " I n c orre c t a rgs for p a ir c oeffi c ients "); if (! a llo ca ted ) a llo ca te (); int ilo , ihi , jlo , jhi ; utils :: b ounds ( FLERR , a rg [ 0 ], 1 , a tom - > ntypes , ilo , ihi , error ); utils :: b ounds ( FLERR , a rg [ 1 ], 1 , a tom - > ntypes , jlo , jhi , error ); dou b le epsilon _ one = utils :: numeri c ( FLERR , a rg [ 2 ], f a lse , lmp ); dou b le sigm a _ one = utils :: numeri c ( FLERR , a rg [ 3 ], f a lse , lmp ); dou b le c ut _ lj _ one = c ut _ lj _ glo ba l ; dou b le c ut _ c oul _ one = c ut _ c oul _ glo ba l ; epsilon _ rf _ one = utils :: numeri c ( FLERR , a rg [ 4 ], f a lse , lmp ); if ( n a rg >= 6 ) { c ut _ c oul _ one = c ut _ lj _ one = utils :: numeri c ( FLERR , a rg [ 4 ], f a lse , lmp ); MultiXscale Deliverable 2.6 Page 44 The completeness of the potential parameters is verified in the init _ one () method. This method invokes the Pa ir :: mix _ energy () and Pa ir : mix _ dist a n c e () member functions if the setfl a g for particle pairs is 0 or if the command p a ir _ modify mix { geometri c , a rithemi c , sixthpower } is used. Additionally, the quantities c ut _ ljsq , c ut _ c oulsq , lj 1 , lj 2 , lj 3 , and lj 4 are defined and stored as 2D arrays. The offset to shift the LJ potential value to zero at the cutoff distance is computed when the p a ir _ modify shift yes command is employed. Finally, the potential parameters arrays are symmetrized. epsilon _ rf _ one = utils :: numeri c ( FLERR , a rg [ 5 ], f a lse , lmp ); } if ( n a rg == 7 ) { c ut _ lj _ one = utils :: numeri c ( FLERR , a rg [ 4 ], f a lse , lmp ); c ut _ c oul _ one = utils :: numeri c ( FLERR , a rg [ 5 ], f a lse , lmp ); epsilon _ rf _ one = utils :: numeri c ( FLERR , a rg [ 6 ], f a lse , lmp ); } int c ount = 0 ; for ( int i = ilo ; i <= ihi ; i ++ ) { for ( int j = MAX ( jlo , i ); j <= jhi ; j ++ ) { epsilon [ i ][ j ] = epsilon _ one ; sigm a [ i ][ j ] = sigm a _ one ; c ut _ lj [ i ][ j ] = c ut _ lj _ one ; c ut _ c oul [ i ][ j ] = c ut _ c oul _ one ; // re ac tion field p a r a m epsilon _ rf [ i ][ j ] = epsilon _ rf _ one ; setfl a g [ i ][ j ] = 1 ; c ount ++ ; } } if ( c ount == 0 ) error - > a ll ( FLERR , " I n c orre c t a rgs for p a ir c oeffi c ients "); } /* ---------------------------------------------------------------------- init for one type p a ir i , j a nd c orresponding j , i ------------------------------------------------------------------------- */ dou b le Pa ir LJC ut RF :: init _ one ( int i , int j ) { if ( setfl a g [ i ][ j ] == 0 ) { epsilon [ i ][ j ] = mix _ energy ( epsilon [ i ][ i ], epsilon [ j ][ j ], sigm a [ i ][ i ], sigm a [ j ][ j ]); sigm a [ i ][ j ] = mix _ dist a n c e ( sigm a [ i ][ i ], sigm a [ j ][ j ]); c ut _ lj [ i ][ j ] = mix _ dist a n c e ( c ut _ lj [ i ][ i ], c ut _ lj [ j ][ j ]); c ut _ c oul [ i ][ j ] = mix _ dist a n c e ( c ut _ c oul [ i ][ i ], c ut _ c oul [ j ][ j ]); } dou b le c ut = MAX ( c ut _ lj [ i ][ j ], c ut _ c oul [ i ][ j ]); c ut _ ljsq [ i ][ j ] = c ut _ lj [ i ][ j ] * c ut _ lj [ i ][ j ]; c ut _ c oulsq [ i ][ j ] = c ut _ c oul [ i ][ j ] * c ut _ c oul [ i ][ j ]; lj 1 [ i ][ j ] = 48 . 0 * epsilon [ i ][ j ] * pow ( sigm a [ i ][ j ], 12 . 0 ); lj 2 [ i ][ j ] = 24 . 0 * epsilon [ i ][ j ] * pow ( sigm a [ i ][ j ], 6 . 0 ); lj 3 [ i ][ j ] = 4 . 0 * epsilon [ i ][ j ] * pow ( sigm a [ i ][ j ], 12 . 0 ); lj 4 [ i ][ j ] = 4 . 0 * epsilon [ i ][ j ] * pow ( sigm a [ i ][ j ], 6 . 0 ); if ( offset _ fl a g && ( c ut _ lj [ i ][ j ] > 0 . 0 )) { dou b le r a tio = sigm a [ i ][ j ] / c ut _ lj [ i ][ j ]; offset [ i ][ j ] = 4 . 0 * epsilon [ i ][ j ] * ( pow ( r a tio , 12 . 0 ) - pow ( r a tio , 6 . 0 )); } else offset [ i ][ j ] = 0 . 0 ; MultiXscale Deliverable 2.6 Page 45 The pairwise forces and potential between atoms i and j are computed in the c ompute () method. The charge ( q ), position vector components (in x -, y -, z -directions, i.e., x [][ 0 ] , x [][ 1 ] , x [][ 2 ] , respectively), and neighbors of the central atom i are extracted by applying a for loop over all central atoms. Applying an additional for loop over neighbors j of atom i , the square of the distance between atoms i and j is computed and stored in variable rsq . The computed distance is compared with the cutoff distance, and if it is smaller, force between atoms i and j is evaluated. Focusing only on the computation of the rf contribution, the following variables are defined: c ut _ ljsq [ j ][ i ] = c ut _ ljsq [ i ][ j ]; c ut _ c oul [ j ][ i ] = c ut _ c oul [ i ][ j ]; c ut _ c oulsq [ j ][ i ] = c ut _ c oulsq [ i ][ j ]; lj 1 [ j ][ i ] = lj 1 [ i ][ j ]; lj 2 [ j ][ i ] = lj 2 [ i ][ j ]; lj 3 [ j ][ i ] = lj 3 [ i ][ j ]; lj 4 [ j ][ i ] = lj 4 [ i ][ j ]; epsilon _ rf [ j ][ i ] = epsilon _ rf [ i ][ j ]; offset [ j ][ i ] = offset [ i ][ j ]; // c ompute I , J c ontri b ution to long - r a nge t a il c orre c tion // c ount tot a l # of a toms of type I a nd J vi a A llredu c e if ( t a il _ fl a g ) { int * type = a tom - > type ; int nlo ca l = a tom - > nlo ca l ; dou b le c ount [ 2 ], a ll [ 2 ]; c ount [ 0 ] = c ount [ 1 ] = 0 . 0 ; for ( int k = 0 ; k < nlo ca l ; k ++ ) { if ( type [ k ] == i ) c ount [ 0 ] += 1 . 0 ; if ( type [ k ] == j ) c ount [ 1 ] += 1 . 0 ; } MPI _ A llredu c e ( c ount , a ll , 2 , MPI _ DOUBLE , MPI _ SUM , world ); dou b le sig 2 = sigm a [ i ][ j ] * sigm a [ i ][ j ]; dou b le sig 6 = sig 2 * sig 2 * sig 2 ; dou b le r c 3 = c ut _ lj [ i ][ j ] * c ut _ lj [ i ][ j ] * c ut _ lj [ i ][ j ]; dou b le r c 6 = r c 3 * r c 3 ; dou b le r c 9 = r c 3 * r c 6 ; et a il _ ij = 8 . 0 * MY _ PI * a ll [ 0 ] * a ll [ 1 ] * epsilon [ i ][ j ] * sig 6 * ( sig 6 - 3 . 0 * r c 6 ) / ( 9 . 0 * r c 9 ); pt a il _ ij = 16 . 0 * MY _ PI * a ll [ 0 ] * a ll [ 1 ] * epsilon [ i ][ j ] * sig 6 * ( 2 . 0 * sig 6 - 3 . 0 * r c 6 ) / ( 9 . 0 * r c 9 ); } return c ut ; } rf _ f c tr _ 0 rf _ fctr _0 = 1 r 3 c rf _ f c tr _ 1 rf _ fctr _1 = εrf − 1 rf _ f c tr _ 2 MultiXscale Deliverable 2.6 Page 46 and the for c e c oul variable is calculated using Frf ( riαjβ ) = qiαqjβ 4 πε 0 [ 1 riαjβ 3−1 rc 3 2( εrf − 1) 2 εrf + 1 ] r iαjβ . The obtained force is scaled by f ac tor _ c oul and added to the fp a irlj , defining the fp a ir variable. The fp a ir variable is further multiplied with the displacement vector components ( delx , dely , and delz ), and forces acting in x -, y -, and z -directions are assigned to atom i . If command newton _ p a ir on is specified, forces acting in the opposite directions are applied to the atom j . When energy contributions are requested (i.e., when efl a g = 1 ), the energy interaction between atoms i and j is also computed. rf _ fctr _2 = 1 + 2 εrf /* ---------------------------------------------------------------------- */ void Pa ir LJC ut RF :: c ompute ( int efl a g , int vfl a g ) { std :: c out << " c ompute " << "\ n "; int i , j , ii , jj , inum , jnum , itype , jtype ; dou b le qtmp , xtmp , ytmp , ztmp , delx , dely , delz , evdwl , e c oul , fp a ir ; dou b le rsq , r 2 inv , r 6 inv , for c e c oul , for c elj , f ac tor _ c oul , f ac tor _ lj ; int * ilist , * jlist , * numneigh , ** firstneigh ; dou b le fp a irlj , fp a ir c oul ; evdwl = e c oul = 0 . 0 ; ev _ init ( efl a g , vfl a g ); dou b le ** x = a tom - > x ; dou b le ** f = a tom - > f ; dou b le * q = a tom - > q ; int * type = a tom - > type ; int nlo ca l = a tom - > nlo ca l ; dou b le * spe c i a l _ c oul = for c e - > spe c i a l _ c oul ; dou b le * spe c i a l _ lj = for c e - > spe c i a l _ lj ; int newton _ p a ir = for c e - > newton _ p a ir ; dou b le qqrd 2 e = for c e - > qqrd 2 e ; inum = list - > inum ; ilist = list - > ilist ; numneigh = list - > numneigh ; firstneigh = list - > firstneigh ; // loop over neigh b ors of my a toms for ( ii = 0 ; ii < inum ; ii ++ ) { i = ilist [ ii ]; qtmp = q [ i ]; xtmp = x [ i ][ 0 ]; ytmp = x [ i ][ 1 ]; ztmp = x [ i ][ 2 ]; itype = type [ i ]; jlist = firstneigh [ i ]; jnum = numneigh [ i ]; for ( jj = 0 ; jj < jnum ; jj ++ ) { j = jlist [ jj ]; f ac tor _ lj = spe c i a l _ lj [ s b m a sk ( j )]; f ac tor _ c oul = spe c i a l _ c oul [ s b m a sk ( j )]; j & = NEIGHMASK ; delx = xtmp - x [ j ][ 0 ]; MultiXscale Deliverable 2.6 Page 47 dely = ytmp - x [ j ][ 1 ]; delz = ztmp - x [ j ][ 2 ]; rsq = delx * delx + dely * dely + delz * delz ; jtype = type [ j ]; if ( rsq < c utsq [ itype ][ jtype ]) { r 2 inv = 1 . 0 / rsq ; // terms for re ac tion field dou b le rf _ f c tr _ 0 = 1 . 0 / pow ( sqrt ( rsq ), 3 . 0 ); dou b le rf _ f c tr _ 1 = ( epsilon _ rf [ itype ][ jtype ] - 1 . 0 ); dou b le rf _ f c tr _ 2 = 1 . 0 + 2 . 0 * epsilon _ rf [ itype ][ jtype ]; if ( rsq < c ut _ c oulsq [ itype ][ jtype ]) { for c e c oul = ( qqrd 2 e * qtmp * q [ j ]) * (( r 2 inv * sqrt ( r 2 inv )) - ( 1 . 0 / pow ( c ut _ c oul [ itype ][ jtype ], 3 . 0 ) * ( 2 . 0 * rf _ f c tr _ 1 / rf _ f c tr _ 2 ))); } else { for c e c oul = 0 . 0 ; } if ( rsq < c ut _ ljsq [ itype ][ jtype ]) { r 6 inv = r 2 inv * r 2 inv * r 2 inv ; for c elj = r 6 inv * ( lj 1 [ itype ][ jtype ] * r 6 inv - lj 2 [ itype ][ jtype ]); } else for c elj = 0 . 0 ; fp a irlj = f ac tor _ lj * for c elj * r 2 inv ; fp a ir c oul = f ac tor _ c oul * for c e c oul ; fp a ir = fp a irlj + fp a ir c oul ; f [ i ][ 0 ] += delx * fp a ir ; f [ i ][ 1 ] += dely * fp a ir ; f [ i ][ 2 ] += delz * fp a ir ; if ( newton _ p a ir || j < nlo ca l ) { f [ j ][ 0 ] - = delx * fp a ir ; f [ j ][ 1 ] - = dely * fp a ir ; f [ j ][ 2 ] - = delz * fp a ir ; } if ( f [ i ][ 0 ] > 10000 || f [ i ][ 1 ] > 10000 || f [ i ][ 2 ] > 10000 ) { std :: c out << " Pa ir LJC ut RF :: c ompute " << "\ n "; std :: c out << " f [ i ][ 0 ]: " << f [ i ][ 0 ] << "\ n "; std :: c out << " f [ i ][ 1 ]: " << f [ i ][ 1 ] << "\ n "; std :: c out << " f [ i ][ 2 ]: " << f [ i ][ 2 ] << "\ n "; } if ( efl a g ) { if ( rsq < c ut _ c oulsq [ itype ][ jtype ]) { e c oul = ( qqrd 2 e * qtmp * q [ j ]) * sqrt ( r 2 inv ) * ( 1 . 0 + ( rf _ f c tr _ 1 / rf _ f c tr _ 2 ) * ( pow ( sqrt ( rsq )/ c ut _ c oul [ itype ][ jtype ], 3 . 0 ))) - ( qqrd 2 e * qtmp * q [ j ]) * ( 1 . 0 / c ut _ c oul [ itype ][ jtype ]) * ( 3 . 0 * epsilon _ rf [ itype ][ jtype ] / rf _ f c tr _ 2 ); e c oul * = f ac tor _ c oul ; } else { e c oul = 0 . 0 ; } if ( rsq < c ut _ ljsq [ itype ][ jtype ]) { evdwl = r 6 inv * ( lj 3 [ itype ][ jtype ] * r 6 inv - lj 4 [ itype ][ jtype ]) - offset [ itype ][ jtype ]; evdwl * = f ac tor _ lj ; MultiXscale Deliverable 2.6 Page 48 Other classes can access potential energy by calling single () member function. Again, focusing only on the rf method, the energy is computed as Urf ( riαjβ ) = qiαqjβ 4 πε 0 riαjβ [ 1 + εrf − 1 2 εrf + 1 ( riαjβ rc ) 3 ] − qiαqjβ 4 πε 0 rc 3 εrf 2 εrf + 1 and scaled by factor f ac tor _ c oul . The single () method returns the interaction energy stored in variable eng , which represents a sum of the LJ and rf energy terms. } else evdwl = 0 . 0 ; } if ( evfl a g ) ev _ t a lly ( i , j , nlo ca l , newton _ p a ir , evdwl , e c oul , fp a ir , delx , dely , delz ); } } } if ( vfl a g _ fdotr ) viri a l _ fdotr _ c ompute (); } /* ---------------------------------------------------------------------- */ dou b le Pa ir LJC ut RF :: single ( int i , int j , int itype , int jtype , dou b le rsq , dou b le f ac tor _ c oul , dou b le f ac tor _ lj , dou b le & ffor c e ) { dou b le r 2 inv , r 6 inv , for c e c oul , for c elj , phi c oul , philj ; dou b le rf _ f c tr _ 1 = ( epsilon _ rf [ itype ][ jtype ] - 1 . 0 ); dou b le rf _ f c tr _ 2 = 1 . 0 + 2 . 0 * epsilon _ rf [ itype ][ jtype ]; r 2 inv = 1 . 0 / rsq ; if ( rsq < c ut _ c oulsq [ itype ][ jtype ]) { for c e c oul = ( for c e - > qqrd 2 e * a tom - > q [ i ] * a tom - > q [ j ]) * ( r 2 inv * sqrt ( r 2 inv ) - ( 1 . 0 / pow ( c ut _ c oul [ itype ][ jtype ], 3 . 0 ) * ( 2 . 0 * rf _ f c tr _ 1 / rf _ f c tr _ 2 ))); } else { for c e c oul = 0 . 0 ; } if ( rsq < c ut _ ljsq [ itype ][ jtype ]) { r 6 inv = r 2 inv * r 2 inv * r 2 inv ; for c elj = r 6 inv * ( lj 1 [ itype ][ jtype ] * r 6 inv - lj 2 [ itype ][ jtype ]); } else for c elj = 0 . 0 ; ffor c e = f ac tor _ c oul * for c e c oul + f ac tor _ lj * for c elj * r 2 inv ; dou b le eng = 0 . 0 ; if ( rsq < c ut _ c oulsq [ itype ][ jtype ]) { phi c oul = for c e - > qqrd 2 e * a tom - > q [ i ] * a tom - > q [ j ] * sqrt ( r 2 inv ) * ( 1 . 0 + ( rf _ f c tr _ 1 / rf _ f c tr _ 2 ) * ( pow ( sqrt ( rsq )/ c ut _ c oul [ itype ][ jtype ], 3 . 0 ))) - for c e - > qqrd 2 e * a tom - > q [ i ] * a tom - > q [ j ] * ( 1 . 0 / c ut _ c oul [ itype ][ jtype ]) * ( 3 . 0 * epsilon _ rf [ itype ][ jtype ] / rf _ f c tr _ 2 ); eng += f ac tor _ c oul * phi c oul ; } if ( rsq < c ut _ ljsq [ itype ][ jtype ]) { philj = r 6 inv * ( lj 3 [ itype ][ jtype ] * r 6 inv - lj 4 [ itype ][ jtype ]) - offset [ itype ][ jtype ]; eng += f ac tor _ lj * philj ; } MultiXscale Deliverable 2.6 Page 49 To compute the energy of a particle to be inserted during OBMD simulations, an additional method called single _ a tomisti c _ o b md () has been implemented. The key difference between single _ a tomisti c _ o b md () and the single () method lies in the arguments that these member functions accept. While the single () member function takes (local) indices of atoms i and j , the single _ a tomisti c _ o b md () function requires the charge of particle i ( qi ) to be inserted instead of index i . It is important to note that in this case, the particle " i " does not yet exist. return eng ; } /* ---------------------------------------------------------------------- */ dou b le Pa ir LJC ut RF :: single _ a tomisti c _ o b md ( dou b le qi , int j , int itype , int jtype , dou b le rsq , dou b le f ac tor _ c oul , dou b le f ac tor _ lj , dou b le & ffor c e ) { dou b le r 2 inv , r 6 inv , for c e c oul , for c elj , phi c oul , philj ; dou b le rf _ f c tr _ 1 = ( epsilon _ rf [ itype ][ jtype ] - 1 . 0 ); dou b le rf _ f c tr _ 2 = 1 . 0 + 2 . 0 * epsilon _ rf [ itype ][ jtype ]; r 2 inv = 1 . 0 / rsq ; if ( rsq < c ut _ c oulsq [ itype ][ jtype ]) { for c e c oul = ( for c e - > qqrd 2 e * qi * a tom - > q [ j ]) * ( r 2 inv * sqrt ( r 2 inv ) - ( 1 . 0 / pow ( c ut _ c oul [ itype ][ jtype ], 3 . 0 ) * ( 2 . 0 * rf _ f c tr _ 1 / rf _ f c tr _ 2 ))); } else { for c e c oul = 0 . 0 ; } if ( rsq < c ut _ ljsq [ itype ][ jtype ]) { r 6 inv = r 2 inv * r 2 inv * r 2 inv ; for c elj = r 6 inv * ( lj 1 [ itype ][ jtype ] * r 6 inv - lj 2 [ itype ][ jtype ]); } else for c elj = 0 . 0 ; ffor c e = f ac tor _ c oul * for c e c oul + f ac tor _ lj * for c elj * r 2 inv ; dou b le eng = 0 . 0 ; if ( rsq < c ut _ c oulsq [ itype ][ jtype ]) { phi c oul = for c e - > qqrd 2 e * qi * a tom - > q [ j ] * sqrt ( r 2 inv ) * ( 1 . 0 + ( rf _ f c tr _ 1 / rf _ f c tr _ 2 ) * ( pow ( sqrt ( rsq )/ c ut _ c oul [ itype ][ jtype ], 3 . 0 ))) - for c e - > qqrd 2 e * qi * a tom - > q [ j ] * ( 1 . 0 / c ut _ c oul [ itype ][ jtype ]) * ( 3 . 0 * epsilon _ rf [ itype ][ jtype ] / rf _ f c tr _ 2 ); eng += f ac tor _ c oul * phi c oul ; } if ( rsq < c ut _ ljsq [ itype ][ jtype ]) { philj = r 6 inv * ( lj 3 [ itype ][ jtype ] * r 6 inv - lj 4 [ itype ][ jtype ]) - offset [ itype ][ jtype ]; eng += f ac tor _ lj * philj ; } return eng ; } MultiXscale Deliverable 2.6 Page 50 In LAMMPS, a fix is a modular component that applies specific operations to a group of atoms during a simulation. Fixes are used to impose constraints, apply forces, modify atom properties, or perform other tasks at each timestep. The purpose of a fix in LAMMPS is to be flexible and extensible, allowing developers to add a custom functionality, and implement custom behavior by overriding relevant methods. The o b md extension is implemented as a LAMMPS fix, enabling grand-canonical molecular dynamics simulations with special boundary conditions and external forces. It operates by: Key Components of the obmd Fix A fix is defined in the input script using the fix command. For example: This applies the o b md fix to all atoms. MPI is used to distribute work across CPUs in a parallel simulation environment: Managing particle insertion and deletion in buffer regions, Applying external forces to maintain user-defined boundary conditions, Communicating particle data across processors using MPI. 1. Fix Class: The fix is implemented as a C++ class derived from the base F ix class. The base class provides a common interface and utility functions. 2. Initialization: Fixes are initialized during the setup phase of a simulation. This includes parsing input parameters, allocating memory, and preparing any data structures required for the fix. fix 2 a ll o b md 3. Interaction with LAMMPS MD Procedure: Fixes are executed at specific points during the simulation. This fix overrides two key methods of the F ix base class to perform its operations, as specified with the setm a sk () function. pre _ ex c h a nge () : Handles particle insertion and deletion. This function is invoked before atoms are exchanged between processors. It: Deletes particles that cross open boundaries. Computes the number of particles to be inserted based on buffer conditions. Inserts new particles into specified regions using the USHER or near algorithms. Deletes particles that may overlap or violate boundary conditions. post _ for c e () : Applies external forces to atoms in specified regions. This function is called after force calculations in the MD timestep. It: Computes momentum and shear forces based on user-defined parameters. Distributes forces across atoms using weighting functions (e.g., smooth or step distributions). Particle Deletion: The try _ deleting () function identifies particles to be deleted based on their positions. MPI is used to: Reduce deletion counts across processors. Ensure consistent deletion of particles across the distributed domain. Particle Insertion: The try _ inserting () function attempts to insert particles into specified regions. MPI is used to: Broadcast insertion attempts and results. Ensure particles are inserted without overlap across processors. Force Distribution: Functions like reg _ for c e () and reg _ for c e _ perp () compute and apply forces to atoms. MPI is used to: Sum forces across processors. A.3 OBMD Architectural Document Distribute forces proportionally to atom properties. MultiXscale Deliverable 2.6 Page 52 MultiXscale Deliverable 2.6 Page 53 Acknowledgements The authors are grateful to Dr. Staš Bevc (NIC) for useful discussion and assistance in optimizing the obmd extension. Funding from the European High Performance Computing Joint Undertaking (JU) under Grant agreement No. 101093169 and Ministry of Higher Education, Science and Innovation from Republic of Slovenia under Grant agreement No. C3360-24-452013 is also gratefully acknowledged. References Acronyms used GPU graphical processing unit OBMD open-boundary molecular dynamics MD molecular dynamics PBCs periodic boundary conditions ROI region of interest rf reaction field DPD dissipative particle dynamics SPC simple point charge LJ Lennard-Jones EOS equation of state RDF radial distribution function Software mentioned LAMMPS Large-scale Atomic/Molecular Massively Parallel Simulator OBMD open-boundary molecular dynamics URLs referenced Page ii https://www.multixscale.eu ... https://www.multixscale.eu https://www.multixscale.eu/deliverables ... https://www.multixscale.eu/deliverables Internal Project Management Link ... https://github.com/multixscale/planning/issues/149 [email protected] .. . mailto:[email protected] http://creativecommons.org/licenses/by/4.0 ... http://creativecommons.org/licenses/by/4.0 Page 7 OBMD GitHub repository ... https://github.com/chocolatesoup/OBMD-LAMMPS-extension/tree/main Page 13 dev.eessi.io ... https://www.eessi.io/docs/repositories/dev.eessi.io/ EESSI ... https://www.eessi.io/docs/ Page 14 OBMD GitHub repository ... https://github.com/chocolatesoup/OBMD-LAMMPS-extension/tree/main Citations [1] R. Delgado-Buscalioni, J. Sabli´c, and M. Praprotnik, “Open boundary molecular dynamics,” Eur. Phys. J. Spec. Top., vol. 224, pp. 2331–2349, 2015. [2] L. Delle Site and M. Praprotnik, “Molecular systems with open boundaries: Theory and simulation,” Phys. Rep., vol. 693, pp. 1–56, 2017. [3] E. G. Flekkøy, R. Delgado-Buscalioni, and P. V. Coveney, “Flux boundary conditions in particle simulations,” Phys. Rev. E, vol. 72, p. 026703, 2005. [4] R. Delgado-Buscalioni and P. V. Coveney, “USHER: An algorithm for particle insertion in dense fluids,” J. Chem. Phys., vol. 119, pp. 978–987, 2003. [5] M. Heidari, R. Cortes-Huerto, D. Donadio, and R. Potestio, “Accurate and general treatment of electrostatic interaction in hamiltonian adaptive resolution simulations,” Euro. Phys. J. Spec. Top., vol. 225, pp. 1505–1526, 2016.