Algorithms

Partial Saddle-Point Approximation (Prev)         Parameter File (Next)        

This page discusses several types of algorithm that are used to perform and analyze PS-FTS calculations, but that are not needed for SCFT calculations. We focus here on explaining the different types of algorithm that are required to perform a PS-FTS calculation, and on listing the available choices among sets of algorithms that serve the same or similar purposes (such as different Brownian dynamics step algorithms). The following page focuses instead on describing the format of parameter file required to initialize a PS-FTS calculation, including the sections that specify choices of algorithms.

Types of Algorithm

The pscf_pc and pscf_pg programs can perform PS-FTS calculations using either a Brownian dynamics (BD) or Monte Carlo (MC) method to stochastically sample a fluctuating field or fields. BD and MC sampling algorithms for PS-FTS have a similar structure. A single BD step or MC move generally involves generation of a random change in one or more exchange fields, followed by use of an iterative algorithm to identify a new partial-saddle point field configuration. Throughout the PSCF documentation and source code, iterative algorithms that are used to identify partial saddle-point configurations, and thus satisfy the mean-field incompressibility constraint, are referred to as "compressor" algorithms.

For example, consider a PS-FTS simulation of a standard AB system with two monomer types and a positive Flory-Huggins interaction parameter. Such a system has a real exchange field \( W_{-}({\bf r}) = \omega_{-}({\bf r}) \) that undergoes random fluctuations and a real pressure-like field \( W_{+}({\bf r}) = i \omega_{+}({\bf r})\) that must be chosen to satisfy the partial saddle-point condition. At the beginning of either a BD step or MC move for such as system, a small random change is added to the exchange field \( W_{-} \). Next, an iterative compressor algorithm is used to adjust the presssure-like field \( W_{+}({\bf r}) \) so as to satisfy the partial saddle-point condition, thus imposing incompressibility at a mean-field level. Finally, for an MC simulation, a decision is made about whether to accept or reject the resulting MC move, based on the value of the Hamiltonian in the new partial saddle-point configuration.

PSCF also provides tools for analysis of simulation data. Users may use the parameter file to enable one or more of a set classes that can periodically compute specific physical quantitites, output data to a file, and/or perform statistical analyses during a simulation. We refer to C++ classes that perform these sorts of periodic data output and analysis operations "analyzer" classes.

The parameter file for a typical PS-FTS calculation will thus normally contain a choice of algorithms for the following three purposes:

• Sampling : A BD step algorithm or a set of MC step algorithms that is used to update the fluctuating field components.
• Compressor : A compressor algorithm that is used to adjust the pressure-like field \( W_{+} \) so as to always satisfy the partial saddle-point condition.
• Analyzers : File IO and/or data analysis operations that are performed at regular intervals during a simulation.

In addition to these basic building blocks, users may also enable two other types of optional elements:

• Ramp : A "ramp" is a optional feature that allows a user to perform a simulation in which one or more parameters of the simulation model are varied continuously over the course of a simulation.
• Perturbation : A "perturbation" is an optional feature that allows the user to modify the form of the field theoretic Hamiltonian by adding a perturbation to the standard form

All of these features are discussed separately below.

In the source code of PSCF, different algorithms of the same type are generally implemented by C++ classes that are subclasses of a common base class. The generic names that we use in this documentation to describe different types of algorithm generally correspond to the names of the corresponding base classes. Specifically:

• Brownian dynamics step algorithms are implemented by subclasses of a base class named BdStep
• Monte Carlo moves are implemented by subclasses of McMove
• File Io and data analysis operations are implemented by subclasses of Analyzer.
• Algorithms for continuous variation of parameters are implemented by subclasses of Ramp
• Classes that define additive changes to the field theoretic Hamiltonian are subclasses of Perturbation

The PSCF parameter file format uses selectable parameter file blocks to allow a user to either choose one from a set of different algorithms of the same type (as is appropriate for BdStep an Compressor algorithms) or to choose a set one or more algorithms from among a list of available algorithms of the same type (as is appropriate for McMove and Analyzer algorithms).

Brownian Dynamics (BD) Step Algorithms

A Brownian dynamics (BD) algorithm for a particle based simulation is usually based on discretization of a stochastic differential equation that generates diffusive dynamics. In the context of a field theoretic simulation, however, there is no natural way to define dynamics or units of time, since the BD algorithm is simply being used as a way to sample an equilibrium distribution. The description of BD step algorithms used to define parameters in PSCF based on units in which the fictitious "time" is thus measured in units of BD steps, while energy is measured in thermal units, in which \( kT = 1 \). The magnitude of changes in the exchange field(s) per time step is controlled by a parameter that is denoted by the label "mobility" in the parameter file, and denoted by the symbol \( \gamma \) in these notes.

As an example, consider the explicit Euler BD step algorithm for a standard AB system. This algorithm yields a change \( \Delta W_{-}({\bf r}) \) per time step in the exchange field \( W_{-}({\bf r}) \) at grid point \( {\bf r} \) that is given by

\[ \Delta W_{-}({\bf r}) = - \gamma \frac{\delta H}{\delta W_{-}({\bf r})} + g({\bf r}) \]

where \( \gamma \) is the BD mobility parameter, \( g({\bf r}) \) is a random field displacement, and

\[ \frac{\delta H}{\delta W_{-}({\bf r})} = \frac{1}{v} \left [ \frac{2 W_{-}({\bf r})}{\chi_{b}} - \Phi_{-}({\bf r}) \right ] \quad. \]

The random displacement \( g({\bf r}) \) is chosen independently at each node of the computational mesh and each BD step from a Gaussian distribution with a vanishing mean value, \( \overline{g({\bf r})} = 0\), and with a variance

\[ \overline{g^{2}({\bf r})} = \gamma \Delta V \quad. \]

Here, \( \Delta V \) is the volume per node (or grid point) of the computational mesh, given by \( \Delta V = V/G \), where \( V \) is the total system volume and \( G \) is the total number of nodes in the mesh. Other authors have sometimes referred to what we call a "mobility" as a time step, sometimes indicated by a symbol such \( \Delta t \) or \( \delta t \).

The mobility parameter has an analogous meaning in more sophisticated BD step algorithms. It always appears as a prefactor of the "force" (or functional derivative) \( \delta H/\delta W_{-}({\bf r}) \) in expressions for changes in \( W_{-}({\bf r}) \) over a single time step.

PSCF currently allows users to choose from among three BD step algorithms, which are labelled below by the names of the C++ classes that implement them:

• LMBdStep : A "Leimkuhler-Matthews" (LM) BD step algorithm
• PredCorrBdStep: A predictor corrector algorithm,
• ExplicitBdStep : An explicit Euler algorithm like that described above,

More details about each of these algorithms can be obtained by clicking on the associated link.

We recommend the Leimkuhler-Matthews (LM) algorithm for general use in all BD simulations. The LM step algorithm is found to provide much more accurate results than the competing predictor-corrector algorithm, allowing the use of larger time steps. The explicit Euler algorithm is signficantly less accurate than the other two alternatives.

Monte Carlo (MC) Move Algorithms

The PSCF implementation of MC sampling allows the user to specify a set of Monte Carlo move types to be used during an MC simulation. When more than one type of move is specified, one of the move types is chosen at random at the beginning of each step, using probabilities that are specified by the user. The algorithm for each attempted step of a MC simulation is thus as follows:

• Choose a type of MC move at random from among a set of choices specified in the parameter file.
• Use the chosen MC move to generate a random change in \( W_{-} \).
• Use the chosen compressor algorithm to adjust the real pressure-like field \( W_{+}({\bf r}) \) so as to satisfy the partial saddle-point condition.
• Evaluate the Hamiltonian in the new partial saddle point state.
• Use a Metropolis criterion that is designed to satisfy detailed balance to decide whether to accept or reject the move.
• If the move is rejected, restore the field configuration that existed at the beginning of the attempted move.

The block of the parameter file for an MC simulation that contains information about MC moves, which is labelled "McMoveManager", contains one or more subblocks that are associated with different types of move. The parameter file block associated with each such MC move type contains a parameter named "probability" that specifies the probability of that type of move being chosen at the start of each step.

PSCF currently provides the following two MC step algorithms:

• ForceBiasMove: A "smart MC" or "force bias" MC move that uses an explicit Euler BD step to generate a proposed MC move.
• RealMcMove : A simple MC move that generates spatially uncorrelated random changes of \( W_{-} \) at different grid points

The ForceBiasMove algorithm generally leads to much more efficient sampling than the simpler real move, and is recommended as the workhorse method for MC simulations.

Acceptance rates for MC moves generally decrease rapidly with increasing system size. This effect is particularly severe with the simple real-space move (RealMcMove), but also occurs for the force-bias move. Because BD simulation does not suffer from any analogous problem, this phenomena causes MC to always become less efficient than BD for sufficiently large systems. We thus generally recommend the use of BD rather than MC simulation for 3D systems with roughly \( 10^{4} \) or more grid points.

One virtue of MC over BD simulation is that, because MC algorithms are designed to satisfy a detailed balance condition, they do not exhibit any systematic sampling error. An infinitely long MC simulation would, in in the absence of round-off error, yield exact results for any average value obtained from a particular model on a specified computational mesh. The only sampling errors in MC simulations are true statistical errors arising from the finite length of a simulation, which can be estimated by standard methods, and errors arising from round off of floating point numbers, which are usually negligible.

BD simulations instead exhibit a systematic error arising from discretization of a time-like coordinate in a diffusion model. The magnitude of this error increases with increases in the mobility parameter. Results of MC simulations of small systems can thus sometimes be useful for comparison when testing new BD algorithms, and when evaluating systematic errors in BD algorithms. Systematic errors in BD simulations can also be estimated by running simulations with several different values of the mobility parameter.

Compressor Algorithms

The PSCF paramater file format allows the user to choose from among several "compressor" algorithms. As discussed above, a compressor is an iterative algorithm that adjusts the \( W_{+} \) field so as to satisfy a partial saddle-point condition. This algorithm must applied at last once per BD or MC stochastic step, after each change in \( W_{-} \). Users may choose a compressor algorithm independently of their choice for a BD step or MC move algorithms: Any compressor algorithm may be used in the context of any BD or MC step algorithm.

Solution of the saddle-point condition normally requires multiple iterations of \( W_{+} \). Each iteration normally involves a proposed adjustment of \( W_{+} \) followed by solution of the modified diffusion equation (MDE) in the resulting state, which is necessary to evaluate the remaining error. Because repeated solution of the MDE is the most expensive operation in a FTS, and the compressor algorithm requires multiple solutions of the MDE for each BD or MC step, almost all of the computational time of a PS-FTS calculation is actually spent within the iterative loop of a compressor algorithm.

The available compressor algorithms are listed below. Users may click on the link associated with each algorithm label to access more a detailed description of each algorithm, and of the format of the associated parameter file block.

• LrAmCompressor: A linear-response Anderson mixing algorithm
• LrCompressor: A simple linear-response algorithm
• AmCompressor: An Anderson-mixing algorithm

The LrAmCompressor is the default choice, and is recommended for general use.

The AmCompressor algorithm is similar to the Anderson mixing algorithm used in most previously published PS-FTS calculations. This algorithm uses information from previous iterations to guide the choice of a new field configuration during each iteration after the first.

The LrAmCompressor and LrCompressor algorithms both use a quasi-Newton approximation that uses the analytically calculated linear response of a hypothetical homogeneous state as an approximation to the actual linear response of the monomer concentration in an inhomgeneous state to small changes in the pressure-like field. The Lr algorithm is a straightforward quasi-Newton algorithm based on this approximation. The LrAm uses this quasi-Newton step as one stage of a modified Anderson mixing algorithm, which also uses information from previous steps.

Analyzer Algorithms

PSCF provides an extensible set of C++ classes for analysis of field theoretic simulations. These are referred to here as analyzers classes, or simply analyzers. Each analyzer class performs a specified analysis and/or file IO operation at some regular interval during a simulation (e.g., every 10 or 100 simulation steps). Analyzer classes can be designed to periodically compute physical quantities of interest, output data to file, or compute statistical properties such as averages and variances, or any combination of these types of operation. The PSCF parameter file format allows the user to enables a set of zero or more analyzers in order to specify which such operations should be performed during a simulation that is run with the SIMULATE command, and how often each should be performed. The same set of analyzer classes may also be used during a postprocessing operation that is performed by the ANALYZE command, as discussed in more detail here.

All of the analyzer classes that are provided with the current version of PSCF are listed below. A brief description of the purpose of each class is provided below, and further details can be obtained by clicking on the class name.

• StepLogger: Record progress of a simulation by periodically writing the number of completed BD or MC steps to standard output (e.g., the screen).
• TrajectoryWriter : Periodically write the current configuration of monomer chemical potential fields (w-fields) to field trajectory file in r-grid format, for later postprocessing.
• ConcentrationWriter : Periodically write all monomer concentration fields (c-fields) to a file, in r-grid format.
• HamiltonianAnalyzer : Compute average value of the field theoretic Hamiltonian, and optionally output Hamiltonian values.
• BinaryStructureFactorGrid: Compute the structure factor of a system with two monomer types.
• MaxOrderParameter : Compute an order parameter that is useful for identifying spontaneous phase transitions.
• FourthOrderParameter: Compute another order parameter that is useful for identifying spontaneous phase transitions.
• ChiDerivative: Compute the derivative of free energy with respect to changes in a chi parameter, for use in thermodynamic integration.
• ConcentrationDerivative Compute the derivative of free energy with respect to changes in total monomer concentation (the inverse of vMonomer), for use in thermodynamic integration.
• PerturbationDerivative: Compute the derivative of free energy with respect to the strength of a perturbation (the lambda parameter), for use in thermodynamic integration.

Note that the first three classes listed above (StepLogger, TrajectoryWriter, and ConcentrationWriter) simply write readily available information to a file or standard output (e.g., the screen) at regular intervals, rather than doing any nontrivial data analysis.

The analyzers named MaxOrderParameter and FourthOrderParameter) are designed specifically as tools for identifying spontaneous phase transitions. Specifically, these analyzers may be used to analyze a parameter ramp that is designed to induce a spotaneous disorder-to-order (crystallization) or order-to-disorder (melting) transition.

The last three analyzers listed above (ChiDerivative, ConcentrationDerivative, and PerturbationDerivative) are designed specifically for use in thermodynamic integration calculations of free energy differences, which are discussed in the next section. Specifically, a PerturbationDerivative analyzer may be used to implement Einstein crystal free energy calculations.

Thermodynamic Integration

Thermodynamic integration is a technique for computing a difference in values of the free energy in two states characterized by different values of a parameter. This free energy difference is computed by numerically integrating the derivative of the free energy \( F \) with respect to the relevant model parameter. Suppose \( \lambda \) is a parameter that appears in the field theoretic Hamiltonian \( H \). The method is based on a theorem which states that the derivative of free energy \( F \) with respect to \( \lambda \) is given by an equilibrium average value

\[ \frac{\partial F}{\partial \lambda} = \left \langle \frac{\partial H}{\partial \lambda} \right \rangle \quad, \]

where \( \partial H/\partial \lambda \) is a derivative evaluated in a fixed field configuration, and \( \langle \cdots \rangle \) denotes an equilibrium ensemble average evaluated using the specified value of \( \lambda \).

The change in \( F \) associated with a change in \( \lambda \) from an initial value \( \lambda_{0} \) to a final value \( \lambda_{1} \) is obtained by integrating the derivative \( \partial F/ \partial \lambda \) over this range. This sort of integral can be computed by either of two methods:

• Separate simulations : Equilibrium simulations can be performed at a set of values of \( \lambda \) within the range of interest, each of which yields a value of \( \partial F/\partial \lambda \) at a single value as an ensemble average. The required integral is then approximated by numerical integration of the resulting sequence of values.
• Continuous integration : The integral may be computed by using a ramp to continuously vary the parameter of interest over the course of a simulation.

Continuous integration may be used either to compute the change in free energy associated with a change in a parameter of the standard field theoretic Hamiltonian (such as an interaction parameter or the monomer concentration) or as change in free energy associated with changes in the strength of a perturbation that modifies the form of the Hamiltonian (as done in Einstein crystal integration).

Parameter Ramp

A "ramp" is an optional feature for FTS that allows a user to continuously vary one or more model parameters during a simulation. Parameter ramps are discussed in more detail on a separate manual page.

The only ramp algorithm currently provided by PSCF is implemented by the class LinearRamp. The LinearRamp class performs a ramp in which one or more of the parameters of a model are taken to vary linearly with the BD or MC simulation step counter during a simulation. The interface for a LinearRamp is closely analogous to that of the LinearSweep class that is used for SCFT sweep calculations, as discussed here.

Perturbation

A "perturbation" is an optional feature that allows a user to modify the form of the field theoretic Hamiltonian by adding an arbitrary additional functional of the w-fields to the standard Hamiltonian. Perturbations are discussed on a separate manual page.

The only type of perturbation that is available in the current version of PSCF is implemented by a class named EinsteinCrystalPerturbation, which is designed to enable the Einstein crystal method of computing free energies. The parameter file format is discussed in detail on a separate page.

Partial Saddle-Point Approximation (Prev)         Partial Saddle-Point Field Theoretic Simulation (PS-FTS) (Up)         Parameter File (Next)