# MATILDA.FT, a Mesoscale Simulation Package for Inhomogeneous Soft Matter

Zuzanna M. Jedlinska  Department of Physics and Astronomy, University of
Pennsylvania, Philadelphia, USA  Christian Tabedzki  Department of Chemical
and Biomolecular Engineering, University of Pennsylvania, Philadelphia, USA
Colin Gillespie  Department of Chemical and Biomolecular Engineering,
University of Pennsylvania, Philadelphia, USA  Nathaniel Hess  Department of
Chemical and Biomolecular Engineering, University of Pennsylvania,
Philadelphia, USA  Anita Yang  Department of Chemical and Biomolecular
Engineering, University of Pennsylvania, Philadelphia, USA  Robert A.
Riggleman<EMAIL_ADDRESS>Department of Chemical and Biomolecular
Engineering, University of Pennsylvania, Philadelphia, USA

###### Abstract

In this paper we announce the public release of a massively-parallel, GPU-
accelerated software, which is the first to combine both coarse-grained
molecular dynamics and field-theoretical simulations in one simulation
package. MATILDA.FT (Mesoscale, Accelerated, Theoretically-Informed, Langevin,
Dissipative particle dynamics, and Field Theory) was designed from the ground-
up to run on CUDA-enabled GPUs, with the Thrust library acceleration, enabling
it to harness the possibility of massive parallelism to efficiently simulate
systems on a mesoscopic scale. MATILDA.FT is a versatile software, enabling
the users to use either Langevin dynamics or Field Theory to model their
systems - all within the same software. It has been used to model a variety of
systems, from polymer solutions, and nanoparticle-polymer interfaces, to
coarse-grained peptide models, and liquid crystals. MATILDA.FT is written in
CUDA/C++ and is object oriented, making its source-code easy to understand and
extend. The software comes with dedicated post-processing and analysis tools,
as well as the detailed documentation and relevant examples. Below, we present
an overview of currently available features. We explain in detail the logic of
parallel algorithms and methods. We provide necessary theoretical background,
and present examples of recent research projects which utilized MATILDA.FT as
the simulation engine. We also demonstrate how the code can be easily
extended, and present the plan for the future development. The source code,
along with the documentation, additional tools and examples can be found on
GitHub MATILDA.FT repository.

GPU, Nvidia, CUDA, Coarse-grained simulations, soft matter, polymers, self-
assembly

††preprint: AIP/123-QED

##  I Introduction

Poymers are a ubiquitous type of material, important both in biological and
industrial settings. Polymers are an umbrella term, gathering macromolecules
composed of smaller, repeating monomers. In the industrial settings, polymers
are used extensively in the tire industry, in flexible composite material
production, and utilized as durable adhesives. In addition, they have been
exploited in more precise applications, such as drug delivery Indulkar  _et
al._ (2016) and the design of artificial catalysis centers Sureka  _et al._
(2019). This is possible due to the propensity of polymers to self-assemble
into higher-order structures, and their ability to undergo phase separation in
solution. Controlled phase separation has been exploited to create nano-
capsules with well-defined pore sizes, by first inducing phase separation in
the capsule shell and then flushing-out one of the components Shi and Weitz
(2017). Similar approaches using non-solvents to induce phase separation in a
polymer solution are common methods to produce polymer membranesHaase  _et
al._ (2017); Tree  _et al._ (2017). Design of these materials requires precise
knowledge of the thermodynamics and microstructure of polymer materials under
a variety of conditions.

Within a biological context, the most important natural polymers are nucleic
acids and proteins. The former are responsible for storing and propagating the
genetic information, whereas the latter act as enzymes to facilitate multiple
biochemical reactions within the cell, participate in active transport of
intracellular components, or serve as structural elements of the cytoskeleton.
A special class of proteins, called Intrinsically Disordered Proteins (IDPs)
have been shown to be a main constituent of the membraneless organelles Wei
_et al._ (2017). The functions of these intracellular compartments vary from
nucleic-acid biosynthesis and organization, to participation in stress, and
immune responses Banani  _et al._ (2017). Just like regular organelles, they
have a distinct chemical composition, and thus spatially separate specific
biological processes. Unlike regular organelles, however, they are not
enclosed by a lipid membrane, and thus membraneless organelles can be easily
formed and dissolved as the need arises. The deregulation of these process has
been implicated in cancer development and the onset of neurodegenerative
disease Kanaan  _et al._ (2020). Efficient prediction of the equilibrium
thermodynamics of these systems remains an ongoing challenge, although
significant recent progress has been madeChoi  _et al._ (2020).

Given their wide-ranging applications, inhomogeneous (macro- or micro-phase
separated) polymers remain a topic of intense investigation. Experimental
methods are available to study their properties and aid in the novel material
design, and the polymer science community has a rich history of close
collaboration between experiment, theory, and simulation. Computer-aided
approaches have had an ever increasing importance as a research tool as
computational power has become more ubiquitous. Due to this continuously
increasing computational speed and processing power, progressively bigger and
more complex systems can be simulated, allowing mesoscopic material
propertiers to be predicted before synthesizing physical samples.

Various open-source Molecular Dynamics (MD) codes have been released, which
are capable of simulating polymeric species on an atomistic or coarse-grained
level. Some notable examples include LAMMPS Thompson  _et al._ (2022), NAMD
Phillips (2022), and GROMACS Abraham  _et al._ (2015). LAMMPS can perform all-
atom simulations on polymer chains, using available force fields. It is also
equipped with biologically-oriented force fields, which enable coarse grained
simulations of biomolecules. In addition, the user can define their own
coarse-grained polymer model, and expand it to include the required
potentials. On the other hand, both GROMACS and NAMD have been specifically
designed to model biological molecules, such as proteins and nucleic acids on
a fully atomistic level. Polymer simulations are challenging in general due to
the wide range of length- and time-scales required for accurate simulation.

In many soft matter fields, particularly those involving the design materials
using polymers, polymer field theory and related techniques have played a
crucial role in the design of new materials and in the interpretation of
experimental resultsFredrickson (2006); Edwards (1965); Leibler (1980); Zheng
and Wang (1995); Milner  _et al._ (1988). Polymer field theories are developed
by beginning with a description of the system in terms of coarse-grained
potentials, such as chains obeying Gaussian statistics, Flory contact
repulsions governed by a $\chi$ parameter, partial charges on the various
species, etc. One first writes down the partition function for this particle
model, then using one of a variety of transformation techniquesFredrickson
(2006); Edwards and Freed (1970); Edwards (1965), decouples the particle
interactions and transforms the model to one where one molecule of each type
interacts with chemical potential fields generated by the various interaction
potentials. The field-theoretic approach is attractive because it enables a
variety of analytic analysis, such as the mean-field approximation which gives
rise to self-consistent field theory (SCFT) or a variety of loop expansions.
The particle-to-field transformation is formally exact, and there are examples
in the literature showing quantitative agreement between the particle and
field version of the modelKoski  _et al._ (2013). In more recent years,
several methods have developed to sample the original particle model
efficientlyDetcheverry  _et al._ (2008); Pike  _et al._ (2009); Daoulas and
Müller (2006); Ganesan and Pryamitsyn (2003); Chao  _et al._ (2017a) that
generally fall under the umbrella of theoretically-informed coarse-grained
models; in these methods, the underlying particle coordinates are retained,
and the particles are mapped to density fields to efficiently calculate the
non-bonded forces and energies. However, with a few notable
exceptionsSchneider and Müller (2019); Arora  _et al._ (2016); Cheong  _et
al._ (2020), there is a dearth of simulation packages available to perform
general simulations on either particle or field-based models, particularly
those that are designed specifically to run on Graphics Processing Units
(GPUs). A single code base that can simulate both particle- and fluctuating
field-theoretic simulations of identical (or very similar) molecular models
could allow readily switching between dynamic and equilibrium simulations as
well as assessment of the importance of fluctuations in any calculating
properties.

In this work, we present a first version of our code for Mesoscale,
Accelerated, Theoretically-Informed, Langevin, Dissipative particle dynamics,
and Field Theory, MATILDA.FT. MATILDA.FT is written from the ground-up
intended to be run on GPUs, and the bulk of the code is written using the CUDA
programming language. MATILDA.FT is capable of modeling both the systems
consisting of a few molecules, as well as those containing millions of
particles. Its strength lies specifically in being able to efficiently
simulate polymeric and other soft materials (e.g., liquid crystalline systems)
on a mesoscopic scale. On this scale, the coarse-grained interactions are
typically “soft” (finite at overlap) and the particle density high; in this
limit, it becomes more efficient to evaluate the non-bonded interactions using
density fields. These large-scale molecular assemblies of interest can
correspond to biomolecular coacervates in explicit solvent, artificially
synthesized ionomers, block copolymer melts, side-chain liquid crystalline
polymers, or polymer-infiltrated nanoparticle packings.

The outline of this paper is as follows. In Section I.1 we begin with a brief
history of the code development. Next, in Section II we outline the main
features of the code and available functionalities. In Section III, we
describe the structure of the models being used in the molecular dynamics and
field-theoretical simulations. Subsequently, in Section IV, we show how the
code is optimized for parallel execution on CUDA-enabled GPUs. In the Methods
Section V we present a more in-depth description of selected features
available in MATILDA.FT. In Section VI we outline how the code is organized,
with focus on its class structure, and extensibility options. Next, in Section
VII we show results for selected example systems. We end with the planned
developments in Section VIII, and a conclude with Section IX.

###  I.1 Brief History

The GPU-TILD code evolved out of a series of programs originally referred to
as dynamic mean-field theories (DMFT) developed by Chao, Koski, and
RigglemanChao  _et al._ (2017b) and augmented by several students in
subsequent years. Over time, it became clear that the approach was strongly
connected to existing methods that went by the moniker theoretically-informed
Monte Carlo (TIMC). The technical differences between the two is that our
approach tends to use higher-order particle-to-mesh schemes to ensure UV-
convergence and that time evolution is through Langevin dynamics in lieu of
Monte Carlo dynamics. As a result we have subsequently referred to the method
as theoretically-informed Langevin dynamics (TILD).Our group had several
versions of an internal code that was developed using openMP and later
openMPI. However, all of these codes suffered from a fairly rigid structure
that was tied to simulations of specific systems associated with specific
research projects and an unintuitive input file format.

During the first summer of the Coronavirus pandemic, RAR began developing a
basic TILD simulation package. The goal was to develop code intended to be
used on GPUs from the ground up, creating a more modular, object-oriented
structure. Once the utility of the code became clear, its use spread
throughout our research group, with substantial additions to the code base
being made by all of the co-authors of this work, most notably CT and ZMJ. The
addition of the field-theoretic simulations portion of the code began in late
2022 by RAR.

##  II Feature overview

In this section, we provide a brief overview of the features available in
MATILDA.FT. They are later described in more detail in the following sections,
IV and V.

MATILDA.FT can performed two types of simulations: TILD or field-theoretic
(FT) simulations. The TILD method is a hybrid particle/field approach where
explicit coordinates of the molecules are retained and used to calculate
bonded interactions while non-bonded forces (e.g. excluded volume and
electrostatics) are calculated by mapping the particles to a density field.
The highly coarse-grained nature of the interactions leads to a significant
speed-up compared to particle-only implementations of the same models, due to
the relatively high particle density in such models. The method is closely
related to the theoretically-informed coarse-grained modeling techniques
developed in the de Pablo groupDetcheverry  _et al._ (2008); Pike  _et al._
(2009) and other related techniquesDaoulas and Müller (2006); Ganesan and
Pryamitsyn (2003). In a FT simulation, on the other hand, particle coordinates
are integrated out completely and one must compute the statistics of the
molecular conformations in external fields generated by the other molecules.
FT simulations are especially powerful when the equilibrium properties of high
molecular weight polymers are of interest, specifically a large dimensionless
polymer concentration $C=n/(V/R_{g}^{3})$, where $n$ is the number of chains
in a system, $V$ is the volume, and $R_{g}$ is the radius of gyration of a
chain. The FT simulations and TILD methods are complimentary in this sense, as
the particle-based TILD approach is much more efficient at lower polymer
concentrations, and TILD can capture some aspects of polymer dynamics that FT
cannot.

### TILD Branch

Next, we move to the overview of the TILD module. Here, the simulations are
performed in the NVT ensemble, in a fully-periodic orthogonal box, either in
two or three dimensions. Although MATILDA.FT can perform simulations of free
particles, it has been designed specifically to efficiently model systems of
polymer melts and solutions. Polymers are modeled as discrete Gaussian chains
with monomers that are connected through harmonic springs, and the density of
each monomer is spread around its center through a convolution with a unit
Gaussian. The strength of the repulsive interactions between the species is
mediated through the Flory-Huggins $\chi$ parameter. Monomers can be either
neutral or charged. If they carry a net charge, then in addition to the
repulsive potential, they also interact through Coulombic electrostatic
forces. Regardless of their net charge, the monomers can be made polarizable
through the use of (classical) Drude Oscillators as detailed in Section V.1.

User-defined groups of particles are the basic structure to which operations
in MATILDA.FT are assigned. These operations can be either static or dynamic.
Static groups maintain the same set of members over the entire course of
simulation whereas dynamic groups periodically update their member lists,
based on a specific membership criterion. The basic set of operations acting
on particle groups are applying the repulsive and electrostatic force,
thermostating, and integration, in order to propagate coordinates in time.
Other operations that can be applied to particles include external/biasing
forces, spatial confinement, creation and breaking of dynamic bonds, or
performing on-the-fly property calculations. Some of these operations require
a distance-based neighbour list, MATILDA.FT provides different styles of
neighbour lists depending on the application.

The user interacts with the code through input scripts, written in a plain-
text format. Before the simulation is started, the entire script is read and
appropriate variables and data structures are initialized. The maximum number
of time steps and the time step size are parameters specified in the input
script. Then the system can undergo an optional equilibration period before
subsequently entering the production run stage. Different input formats are
required by the TILD and FT simulations. For the TILD branch, two files need
to be provided. The first one is the main input file, providing information
about simulation dimensionality, box size, density grid spacing, and
interaction potentials between particle types. It also defines the particle
groups and assigns integrators and forces to act on them during each time
step. The second file contains information about particle coordinates, types,
molecules they belong to, and optionally, their charge. Currently, this file
can be provided in the format consistent with the LAMMPS data file, in either
atomic or charge atom style. To allow the use of data generated by other
codes, the initial configuration can also be read from a GSD-format file,
developed by the Glotzer Lab Lab (2023). On the other hand, FT simulations
only require a single input file, as all the fields and molecular information
are initialized from within the input script.

A sample input script, providing basic TILD functionality is shown in Listing
1. Here, the simulation is performed in 2 dimensions, for 10,000 time steps.
The script specifies various log and output frequencies. The variable pmeorder
controls the order of interpolation for constructing the density fields. The
dimensionless time step is set to 0.005. The configuration data is read from
the input.data file. An external force midpush (along the y-axis, with the
dimensionless magnitude of 0.5) is applied to all particles. The integrator is
chosen to be GJF. The lines beginning with pair_style specify the repulsive
interactions between the selected particles types. A more detailed description
of all available options can be found in the documentation and in Section VI.

⬇

1Dim 2

2

3max\_steps 10001

4log_freq 1000

5binary_freq 1000

6pmeorder 1

7traj_freq 1000

8

9delt 0.005

10read_data input.data

11

12Nx 65

13Ny 65

14

15extraforce all midpush 1 0.5

16

17integrator all GJF

18

19pair_style gaussian 1 1 1.5625 1.0

20pair_style gaussian 2 2 1.5625 1.0

21pair_style gaussian 1 2 3.00 1.0

Listing 1: Example input script

### FT Branch

As the development of the FTS features of the code only began relatively
recently, the feature set is currently somewhat limited, but will expand
significantly in the coming months. Currently, FTSs are limited to mean-field
calculations as in self-consistent field theory (SCFT) with linear, discrete
Gaussian chain models. The molecules can be of arbitrary blockiness with an
arbitrary number of components, and the potentials implemented include the
Flory contact repulsion and the Helfand weak compressibilityHelfand (1975a).
More details about the interactions between the species are provided in
Section III. As detailed below, the key elements of the FTS implementation are
three classes: Potentials, which govern the non-bonded interactions and act on
Species. The Species class stores the total density of each chemical component
and is populated by individual Molecule classes. For example, an
A-homopolymer/B-homopolymer/AB-diblock copolymer blend would have two species
(A and B) and three molecules.

### Units

All simulations in MATILDA.FT are performed in reduced units. The base units
for TILD simulations are energy $[k_{b}T]$, mass $m$ of a monomer, and length
is normalized by the statistical segment size $b$. For the FT simulations, the
unit length scale is the radius of gyration computed for an ideal chain of a
specified reference chain length, $N_{r}$, $R_{g}=\sqrt{(N_{r}-1)/6}$. In
addition, the field variables are also scaled by this reference chain length
so that the natural interaction parameters become $\chi N_{r}$, $\kappa
N_{r}$, etc. By substituting an appropriate value, other derived units can be
obtained, such as time and force.

### Output and Additional Features

The user might want to retrieve the information generated during the
simulation. By default, a frame of the trajectory file is written with a
specified frequency. This trajectory file is consistent with the LAMMPS traj
file style. Since this operation slows down the program execution, an
alternative binary file can also be saved. Two binary outputs are written out
- one storing particle positions and another storing the density fields. These
can be later converted to a human-readable format with the post-processing
tools that are distributed with the code. The simulation can be initialized
from a configuration file, following LAMMPS atomic or charge input styles. It
can also be resumed from the previously generated trajectory file by using a
resume command in the input script. Each job run can also be split into an
equilibration and production sections by explicitly setting the number of
steps for each section in the input script. This writes the data to two
separate files, each of which can have their write rates modified
independently of other sections. This allows for data to be collected during
and post equilibration within a single run, without a need to manually resume
the simulation. The trajectory data can also be saved in the GSD format. More
details about the specific output options can be found in thedocumentation.

##  III Structure of theoretically-informed coarse-grained and field-theoretic
models

In this section we provide the necessary theoretical background to understand
the models used in MATILDA.FT, and the logic of the simulation workflow. The
starting ingredients for all of the modeling to be handled by MATILDA.FT are
highly coarse-grained models for soft-matter systems. For simplicity, we will
describe the basic structure in terms of a simple A-B Gaussian chain diblock
copolymer melt, though the generalization to other systems will become
apparent below. For a polymer melt with $n$ polymer chains each containing
$N_{A}+N_{B}=N$ monomers, the microscopic polymer densities are

$\hat{\rho}_{K}(\mathbf{r})=\sum_{j}^{n}\sum_{s}^{N_{K}}\delta(\mathbf{r}-\mathbf{r}_{j,s}),$
(1)

where $K$ is either species A or B, and $\mathbf{r}_{j,s}$ is the position of
the $s^{\mathrm{th}}$ bead on the $j^{\mathrm{th}}$ chain. The monomers on
each chain are typically connected via harmonic bonds,

$\beta
U_{0}=\sum_{j}^{n}\sum_{s}^{N-1}\frac{3}{2b^{2}}|\mathbf{r}_{j,s}-\mathbf{r}_{j,s+1}|^{2},$
(2)

where $b$ is the statistical segment size, and we have assumed equal $b$ for
species A and B. The Flory repulsion is written in one of two equivalent
formsKoski  _et al._ (2013); Villet and Fredrickson (2014); Weyman  _et al._
(2021) depending on whether the model is implemented as a field- or particle-
based model. In the particle-based approaches, we make the potential non-local
as

$\beta U_{1}=\frac{\chi}{\rho_{0}}\int d\mathbf{r}\int
d\mathbf{r}^{\prime}\,\hat{\rho}_{A}(\mathbf{r})\,u_{G}(|\mathbf{r}-\mathbf{r}^{\prime}|)\,\hat{\rho}_{B}(\mathbf{r}^{\prime}),$
(3)

where $u_{G}$ is a unit Gaussian potential,
$u_{G}(r)=(2\pi\sigma^{2})^{-\mathbb{D}/2}e^{-r^{2}/2\sigma^{2}}$,
$\mathbb{D}$ is the dimensionality of the system, and $\sigma$ controls the
range of the interactions. The standard Flory-Huggins model is recovered in
the limit $\sigma\rightarrow 0$. The final potential penalizes deviations of
the local density from the averageHelfand (1975b) $\rho_{0}=nN/V$,

$\beta U_{2}=\frac{\kappa}{2\rho_{0}}\int d\mathbf{r}\int
d\mathbf{r}^{\prime}\;[\hat{\rho}_{+}(\mathbf{r})-\rho_{0}]\,u_{G}(\mathbf{r}-\mathbf{r}^{\prime})\,[\hat{\rho}_{+}(\mathbf{r}^{\prime})-\rho_{0}]$
(4)

With these ingredients in hand, we can write the partition function as

$\mathcal{Z}=z_{0}\int d\mathbf{r}^{nN}\;e^{-\beta U},$ (5)

where $z_{0}$ is a prefactor that contains all of the self-energy terms,
factors accounting for molecular indistinguishability, and the thermal de
Broglie wavelengths. In equilibrium particle-based simulations, one is
primarily interested in calculating averages of quantities that can be
expressed as functions of the particle coordinates, $M(\mathbf{r}^{nN})$, as

$\langle M\rangle=\frac{1}{\mathcal{Z}}\int
d\mathbf{r}^{nN}\,M(\mathbf{r}^{nN})\,e^{-\beta U},$ (6)

and expressions for the usual thermodynamic quantities of interest, such as
the average density, energy, and pressure, can be readily obtained from
expressions commonly used in molecular dynamics simulationsAllen and Tildesley
(2017).

For field-theoretic approaches, we use a local potential but render the
densities non-local by distributing the point particles over a Gaussian
distribution $h(r)=(2\pi a^{2})^{-\mathbb{D}/2}e^{-r^{2}/2a^{2}}$, and the
Gaussian-distributed particle density is given by

$\breve{\rho}_{K}(\mathbf{r})=\int
d\mathbf{r}^{\prime}\;h(\mathbf{r}-\mathbf{r}^{\prime})\,\hat{\rho}_{K}(\mathbf{r}^{\prime})=[h\ast\hat{\rho}_{K}](\mathbf{r}),$
(7)

where the final equality introduces our short-hand notation for a convolution
integral. For the choice $\sigma^{2}=2a^{2}$, we can exactly re-writeKoski
_et al._ (2013); Villet and Fredrickson (2014); Weyman  _et al._ (2021) the
non-bonded potentials in Eqs. 3 and 4 as

$\beta U_{1}=\frac{\chi}{\rho_{0}}\int
d\mathbf{r}\;\breve{\rho}_{A}(\mathbf{r})\,\breve{\rho}_{B}(\mathbf{r}),$ (8)

and

$\beta U_{2}=\frac{\kappa}{2\rho_{0}}\int
d\mathbf{r}\;[\breve{\rho}_{+}(\mathbf{r})-\rho_{0}]^{2}.$ (9)

Using known Gaussian functional integrals Edwards (1965); Fredrickson (2006),
one can then exactly transform the particle-partition function in Eq. 5 to a
field-theoretic one of the form

$Z=z_{1}\int\mathcal{D}\\{w\\}\;e^{-\mathcal{H}[\\{w\\}]},$ (10)

where $z_{1}$ contains the constants from $z_{0}$ as well as the normalizing
factors from the Gaussian functional integrals,
$\\{w\\}=\\{w_{+},w_{AB}^{(+)},w_{AB}^{(-)}\\}$ is the set of chemical
potential fields, and $\mathcal{H}$ is the effective Hamiltonian governing the
weights of the microstates. For the diblock copolymer model considered here,
$\mathcal{H}$ takes the form

$\begin{array}[]{rcl}\mathcal{H}&=&\frac{C}{\chi N_{r}}\int
d\mathbf{r}\left([w_{AB}^{(+)}(\mathbf{r})]^{2}+[w_{AB}^{(-)}(\mathbf{r})]^{2}\right)\\\
&&+\frac{C}{2\kappa N_{r}}\int d\mathbf{r}[w_{+}(\mathbf{r})]^{2}-iC\int
d\mathbf{r}\;w_{+}(\mathbf{r})\\\ &&-n_{D}\log
Q_{D}[\mu_{A},\mu_{B}],\end{array}$ (11)

where the first line contains the potential fields that arise due to the Flory
interactionKoski  _et al._ (2013), the second line contains the terms that
arise from the Helfand potential, and the final line contains the excess
chemical potential of the polymers in a given field. The potential fields
$\mu_{A}$ and $\mu_{B}$ experienced by monomers A and B computed using

$\begin{array}[]{rcl}\mu_{A,c}(\mathbf{r})&=&\left\\{i(w_{+}+w_{AB}^{(+)})-w_{AB}^{(-)}\right\\}(\mathbf{r})/N_{r}\\\
&&\\\
\mu_{B,c}(\mathbf{r})&=&\left\\{i(w_{+}+w_{AB}^{(+)})+w_{AB}^{(-)}\right\\}(\mathbf{r})/N_{r},\end{array}$
(12)

with the smeared potential fields appearing in Eq. 11 calculated as
$\mu_{K}(\mathbf{r})=[h\ast\mu_{K,c}](\mathbf{r})$.

While the particle implementation can report qualitatively realistic dynamic
quantities, the FT implementation is strictly interested in equilibrium
quantities. Equilibrium averages are typically expressed as functionals of the
potential fields and calculated as

$\langle
M\rangle=\frac{1}{\mathcal{Z}}\int\mathcal{D}\\{w\\}\,M[\\{w\\}]\,e^{-\mathcal{H}}.$
(13)

Since $\mathcal{H}$ is typically complex-valued, sampling the integral over
the field configurations is non-trivial; this is typically accomplished
through the mean-field approximation, leading to SCFT, through complex
Langevin (CL) samplingFredrickson (2006); Parisi (1983); Klauder (1984), or
Monte Carlo samplingStasiak and Matsen (2013).

To update the chemical potential fields in either a CL or an SCFT calculation,
the effective “forces” on the fields must be obtained as functional
derivatives of $\mathcal{H}$. These can be obtained through explicit
differentiation,

$\mathbf{F}_{w}(\mathbf{r})=-\frac{\delta\mathcal{H}}{\delta w(\mathbf{r})},$
(14)

where $w(\mathbf{r})$ is one of the three fields
$w_{+}(\mathbf{r}),w_{AB}^{(+)},$ or $w_{AB}^{(-)}$. In a CL simulation, the
fields are sampled using an overdamped Langevin equation,

$\frac{\partial w}{\partial
t}=\lambda_{w}\mathbf{F}_{w}(\mathbf{r})+\eta(\mathbf{r},t),$ (15)

where $\eta(t)$ is a stochastic noise term chosen to satisfy the fluctuation-
dissipation theoremFredrickson (2006); Lennon  _et al._ (2008). An algorithm
that drives the system to an SCFT solution is easily obtained from Eq 15 by
simply setting the noise term to zero.

##  IV Parallel Algorithms, Data Structures and Object-oriented Programming
with CUDA and C++

The unique feature which sets MATILDA.FT apart from other popular MD codes,
such as LAMMPS or NAMD is that it has been designed from the beginning to
execute on CUDA-enabled GPUs. The code is intended for highly coarse-grained
models where the non-bonded forces can be evaluated using density fields and
not summing over neighboring pairs of particles. It uses a dedicated CUDA/C++
programming language in order to fully harness parallel capabilities. Its
model of parallelization differs from the conventional CPU domain
decomposition. Whereas on the CPU the groups of particles are assigned to
different processor based on their spatial arrangement, the GPU
parallelization occurs on the particle or individual grid location level,
where each thread is responsible for processing instruction for the selected
particle/grid point. It is simply handled by assigning a separate thread to
individual particles, by filtering the thread IDs. This is handled by the
first statement in each kernel call, provided in Listing 2.

⬇

1int list_ind = blockIdx.x * blockDim.x + threadIdx.x;

2if (list_ind >= ns)

3 return;

Listing 2: Kernel call to assign a thread to a particle id.

MATILDA.FT also makes extensive use of Thrust library, which is an extension
of C++ Standard Template Library (STL) to work with GPUs Bell and Hoberock
(2012). The Thrust Library provides dedicated storage containers (equivalent
to STL vectors in C++), which enable easier host-device communication and
avoiding the requirement for explicit cudaMemcpy calls. The Thrust Library
also makes available dedicated parallel algorithms to operate on these
containers and achieve better performance. In addition, by avoiding
complicated host-device memory transfer syntax, the use of thrust makes it
easy for those who do not have much GPU-programming experience to understand
and expand the MATILDA.FT source code.

###  IV.1 Neighbour lists

Neighbour lists keep track of the other particles present within a certain
distance of the center of the particle of interest. This information is
required, for example, when model reactive particles need to search for
reaction partners. In this way, each time the force is being applied, a
comparison of all possible inter-particle distances, an
$\mathcal{O}(N_{tot}^{2})$ operation, where $N_{tot}$ is the total number of
particles, is avoided. Instead, an $\mathcal{O}(N)$ scaling is achieved during
each function call, since only pre-computed neighbours are checked for
possible interactions. The process of creating the neighbour list is based on
cells, to which the simulation box is divided, and thus it also scales as
$\sim N$.

Each neighbour list requires two parameters: $r_{cutoff}$ and $r_{skin}$. They
correspond to the maximum range of the force using the neighbour list, and the
skin radius $(r_{skin}>r_{cutoff})$ respectively. The skin radius is the
actual radius of the sphere (or a disc in 2D) within which the search for the
possible neighbours occurs. By specifying the skin radius, some particles
which are initially beyond the range of the force are included in the
neighbour list, and thus the list can be reused over multiple time steps. They
neighbour list can be rebuilt with a specified frequency, or, alternatively,
an automatic trigger can be used. Currently, list rebuilding is triggered when
a particle moves a distance larger than half the skin radius.

Each cell is chosen to have a side length equal to the half of the specified
skin radius. The process of neighbour list building is divided into two
stages. In the first phase, particle ”binning” takes place, where each
particle in the group gets assigned to one cell in the box. Finding the
corresponding cell is trivially parallelized over all particles, scaling as
$\mathcal{O}(N)$, whereas their assignment the grid point requires an atomic
operation. Each cell stores the ID of its member particles, and keeps track of
their total count. Subsequently, the program loops over each particle in the
group and checks the distance to other particles present in the cells
overlapping with $r_{skin}$. Each cell is chosen to have a side length equal
to the half of the specified skin radius, resulting in 125 cells in 3D, and 25
in 2D. Only the particles within the specified radius are included as
particle’s neighbours. The full neighbor list is stored as a 1D thrust device
vector, so no host-GPU data transfer is required. Another 1D device vector
stores cell content and their respective occupation count. This method uses a
predefined maximum cell capacity. If the specified capacity is too small, it
gets automatically adjusted during the simulation run. The full neighbour list
stores all neighbours of each particle, effectively double counting each
particle pair. A ”half-list” stores, for each particle, only the neighbours
with index lower than its own index. Other, specialized lists are also
available, and interface with specific particle-particle operations. For
example, to aid in dynamic bonding, where each donor particle only stores the
information about the acceptors present within the skin radius. In the future
release, other, more efficient schemes of neighbour lists will be implemented.
For example, each particle can be responsible for the region covering only
half of its neighbourhood, such that when combined, the particles cover the
search space in a minimally overlapping manner.

###  IV.2 Bonded interactions

Bonded interactions are present whenever polymer chains are modeled.
Currently, they represent the harmonic springs connecting adjacent monomers of
the same molecule. The calculation of the resulting forces has been
parallelized to be performed on the GPU, and Newton’s third law is not used to
avoid the use of atomic operations. The contribution of the bonded
interactions is calculated individually for each particle by assigning it to a
separate tread. MATILDA.FT can also report the resulting bonded energy and the
contribution the bonded interactions make to the pressure virial coefficient.
These calculations are performed on the CPU.

Two common angle potentials are also implemented in MATILDA.FT. To enable
simulations of discrete worm-like chains, we have the cosine form

$u_{wlc}(\theta_{ijk})=\lambda\left[1+\cos(\theta_{ijk})\right],$ (16)

where $\lambda$ controls the stiffness of the potential and $\theta_{ijk}$ is
the inside angle between particles $i,j$ and $k$. The second potential
implements harmonic angles as

$u_{h}(\theta_{ijk})=k_{\theta}\left(\theta_{ijk}-\theta_{0}\right)^{2},$ (17)

with spring constant $k_{\theta}$ and equilibrium angle $\theta_{0}$. The
angle styles are specified in the input script along with the type (“wlc” or
“harmonic”), followed by the force constant (for both styles) and the
equilibrium angle if harmonic angles are used.

###  IV.3 Long-range interactions

Long-range interactions include repulsive interactions mediated by the Flory-
Huggins $\chi$ parameter, and the electrostatic forces acting between the
charged monomers. The distinctive feature of MATILDA.FT is the way in which it
handles these interactions. While bonded interactions use explicit coordinates
to calculate inter-particle distances, long-range interactions use the
mass/charge density field to compute the resulting forces. In this process a
Particle-to-Mesh (PM) method is used. In this scheme, the box is divided into
a discrete grid, with the number of grid-points in each direction being a
user-defined quantity. At the beginning of the simulation, a Fourier-space
representation of the inter-particle potential is calculated using Fast
Fourier Transform (FFT), and is stored for the rest of the simulation. Then,
at each time-step, every particle assigns its density contribution to nearby
grid points, using a spline interpolation scheme with the weights given in the
appendix of Ref. Deserno and Holm (1998). The order of the interpolating
spline can be chosen from 1 to 4, with higher order interpolation requiring
more computation. Regardless of the form of the pair potential $u(r)$, forces
are given in real space by

$f(\mathbf{r})=-\int d\mathbf{r}^{\prime}\nabla
u(\mathbf{r}-\mathbf{r}^{\prime})\rho(\mathbf{r}^{\prime}).$ (18)

The convolution in Eq. 18 is evaluated in Fourier space using FFTs, where it
becomes a simple multiplication. Then an inverse FFT is used to transform the
forces back into real space where the forces are interpolated back onto the
particle centers. The process is illustrated schematically in Figure 1. The
current pair potentials implemented in MATILDA.FT include Gaussian forms as
well as the nanoparticle-nanoparticle and nanoparticle-monomer forms
demonstrated in previous work by some of usKoski  _et al._ (2013).

Figure 1: Schematic illustration of the Particle-to-Mesh (PM) scheme.
Specifically, shown here is a first-order spline interpolation, where the
particle density is mapped to the two nearest grid points in each dimension.
The same spline weights are used to map the forces back to the particles.

###  IV.4 Liquid Crystalline Systems

We model liquid crystalline interactions through a modified Maier-Saupe (MS)
potential that is a discrete version of the McMillan modelMcMillan (1971). In
our implementation, the MS interactions involve two pairs of particles: one of
each pair is the “center” of the interaction $i$, and the other becomes a
partner particle $j$ that is used to define the local molecular orientation on
particle $i$,
$\mathbf{u}_{i}=\frac{\mathbf{r}_{i}-\mathbf{r}_{j}}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}$
(see schematic in Fig. 9 below). The local orientation vector is used to
define an orientation tensor for each particle
$\mathbf{S}_{i}=\mathbf{u}_{i}\mathbf{u}_{i}-\mathbf{I}/\mathbb{D}$, which is
mapped onto an orientation field similar to the density fields,

$\mathbf{S}(\mathbf{r})=\sum_{i}\mathbf{S}_{i}\,\delta(\mathbf{r}-\mathbf{r}_{i}).$
(19)

In Equation 19, the sum is over the LC centers. The orientation field
$\mathbf{S}(\mathbf{r})$ is then used to compute the Maier-Saupe potential
energyPryamitsyn and Ganesan (2004),

$\beta U_{MS}=-\frac{\mu}{\rho_{0}}\int d\mathbf{r}\int
d\mathbf{r}^{\prime}\;\mathbf{S}(\mathbf{r}):\mathbf{S}(\mathbf{r}^{\prime})\,u_{G}(|\mathbf{r}-\mathbf{r}^{\prime}|),$
(20)

where $\mu$ is the Maier-Saupe potential parameter and $u_{G}(r)$ is the
Gaussian potential that renders the interactions non-local. The forces are
derived by explicit differentiation and are presented in the documentation of
the code, and as we show below in Section VII, this model captures both
nematic and smectic A phases.

Particles that carry an orientation vector $\mathbf{u}_{i}$ are specified in
an additional input file that is similar in nature to the lists of bonded
partners. When specifying that the MS potential is to be used, the name of the
additional input file is also provided; contained in this file is a list of
pairs of particles $i$ and $j$ used to define the orientation vector
associated with particle $i$. This implementation allows for easy creation of
either main-chain liquid crystalline polymers or side-chain liquid crystalline
polymers, a detailed study of which will be the subject of a forthcoming
publication. Furthermore, by making one of the end sites within an LC mesogen
a different site type, one can indirectly control anchoring conditions at
phase boundaries by making making this other type more or less repulsive with
a particular species in the nearby phase.

##  V Methods

In this section, we describe in more detail how selected functionalities are
implemented in MATILDA.FT. When possible, we emphasize how the algorithms
presented below were optimized to take advantage of the GPU architecture to
accelerate performance. The summary of all available functionalities and
parameters can be found in the documentation.

###  V.1 Drude oscillators

Polarization can play an important role in phase behaviour of polymer
solutions, especially in biological context. Molecular polarizability has
influence on polymer solubility, and it can also modify how polymer chains
interact with salt ions present in the environment. In MATILDA.FT,
polarizability effects are introduced through the use of classical Drude
oscillators. In this approach, a “Drude particle” is attached to the parent
particle via a harmonic spring with stiffness $k_{D}$ and zero equilibrium
length. This Drude particle is assigned a partial charge $\delta q_{D}$ and
the partner particle $-\delta q_{D}$, such that the net charge of the two-
particle pair remains unchanged. The magnitude of the spring constant $k_{D}$
can be related to the molecular polarizability, with polarizability decreasing
with increasing stiffness. The Drude particle gets assigned a small mass, so
that it can be integrated with other particles using standard equations of
motion. This simplification circumvents the issue of treating polarizability
effects on the quantum-mechanical level, while still being able to reproduce
spatial variations in polarization. Drude particles do not participate in
excluded volume interactions, and thus the only forces acting on them are
electrostatic in nature. We note that our implementation is different from
those typically used in atomistic or more fine-scale coarse-grained models,
where the Drude particle is typically thermostatted independently at a low
temperature, enabling the polarizability to be estimated with the classical
expression $\alpha=\delta q_{D}^{2}/k_{D}$. Since our charges are distributed
over a unit Gaussian, this expression does not apply, and we have
parameterized our effective dielectric constant as a function of the various
parameters of the Drude oscillators ($q_{D},k_{D}$, and the spread of the
Gaussian, $\sigma$), which is presented below in Example Systems.

###  V.2 Dynamic bonding and Lewis acid-base pairs

In addition to static bonds, which are initialized at the beginning of the
simulation, and remain unchanged throughout its course, MATILDA.FT allows for
dynamic bonds to be formed between particles during the simulation run.
Dynamic bonds are created or destroyed based on the user-defined acceptance
criterion. Currently, Metropolis-Hastings acceptance criterion is used, with
an optional shift in the reaction energy. The energy of the bond is calculated
based on the extension of the harmonic spring assigned to the pair of bonded
atoms.

Dynamic bonds are created between a donor and an acceptor particle. Whether a
particle acts as a donor or an acceptor is specified in an external text file,
which maps these roles onto particle indices. To make the process
computationally efficient while executing on a GPU, MATILDA.FT utilizes a
dedicated neighbour list style. Only the donor particles are allowed to
initialize the process of bond making and breaking. They store information
only about the acceptor particles present in their vicinity. Two separate
lists keep track of bonded and free donor particles. During each
bonding/unbonding step, two kernels are launched in a randomly chosen order.
One kernel attempts to break some of the existing bonds, while the other
creates new ones. These kernels are only dispatched using the indices of
relevant donor particles - free donors for the bonding kernel, and bonded
donors for the bond-breaking one. Dynamic bonds can also be used for
simulating induced dipoles, for example the Lewis acid-base pairs. This
requires only a minimal change in the input script, namely, specifying the
magnitude of the charge to be assigned to each of the bonding partners.
Partial charge of opposite sign gets assigned to donors and acceptors upon
binding, so the net charge of the system remains constant.

###  V.3 Hydrodynamics using DPD Thermostat

Random noise is introduced into the simulation to account for the lost degrees
of freedom that are integrated out during the coarse-graining process.
Although Brownian Dynamics can be used to capture the behaviour of polymer
solutions, the resulting equations of motion do not satisfy the fluctuation-
dissipation theorem. Thus, they do not reproduce correct hydrodynamics. An
alternative method, available in MATILDA.FT, which obeys Navier–Stokes
equations, is the Dissipative Particle Dynamics (DPD) thermostating. A brief
overview of the method is provided below. For more detailed description, we
refer the reader to Groot  _et al._ (1999). In this method, three types of
forces act on particles. A conservative force $\mathbf{f}^{c}$, frictional
(dissipative) force $\mathbf{f}^{d}$, and a random force $\mathbf{f}^{r}$. The
conservative force corresponds to the gaussian repulsions, and to the
electrostatic interactions. The dissipative and random forces are pairwise
additive and thus the local momentum is conserved. The dissipative force on
particles i is given by Eq. 21,

$\mathbf{f}_{i}^{d}=-\frac{1}{2}\frac{\sigma^{2}}{k_{b}T}\sum_{j}(\omega(|\mathbf{r}_{ij}|))^{2}(\mathbf{v}_{ij}\cdot\mathbf{r}_{ij})\,\mathbf{r}_{ij},$
(21)

whereas the random force is given by

$\mathbf{f}^{r}_{i}=\sigma\sum_{j}(\omega(|\mathbf{r}_{ij}|))\gamma_{ij}\frac{1}{\sqrt{\Delta
t}}\hat{\mathbf{r}}_{ij},$ (22)

where $\hat{\mathbf{r}}_{ij}$ indicates a unit vector in the direction of
$\mathbf{r}_{ij}$.

As DPD involves the use of range-limited forces, it requires a neighbour list.
Since the interactions are pairwise-additive and symmetric, only a ”half” list
is needed, where each particle only stores the information about the particles
with the index lower than its own. This reduces the amount of required
computation, and ensures optimal GPU thread utilization. The DPD thermostat
should be used with the Velocity Verlet integrator. Code performance is
affected by the choice of the value of $\sigma$, and the time step. Strategies
for choosing the values and the reasoning behind them is outlined in the
documentation.

##  VI Code structure, extensibility, and flexibility

This section will give a brief overview of the organization of the code and
how its different parts cooperate with each other.

###  VI.1 Input Script

The first step of the simulation is selecting the method to be used, either
TILD or FT. This is passed as a command line argument when the program is
called. Using ./MATILDA.FT -particles will run the TILD simulation, whereas
-ft option will initialize an FT run. Additional command line arguments, like
the name of the input script to be used, are described in the documentation.
Two files are required for a TILD simulation. The main input file is
responsible for setting up simulation parameters, such as the dimensionality
of the system, size of the simulation box, grid density, time step size, and
the number of time steps to perform. The input script also defines the
interaction potentials between selected atom types (Potentials), and the
parameters used for calculating electrostatic forces - Bjerrum length and the
charge spreading length (currently uniform for the entire system, will be
allowed to vary between different types of particles in the future release).
The same file also specifies particle groups, along with the corresponding
neighbour lists, integrators, and additional forces to operate on selected
particles.

The second (data) file provides the initial positions of the particles, their
types, and the molecules they belong to. This data file also initialize the
static bonds and angles used in the simulation. Currently, the initial atom
configuration can be read either from the LAMMPS data file (in angle or charge
style) or from a GSD file.

###  VI.2 Code Organization

Figure 2: Schematic outline of code structure

The code takes advantage of C++ object oriented programming approach. It is
divided into classes, which interact with each other and exchange data as
needed. Each class is responsible for handling a particular functionality. The
base class serves as an interface used to interact with other parts of the
code. Then specialized sub-classes are derived from the base class, to provide
a specific functionality. This organization makes extending the code to
include new functionalities a relatively easy and straightforward process,
with simple integration of the new components into the existing code VI.4. For
example, the NeighbourList class is responsible for constructing and storing
the neighbour list for the selected group of particles. This neighbour list is
then used by the additional forces (created as a subclass of the ExtraForce
class) to accelerate the operations performed on this group. Depending on the
nature of the additional force, specialized neighbour lists can be used, in
order to further accelerate the performance. A diagram of class code
organization in shown in Fig. 2. Below, we provide a brief description of
selected classes. The outline of the code structure, along with the detailed
description of all class functionalities and options are provided in the
documentation.

####  VI.2.1 Global variables space

The global variable space holds the main data structures used in the
simulation. It stores the global arrays containing particle types, positions,
forces acting upon them, velocities, and static bonds, arranged according to
the particle ID. Reading of the input script is also handled at this level.
Before the beginning of the simulation, these arrays are initialized and then
periodically updated by other classes during the time-stepping process. In the
future release these structures will be placed in a separate Box class, to
closely resemble to organization of the FT branch of the code.

####  VI.2.2 Group Class

The Group base class provides data structures which store indices of the
member particles. It sets device-specific variables (BLOCK and GRID sizes)
that are used in kernel calls dispatch on the group particles. All forces and
neighbour lists in MATILDA.FT operate on specific groups. Pointers to each
group object are stored in a globally accessible vector array. Each group is
assigned a unique name, which is used to pass its pointer to their classes.
Groups can be static or dynamic. Static groups are initialized at the
beginning of the simulation and their content remain unchanged over the
simulation course. Dynamic groups, on the other hand, periodically check and
update their members, based on the specified membership criterion. A special
group, named ”all”, is initialized by default at the beginning of the
simulation, and contains all particles in the simulation box. Currently, two
static group types are available - grouping by particle type or by its global
id. Type-based groups collect all the particles with the same type (as
specified in the input.data). Id-based groups require the user to provide an
external plain text file which has contains the indices of the particles to be
included in the group. Currently available dynamic group style, “regions”,
allows the user to define a separate region in space (along all or only
specific axis). Particles found within that region get assigned to the group.
The can be easily extended to incorporate user-defined groups, see Section
VI.4 below.

####  VI.2.3 Neighbour List Class

Like all operations in MATILDA.FT, neighbour lists act on groups. For the
purpose of neighbour list building, the simulation box is assigned a grid
which divides it into discrete, non-overlapping cells. Interfacing with a
specific group of particles and the division of the simulation box into non-
overlapping cells is handled by the NList base class.

Building up on this basic functionality, more specialized sub-classes of
neighbour lists are created. The simplest one, distance, is the ”full”
neighbour list, which stores the entire information about particle proximity,
and thus double-counts each pair of neighbours. A slightly more elaborate sub-
class, half_distance, stores only half of this information. More specifically,
each particle only keeps track of the particles which have a lower id than
their own. This method also avoids unnecessary calculations when pairwise
interactions are present. Otherwise, a branching statement needs to be
included, to perform the computation on only one member of the particle pair.
Alternatively, the same calculation is performed for both members.
Unfortunately, thread divergence is not permitted in the GPU model of
execution, and thus a branching statement forces all the threads to wait until
the other ones have competed the first branch. This waste of resources is
avoided by using an appropriate neighbour list.

A dedicated neighbour list has been designed specifically to support the
dynamic bonding functionality. Briefly, only the donor particles, which
initialize the bonding, are assigned the neighbours. These neighbours are then
filtered to only include the acceptor particles, as only the donor-acceptor
combination can create a valid bonded pair. By pre-calculating the possible
pairs at this step, these checks can later be avoided in the kernel calls
(preventing possible thread divergence).

####  VI.2.4 ExtraForce Class

In addition to electrostatic and repulsive interactions, selected groups of
particles can be subject to additional user-defined forces. These are
specified in the input script using the extraforce command. The ExtraForce
base-class is responsible for assigning the force to the specific group of
particles, and ensuring it is applied to this group at the specified time-
steps (wither each step or user-defined frequency). In addition, some range-
limited forces require a neighbour list to restrict the search space only to
the particles present within the specific range. Currently available forces
are:

  * •

Wall \- which enables the particles to be confined within a specific region or
to simulate surface interactions. The user can chose from available wall-
particle potentials, or specify their own form of interaction, by extending
the source code.

  * •

Langevin \- Adds random noise to the selected group of particles. Can be used
with Velocity Verlet integrator to simulate Brownian dynamics.

  * •

Midpush \- Adds a force to push the selected group of particles towards the
center of the box along a specified axis.

  * •

DPD \- Dissipative Particle Dynamics. This subclass of ExtraForce provides an
alternative way to introduce random noise into the simulation, and should be
used along with the Velocity Verlet integrator. In contrast the Langevin
thermostat, however, it is pair-wise additive and conserves local momentum.
Thus it capable of correctly reproducing hydrodynamic behaviour of the system.
Since the force acts over a limited range, a neighbour list needs to be
constructed for the particles of interest.

  * •

Lewis \- This additional forces can be used to introduce dynamic bonds in the
simulation. Whereas static bonds are initialized at the beginning of the
simulation, and remain unchanged, dynamic bonds can be formed and broken
according to the specified acceptance criterion. This force requires a
specialized neighbour list (bonding) which has been designed to optimize the
required computations.

More details about the ExtraForce class can be found in the documentation.

####  VI.2.5 Compute Class

The Compute class is responsible for performing on-the fly calculations of
properties of the system. This enables the user to monitor the evolution of
the system in real time and also saves time spent on post-processing.

  * •

Average Structure Factor $\langle S(k)\rangle$ \- this compute provides the
information about the average static structure factor of the particle system.
The static structure factor, given by

$S(k)=\frac{1}{N}\langle\hat{\rho}_{k}\hat{\rho}_{-k}\rangle,$ (23)

is defined as the correlation function of the system density represented in
the Fourier space. The density is given by
$\rho(r)=\sum_{i=1}^{N}\delta(r-r_{i})$, and in Fourier space it becomes
$\hat{\rho}_{k}=\sum_{i=0}^{N}e^{ik\cdot r_{i}}$ Zhang (2016). It performs the
calculation and writes the data to an external file according to user-
specified frequency.

  * •

Chemical Potential $\mu$ \- this compute can be used to calculate chemical
potential of the given species in the system using a chain deletion method.
User provides a range (based on global molecule ID) of molecules to operate
on. These particles are removed at random from the system, and the change in
energy upon removal is calculated. This data is written to an external file.

####  VI.2.6 Integrator Class

In order to solve the equations of motion and propagate the particle
coordinates in time, numerical integration is required. In MATILDA.FT, three
different numerical algorithms are available, and are briefly described below:

  * •

Velocity-Verlet (VV). Needs to be coupled with additional thermostat.
Available thermostats include Langevin noise or Dissipative Particle Dynamics,
which of which are part of the ExtraForce class.

  * •

Euler-Maruyama (EM). Generates the thermal noise internally during the update
and serves as the simplest stochastic integration scheme to implement.

  * •

Grønbech-Jensen and Farago (GJF). Generates thermal noise internally during
the update, and We find that this algorithm allows time steps up to 10x larger
than the EM algorithm with no loss of accuracy or stability.

###  VI.3 FTS Implementation

####  VI.3.1 Class structure

Figure 3: Basic actions and roles of the FTS Classes in MATILDA.FT.

As the FTS branch was begun more recently, many planned features are still in
development. There is also a difference in class organization between the
older (TILD), and the newer (FT) branches. The classes that comprise an FT
simulation are more tightly integrated and the scope of object-oriented
organization is larger, as compared to the TILD branch which still uses global
variables.

A field-theoretic simulation lives in an FTS$\\_$Box class, which contains
three key classes: Potentials, Molecules, and Species (see Figure 3 for
graphical outline of FT branch organization). The Potentials class performs
all of the functions that are related to the the various non-bonded
interactions, including updating the potential fields associated with a
particular interaction. The densities that show up in the effective forces are
taken from the Species class, which serves as a container for these densities.
Species generates the unsmeared chemical potential fields,
$\mu_{K,c}(\mathbf{r})$, by looping over the interaction potentials and
accumulating the relevant potential fields. Next, the Molecules class takes
these potential fields, applies any density smearing that may be necessary,
and computes both the center and smeared density fields. The smeared density
fields are then accumulated into the relevant Species class. The general flow
of the code is summarized in Figure 4.

Figure 4: Outline of an FTS simulation as implemented in MATILDA.FT. The
termination condition could be convergence to within a prescribed tolerance in
SCFT or reaching the maximum desired number of time steps in a CL simulation.

####  VI.3.2 Field update schemes

Currently, two schemes have been implemented to update potential fields, in
order to evolve in time equations of motion such as Eq. 15. The
straightforward explicit Euler-Maruyama (EM) integration scheme discretizes
the equation in time and uses the forces at the current time to estimate the
field configurations at a future time as

$w^{t+\delta t}(\mathbf{r})=w^{t}+\delta
t\,\lambda_{w}F_{w}^{t}(\mathbf{r})+\sqrt{2\delta
t\lambda_{w}}\zeta_{t}(\mathbf{r}),$ (24)

where $\delta t$ is the size of the time step and $\zeta_{t}(\mathbf{r})$ is
purely real Gaussian noise with unit variance that is uncorrelated in both
space and time. As mentioned above, simply neglecting the noise term converts
this algorithm to one that drives the system to a mean-field solution.

The other algorithm that is implemented is a 1st-order, semi-implicit updating
scheme, that has been shown to allow for time steps significantly larger than
allowed by the EM schemeLennon  _et al._ (2008); Ceniceros and Fredrickson
(2004); Villet and Fredrickson (2014). In this approach, one derives an
approximate expression for the force $F_{w}^{t,lin}(\mathbf{r})$ that is
linear in the potential field. In real-space, these expressions take the form
of a convolution

$F_{w}^{t,lin}(\mathbf{r})=\int
d\mathbf{r}^{\prime}\;A_{w}(\mathbf{r}-\mathbf{r}^{\prime})\,w(\mathbf{r}^{\prime}),$
(25)

where $A_{w}(\mathbf{r})$ is the linear coefficient. As a result of this
convolution, the 1S updating scheme is most effectively handled in Fourier
space where we have

$F_{w}^{t,lin}(\mathbf{k})=A_{w}(\mathbf{k})\,w(\mathbf{k}).$ (26)

To affect the semi-implicit scheme, Eq. 24 is written in Fourier space and
modified by subtracting the linear term at $t+\delta t$ and adding it at $t$
giving,

$\begin{array}[]{rcl}w^{t+\delta t}(\mathbf{k})&=&\delta
t\,\lambda_{w}\left[F_{w}^{t}(\mathbf{k})+A_{w}(\mathbf{k})\,w^{t}(\mathbf{k})-A_{w}(\mathbf{k})w^{t+\delta
t}(\mathbf{k})\right]\\\ &&\\\ &&+w^{t}(\mathbf{k})+\sqrt{2\delta
t\lambda_{w}}\zeta_{t}(\mathbf{k}).\end{array}$ (27)

We note that $\zeta_{t}(\mathbf{k})$ is generated as a spatially uncorrelated
noise field in real-space that is explicitly Fourier transformed. Equation 27
can be readily solved for the field at $t+\delta t$ giving

$w^{t+\delta t}=\frac{w^{t}+\delta
t\,\lambda_{w}\left[F_{w}^{t}+A_{w}\,w^{t}\right]+\sqrt{2\delta
t\lambda_{w}}\zeta_{t}}{1+\delta t\lambda_{w}\,A_{w}},$ (28)

where we have suppressed the wavevector dependence for brevity.

The functional form of the linear coefficients $A_{w}$ generally contain one
or two contributions that have a stabilizing effect on the time
integrationLennon  _et al._ (2008); Ceniceros and Fredrickson (2004). The
first arises from the terms that are quadratic in the fields in $\mathcal{H}$
(e.g., the first two lines in Eq. 11); this term is included for every type of
interaction potential. During the initialization of an FT simulation, the
Potentials class adds this relevant term to $A_{w}$. The second contributions
are the linear approximates of the density operators, which involve
convolutions of Debye functions with the potential fields; these contributions
are handled by the Molecules class during initialization.

####  VI.3.3 Molecule Types

Currently, the only implemented molecule type is a linear, discrete Gaussian
chain with an arbitrary number of blocks. This class handles the calculation
of the chain propagators, and during initialization the code automatically
checks whether the molecule is symmetric to avoid calculating the
complimentary propagator if possible. The other key step taken in the
initialization is to accumulate the relevant Debye-function-like contributions
to the linear coefficients associated with each potential. From previous
studies Lennon  _et al._ (2008); Ceniceros and Fredrickson (2004), not all
terms that show up in the precise linear expansion of the force are
stabilizing; to that end, we do not include the Debye terms in the
$w_{AB}^{(-)}$ field, but they are included in the $w_{+}$ and $w_{AB}^{(+)}$
fields.

###  VI.4 Extensibility

MATILDA.FT has been designed to be easily extensible by other users, according
to their specific needs. The process of expanding the source code is
simplified thanks to the division into classes, consistent with C++ philosophy
of Object Oriented Programming. Due to this organization, the user needs only
to understand how the base class interfaces with other parts of the program,
and does not need to rewrite the entire logic. Simple functionalities can
easily be added by inheriting the capabilities of provided base classes or
modifying the existing once.

⬇

1/* kernel function to update group members

2 based on their position and type */

3__global__ void d_CheckGroupMembers(

4 const float* x, //position array

5 thrust::device_ptr<float> d_wall_data,

6 const int n_walls, // number of walls

7 thrust::device_ptr<int> d_all_id,

8 const int ns, // group size

9 const int Dim, // Dimensionality

10 const int* tp){

11 // tp[] array stores particle types

12 int list_ind = blockIdx.x * blockDim.x + threadIdx.x;

13 if (list_ind >= ns)

14 return;

15 int ind = list_ind;

16 for(int i = 0; i < n_walls; ++i){

17 int j = int(d_wall_data[3 * i]);

18 float low = d_wall_data[3 * i + 1];

19 float high = d_wall_data[3 * i + 2];

20 float xp = x[ind * Dim + j];

21 d_all_id[ind] = ind;

22 if (xp>=low && xp<=high && tp[ind]==1)

23 // additional type check

24 d_all_id[ind] = ind;

25 else

26 d_all_id[ind] = -1;

27 } // i < n_walls

Listing 3: Code extensibility example

This simplicity is demonstrated in Listing 3, where a new group is produced.
It is based on the dynamic group which only includes particles present in the
specific region. Here, in addition to spatial criterion, only particles of a
particular type (here 1) are considered group members. Here, only the modified
GPU kernel is shown for brevity. Two changes to the original code were
introduced to achieve new functionality. The first one on line 10 (an
additional parameter in the function call), and another on line 22 (type
check). This file can later be incorporated into a new subclass of the Group,
compiled with the rest of the source code, and will work seamlessly in the
simulation. Included with the source-code is also a pre-made make file which
makes compilation of additional components easy, only requiring that the be
added to the source list. The full example, along with a step-by-step
explanation, and make file description is available in the examples/extend
folder in the GitHub repository.

##  VII Example systems

###  VII.1 Coacervate

In this example, a small system consisting of the total of 434 molecules, each
with a degree of polymerization, $N=82$, has been simulated. Half of these
molecules carry positively charged monomers, with the other half having each
monomer with negative charges of the same magnitude. No explicit solvent is
present. The system starts in a random, homogeneous phase and over the course
of simulation phase separates into polymer-rich and polymer-depleted regions
as coacervation occurs. The snapshots from the beginning (left), and the end
(right) are shown in Figure 5. On the Nvidia Quadro RTX 5000 GPU, the
simulation took 1399 seconds to perform 2,000,000 time steps, with the
resulting speed of 1429.6 ts/sec. This example can be found in the GitHub
repository, in the examples/Coacervate directory. The full movie of the time
evolution of the system is available in the supporting information.

Figure 5: Initial (left) and final (right) snapshot from the simulation of
coacervating binary mixture. This simulation is a particle-based
implementation of the model considered previously as a field-theory by one of
usRiggleman  _et al._ (2012) with a dimensionless excluded volume parameter
$B=0.05$ and dimensionless Bjerrum length $E=10000$.

###  VII.2 Hydrodynamics using DPD integrators

To demonstrate the effect of using DPD thermostat, we simulated a simple
binary mixture of uncharged particles in a two-dimensional box. The value of
$\kappa$ in the Helfand potential has been chosen to be equal 1 and
$\chi_{AB}=5$ to induce strong phase separation. The value of $\sigma$ was
arbitrarily set to 0.5, with $r_{cutoff}=1.0$. The neighbour list update
frequency was set to every 6 time-steps, and $r_{skin}=2.5$. On the Nvidia
Quadro RTX 5000 GPU, the simulation took 5259 seconds to perform 1,200,000
time steps, with the resulting speed of 228.18 ts/sec. In Figure 6, the insets
show representative snapshots from the course of the simulation along with the
corresponding structure factor for the red species in the main figure.

Figure 6:  Static structure factor $S(k)$ calculated on the red species
undergoing spinodal decomposition at various time steps (TS) using the DPD
integrator. The insets show representative snapshots of the system at the same
times.

###  VII.3 Polarizable diblock co-polymer with explicit solvent

In this example, 1330 diblock co-polymer chains with $N=74$, and blocks of
equal size, were simulated in explicit solvent. Both the solvent and polymer
chains are were made polarizable through the use of classical Drude
oscillators. Since in this system only the polymer concentration of the
condensate was of interest, all simulations were started with the polymers in
a dense “slab” configuration. In this configuration all particles are biased
to migrate towards the middle of the simulation box, creating a homogeneous,
dense polymer phase. During the production run this bias is removed, and the
slab is allowed to expand.

Figure 7: Plot of the density obtained in the dilute and dense phase as the
value of $\chi$ between the monomers and solvent is varied. Included are also
renders of the three-dimensional structure of the system corresponding to the
selected data points. Positively charged monomers are displayed in red, while
negatively charged ones are colored blue.

The parameters for the Drude oscillator attached to the solvent molecules have
been chosen as to reproduce the dielectric constant of water. In this way, the
Bjerrum length can be set to the value it has in the vacuum, and dielectric
screening is then emergent from the polarizable solvent. The calculated
dielectric constant for chosen combination of parameters are shown in Figure
8. The dielectric constant of water is around 78.4, so that optimal choice of
parameter corresponds to $a_{0}\approx 0.5$ and $k_{D}\approx 1.0$ In Figure 7
we show the plot of the reduced concentration $c^{*}$ of the dense and dilute
phases phase, and the corresponding value of $\chi$ between the polymer
monomers and the solvent molecules. We also include corresponding snapshots of
final structure of the system.

Figure 8: Calculated values of dielectric constant for the solvent molecule,
as a function of charge spreading length, $a_{0}$ and the spring constant of
the Drude oscillator, $k_{D}$

###  VII.4 Liquid crystals

A simple model of a pure liquid crystal was simulated where the mesogen was
discretized into three interaction sites with the Maier-Saupe (MS) interaction
taken from the center of the mesogen, see the inset in Figure 9. Bonds between
adjacent liquid crystal sites used a force constant $k_{b}=100$ and
equilibrium distance $1$, and the orientation was maintained with a worm-like
chain angle potential with prefactor $\lambda=50$. Finally, an additional
Helfand potential was employed to maintain an approximately uniform density
with $\kappa=11$ and the total site density was $\rho_{0}=3$. This combination
of the anisotropic molecular shape with a MS interaction taken from its center
makes the model similar to the McMillan mean-field model McMillan (1971).

Figure 9: Liquid crystalline order parameter as a function of the strength of
the Maier-Saupe parameter $\mu$ calculated during a simulation where $\mu$ was
continuously ramped throughout.

Figure 9 shows the average liquid-crystalline order parameter $\lambda$ that
is calculated on the fly as the MS $\mu$ parameter is increased from 0 to 120.
$\lambda$ is taken as 3/2 times the largest eigenvalue of the average
$\mathbf{S}$ tensor. When $\mu\approx 45$, a sharp increase in
$\langle\lambda\rangle$ indicates the isotropic-to-nematic transition.

##  VIII Planned future developments

In this section we present the near-future and long term goals regarding the
development of MATILDA.FT. For the TILD branch, the main modification that is
currently in progress, will be converting to a fully object-oriented style, to
closely follow the design of the newer, FT branch. This will make the two
parts of code operate more seamlessly, and will facilitate future
modifications and make extensibility easier. Afterwards, we will focus on
optimizing existing algorithms, specifically modifying the calculation of
repulsive interactions based on Gaussian density.

The long-term goals for the TILD branch include introducing enhanced sampling
techniques, such as Umbrella Sampling and Gibbs Ensemble Monte Carlo. This
would enable the study of systems that suffer from being stuck for extended
periods of time in metastable potential minima, such as polypeptides. We are
also planning on expanding the software to be able to perform rigid body
dynamics, which can be used to model patchy particles with specific
interaction. Finally, we plan to harness the power of using multi-GPU
architectures, to aid in simulations with multiple boxes. This feature would
greatly facilitate simulations using enhanced sampling techniques like replica
exchange, parallel tempering or Gibbs Ensemble MC. A “quality of life” feature
we plan to implement is internal routines to generate starting configurations
so that one could simply specify the polymer architecture, length, and
density, and the code provides a starting configuration. Finally, the existing
neighbor lists could be used to implement slip-springs to capture entanglement
effectsRamírez-Hernández  _et al._ (2015, 2017), which would enable studies
of, for example, viscoelastic phase separation. As the user community grows,
we will try to implement features that are in high demand, and would provide
useful tools to the scientific community.

We plan to implement several key features for the FT simulation methods that
should be available in the near term. While the time evolution equations of
the fields described above are presented from the perspective of a complex
Langevin (CL) simulation, the CL equations are not yet implemented, and all of
the equations of motion are currently noise free. In addition, we plan to
implement monomer smearing with unit Gaussians to regularize the models
against ultra-violet divergencesKoski  _et al._ (2013); Riggleman  _et al._
(2012); Villet and Fredrickson (2014). Finally, the inclusion of electrostatic
interactions, including polar and polarizable polymer monomersMartin  _et al._
(2016), is planned in the near future. Further in the future we plan to
implement more exotic polymer arcitectures such as bottlebrush and star
polymers. Finally, a variety of nanoparticles, including anisotropicKoski  _et
al._ (2013) and polymer-grafted particlesChao  _et al._ (2014); Koski  _et
al._ (2015) as both field-based particles and explicit “hybrid” particlesSides
_et al._ (2006) are planned for the future.

##  IX Concluding Remarks

In this paper we presented MATILDA.FT, an open source mesoscale simulation
software, which utilizes GPU architecture and primarily a CUDA/C++ programming
language to achieve massive parallelism over thousands or millions of
particles. The particle-based simulations are primarily designed for highly
coarse-grained potentials that are finite on contact and where the particle
density will be relatively high so there is a gain in overall efficiency by
evaluating the non-bonded interactions using density fields. It has already
proven to be able so simulate a vast array of distinct systems and can be
easily extended to incorporate new ones, thanks to the object-oriented
implementation. As far as we are aware, MATILDA.FT is the first published
open-source software to combine both coarse-grained Langevin dynamics and
field-theoretic simulation frameworks into a single code base. The code has
been written with new users in mind, and its use and basic extensibility do
not require specialized programming knowledge. This will make it a valuable to
the broad part of the scientific community, providing a powerful computational
tool to both theorists and experimentalists across many fields of material
science and bio-engineering. We plan on continuously extending and improving
the existing code, and building a user community where scientists can share
their ideas and experience. As with many other open-source codes, such as
LAMMPS, this approach is extremely beneficial and allows for a faster, and
more impactful software development, to keep up with upcoming scientific
challenges.

##  X Acknowledgements

This work used Bridges-2 GPU at Pittsburgh Supercomputing Center through
allocation DMR150034 from the Advanced Cyberinfrastructure Coordination
Ecosystem: Services & Support (ACCESS) program, which is supported by National
Science Foundation grants #2138259, #2138286, #2138307, #2137603, and
#2138296. This work was supported by the National Science Foundation through
grant MRSEC/DMR-1720530 (R.A.R. and partial Z.M.J.), NSF awards OISE-1545884
(C.T.) and CHE-220375 (C.G.), and the Department of Physics and Astronomy at
the University of Pennsylvania (Z.M.J.). The images used in this publication
were generated using the free visualization software Ovito version 2.9.0
Stukowski (2010).

##  XI Data Availability

The source code is open source and available under the Gnu public license
(GPL2) at www.github.com/rar-ensemble/MATILDA.FT. The “examples” folder
contains all of the input files needed to reproduce the simulations presented
herein.

##  XII Author Contributions

Z.M.J: writing/original draft (lead); software (supporting); Validation
(equal); Visualization (lead); C.T: writing/review and editing (supporting);
software (supporting); C.G: writing/review and editing (supporting); software
(supporting); N.H: software (supporting); A.Y: software (supporting); R.A.R:
writing/original draft (supporting); review and editing (lead); software
(lead); conceptualization (lead); Validation (equal); Visualization
(supporting).

## References

  * Indulkar  _et al._ (2016) A. S. Indulkar, Y. Gao, S. A. Raina, G. G. Z. Zhang, and L. S. Taylor, Molecular Pharmaceutics 13, 2059 (2016), publisher: American Chemical Society.
  * Sureka  _et al._ (2019) H. V. Sureka, A. C. Obermeyer, R. J. Flores, and B. D. Olsen, ACS Applied Materials & Interfaces 11, 32354 (2019), publisher: American Chemical Society.
  * Shi and Weitz (2017) W. Shi and D. A. Weitz, Macromolecules 50, 7681 (2017), publisher: American Chemical Society.
  * Haase  _et al._ (2017) M. F. Haase, H. Jeon, N. Hough, J. H. Kim, K. J. Stebe, and D. Lee, Nature communications 8, 1234 (2017).
  * Tree  _et al._ (2017) D. R. Tree, K. T. Delaney, H. D. Ceniceros, T. Iwama, and G. H. Fredrickson, Soft Matter 13, 3013 (2017).
  * Wei  _et al._ (2017) M.-T. Wei, S. Elbaum-Garfinkle, A. S. Holehouse, C. C.-H. Chen, M. Feric, C. B. Arnold, R. D. Priestley, R. V. Pappu, and C. P. Brangwynne, Nature Chemistry 9, 1118 (2017), number: 11 Publisher: Nature Publishing Group.
  * Banani  _et al._ (2017) S. F. Banani, H. O. Lee, A. A. Hyman, and M. K. Rosen, Nature Reviews Molecular Cell Biology 18, 285 (2017), number: 5 Publisher: Nature Publishing Group.
  * Kanaan  _et al._ (2020) N. M. Kanaan, C. Hamel, T. Grabinski, and B. Combs, Nature Communications 11, 2809 (2020), number: 1 Publisher: Nature Publishing Group.
  * Choi  _et al._ (2020) J.-M. Choi, A. S. Holehouse, and R. V. Pappu, Annual Review of Biophysics 49, 107 (2020), _eprint: https://doi.org/10.1146/annurev-biophys-121219-081629.
  * Thompson  _et al._ (2022) A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Computer Physics Communications 271, 108171 (2022).
  * Phillips (2022) J. C. Phillips, “NAMD - Scalable Molecular Dynamics,” (2022).
  * Abraham  _et al._ (2015) M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindahl, SoftwareX 1-2, 19 (2015).
  * Fredrickson (2006) G. Fredrickson,  _The equilibrium theory of inhomogeneous polymers_ , 134 (Oxford University Press on Demand, 2006).
  * Edwards (1965) S. F. Edwards, Proceedings of the Physical Society 85, 613 (1965).
  * Leibler (1980) L. Leibler, Macromolecules 13, 1602 (1980).
  * Zheng and Wang (1995) W. Zheng and Z.-G. Wang, Macromolecules 28, 7215 (1995).
  * Milner  _et al._ (1988) S. T. Milner, T. A. Witten, and M. E. Cates, Macromolecules 21, 2610 (1988).
  * Edwards and Freed (1970) S. Edwards and K. Freed, Journal of Physics C: Solid State Physics 3, 739 (1970).
  * Koski  _et al._ (2013) J. Koski, H. Chao, and R. A. Riggleman, The Journal of chemical physics 139, 244911 (2013).
  * Detcheverry  _et al._ (2008) F. A. Detcheverry, H. Kang, K. C. Daoulas, M. Müller, P. F. Nealey, and J. J. De Pablo, Macromolecules 41, 4989 (2008).
  * Pike  _et al._ (2009) D. Q. Pike, F. A. Detcheverry, M. Müller, and J. J. de Pablo, The Journal of chemical physics 131, 084903 (2009).
  * Daoulas and Müller (2006) K. C. Daoulas and M. Müller, The Journal of chemical physics 125, 184904 (2006).
  * Ganesan and Pryamitsyn (2003) V. Ganesan and V. Pryamitsyn, The Journal of chemical physics 118, 4345 (2003).
  * Chao  _et al._ (2017a) H. Chao, J. Koski, and R. A. Riggleman, Soft Matter 13, 239 (2017a).
  * Schneider and Müller (2019) L. Schneider and M. Müller, Computer Physics Communications 235, 463 (2019).
  * Arora  _et al._ (2016) A. Arora, J. Qin, D. C. Morse, K. T. Delaney, G. H. Fredrickson, F. S. Bates, and K. D. Dorfman, Macromolecules 49, 4675 (2016).
  * Cheong  _et al._ (2020) G. K. Cheong, A. Chawla, D. C. Morse, and K. D. Dorfman, The European Physical Journal E 43, 1 (2020).
  * Chao  _et al._ (2017b) H. Chao, J. Koski, and R. A. Riggleman, Soft Matter 13, 239 (2017b).
  * Lab (2023) G. Lab, “GSD 2.7.0 documentation,” (2023).
  * Helfand (1975a) E. Helfand, The Journal of Chemical Physics 62, 999 (1975a).
  * Villet and Fredrickson (2014) M. C. Villet and G. H. Fredrickson, The Journal of Chemical Physics 141, 224115 (2014), publisher: American Institute of Physics.
  * Weyman  _et al._ (2021) A. Weyman, V. G. Mavrantzas, and H. C. Öttinger, The Journal of Chemical Physics 155, 024106 (2021).
  * Helfand (1975b) E. Helfand, The Journal of chemical physics 62, 999 (1975b).
  * Allen and Tildesley (2017) M. P. Allen and D. J. Tildesley,  _Computer simulation of liquids_ (Oxford university press, 2017).
  * Parisi (1983) G. Parisi, Physics Letters B 131, 393 (1983).
  * Klauder (1984) J. R. Klauder, Physical Review A 29, 2036 (1984).
  * Stasiak and Matsen (2013) P. Stasiak and M. Matsen, Macromolecules 46, 8037 (2013).
  * Lennon  _et al._ (2008) E. M. Lennon, G. O. Mohler, H. D. Ceniceros, C. J. García-Cervera, and G. H. Fredrickson, Multiscale Modeling & Simulation 6, 1347 (2008).
  * Bell and Hoberock (2012) N. Bell and J. Hoberock, in  _GPU computing gems Jade edition_ (Elsevier, 2012) pp. 359–371.
  * Deserno and Holm (1998) M. Deserno and C. Holm, The Journal of Chemical Physics 109, 7678 (1998), publisher: American Institute of Physics.
  * McMillan (1971) W. L. McMillan, Physical Review A 4, 1238 (1971).
  * Pryamitsyn and Ganesan (2004) V. Pryamitsyn and V. Ganesan, The Journal of chemical physics 120, 5824 (2004).
  * Groot  _et al._ (1999) R. D. Groot, T. J. Madden, and D. J. Tildesley, The Journal of Chemical Physics 110, 9739 (1999), publisher: American Institute of Physics.
  * Zhang (2016) K. Zhang, “On the Concept of Static Structure Factor,” (2016), arXiv:1606.03610 [cond-mat].
  * Ceniceros and Fredrickson (2004) H. D. Ceniceros and G. H. Fredrickson, Multiscale Modeling & Simulation 2, 452 (2004).
  * Riggleman  _et al._ (2012) R. A. Riggleman, R. Kumar, and G. H. Fredrickson, The Journal of chemical physics 136, 024903 (2012).
  * Ramírez-Hernández  _et al._ (2015) A. Ramírez-Hernández, B. L. Peters, M. Andreev, J. D. Schieber, and J. J. de Pablo, The Journal of chemical physics 143, 243147 (2015).
  * Ramírez-Hernández  _et al._ (2017) A. Ramírez-Hernández, B. L. Peters, L. Schneider, M. Andreev, J. D. Schieber, M. Müller, and J. J. de Pablo, The Journal of chemical physics 146, 014903 (2017).
  * Martin  _et al._ (2016) J. M. Martin, W. Li, K. T. Delaney, and G. H. Fredrickson, The Journal of Chemical Physics 145, 154104 (2016), publisher: AIP Publishing LLCAIP Publishing.
  * Chao  _et al._ (2014) H. Chao, B. A. Hagberg, and R. A. Riggleman, Soft Matter 10, 8083 (2014).
  * Koski  _et al._ (2015) J. Koski, H. Chao, and R. A. Riggleman, Chemical Communications 51, 5440 (2015).
  * Sides  _et al._ (2006) S. W. Sides, B. J. Kim, E. J. Kramer, and G. H. Fredrickson, Physical review letters 96, 250601 (2006).
  * Stukowski (2010) A. Stukowski, MODELLING AND SIMULATION IN MATERIALS SCIENCE AND ENGINEERING 18 (2010), 10.1088/0965-0393/18/1/015012.