Tutorial
This tutorial presents the basic usage of Polyorder.jl.
For more examples and details about the usage of Polyorder.jl, consult the User Guide section of the documentation and the testing codes in test folder.
Construct an SCFT model
To use Polyorder, first we create an SCFT model:
julia> using Polymer # for AB_system
julia> using Scattering # for BravaisLattice, UnitCell
julia> using Polyorder
# Create a polymer system with a symmetric AB diblock copolymer: fA=0.5, χN=20.0.
julia> ab = AB_system();
# Create a 1D lattice with size 4.0 Rg.
julia> lat = BravaisLattice(UnitCell(4.0));
# Create a scft model with same Δs for A and B blocks.
julia> scft = NoncyclicChainSCFT(ab, lat, 0.01);Non-orthogonal unit cell is supported (requires v0.10.0+)
# Create a 2D Hexagonal unit cell with side length 4.0 Rg.
julia> uc = UnitCell(Hexagonal2D(), 4.0)
# Create a Bravais Lattice based on the unit cell.
# We can also specify the space group by give an ID.
# Here 17 means 2D P6mm space group.
julia> lat = BravaisLattice(uc, 17);Mixed contour steps (ds) and MDE solvers (mde) are supported (requires v0.5.0+). Since v0.17.0, the mixed contour steps and MDE solvers can be specified for each block separately using block labels. Meanwhile, ds and mde provided for each component as a vector is not allowed anymore. A single real number or a dict with label=>ds pair can be mixed for each component. A single MDEAlogrithm type or a dict with label=>MDEAlgorithm pair can be mixed for each component.
# Even there is only one component, it should be wrapped in a Dict.
julia> ds = Dict(:AB => Dict(:A=>0.01, :B=>0.02))
julia> mde = Dict(:AB => Dict(:A=>OSF, :B=>RQM4))
# mixed contour steps: ds = 0.01 for A block, ds = 0.02 for B block
julia> scft1 = NoncyclicChainSCFT(ab, lat, ds);
# mixed MDE solvers: OSF for A block, RQM4 for B block
julia> scft2 = NoncyclicChainSCFT(ab, lat, 0.01; mde);
# mixed ds and MDE
julia> scft3 = NoncyclicChainSCFT(ab, lat, ds; mde)Since v0.24.0, the following API is available:
julia> scft = NoncyclicChainSCFT(ab, lat; ds, mde)Also starting from v0.24.0, a new type MDESolver is introduced. It combines contour steps and MDE solvers. In addition, an extra options field of MDESolver is used to provide extra options when instantiating the MDE solver with mdesolver.algo(q, w; mdesolver.options...), where mdesolver is a MDESolver. We can construct the same SCFT model as follows:
julia> mdesolvers = Dict(:AB => Dict(:A=>MDESolver(OSF, 0.01), :B=>MDESolver(RQM4, 0.02)))
julia> scft = NoncyclicChainSCFT(ab, lat; mdesolvers)
# Each block of each component share the same MDE solver:
julia> scft = NoncyclicChainSCFT(ab, lat, MDESolver(OSF, 0.01))Solve the SCFT model
Then we can solve the SCFT model with current simulation cell:
julia> Polyorder.solve!(scft)We can control how Polyorder handles the solution process by using Polyorder.Config and SCFTConfig:
julia> scft = NoncyclicChainSCFT(ab, lat, 0.01; updater=PicardMann(0.4));
# Change the default tolerance form 1e-5 to 1e-4, default maximum number of iterations from 2000 to 1000.
julia> scft_config = SCFTConfig(tol=1e-4, max_iter=1000)
# Wrap scft config into Polyorder.Config.
julia> config = Polyorder.Config(scft=scft_config)
# Apply the new configurations
julia> Polyorder.solve!(scft, config)One can provide a new updater rather than scft.updater to drive the actual simulation
julia> updater = AndersonSD(0.2; m=10, warmup=0)
julia> Polyorder.solve!(scft, updater)Many updaters for solving SCFT equations are available, for example: PicardMann, Nesterov, ARDM, NGMRES_S1, Anderson_S1, NGMRES, OACCEL, Anderson (v0.6.0+), SIS, ETD, PO, ETDPEC (v0.7.0+), and SISF, ETDF (v0.7.2+).
# Anderson can be preconditioned by any qualified preconditioner.
julia> precond = ETD(0.4)
julia> updater = Anderson(precond; m=15, αw=0.4, warmup=50)
# Note that we can still provide a `config` object to control the simulation.
julia> Polyorder.solve!(scft, updater, config)The simulation can be terminated by a combination of controls based on IterationControl.
A progress meter based on ProgressMeter.jl can be enabled by setting config.io.progress_solve to be true. (v0.12.0+)
Cell optimization
There are three ways to optimize the simulation cell:
- Variable cell optimization. It updates auxiliary fields and cell simultaneously. It is the most efficient way.
- Stress-guided cell optimization. It uses optimization methods to find the zero-stress unit cell with the aid of stress, i.e. the gradient with respect to cell size and shape.
- Gradient-free cell optimization. It optimizes the unit cell sizes using gradient-free optimization methods. Only orthogonal unit cells are supported.
Variable cell optimization (v0.13.0+)
The variable cell optimization is the recommended way to optimize the cell as long as a good initial guess of the stable phase structure is available. The stress optimization and free energy optimization are performed simultaneously which greatly improved the efficiency.
To invoke variable cell optimization, either NoncyclicChainSCFT object has an VariableCell updater or pass a VariableCell updater to cell_solve!.
# `scft.updater` is a `VariableCell` object.
julia> cell_solve!(scft)
# Provide a `VariableCell` updater explicitly.
julia> updater = VariableCell(BB(1.0), SD(0.2))
julia> cell_solve!(scft, updater)
# An `VariableCellOpt` object can be provided explicitly to invoke the variable cell optimization.
julia> opt = VariableCellOpt(updater)
julia> cell_solve!(opt, scft)
# Provide a `Polyorder.Config` object to further control the cell optimization and individual fixed-cell SCFT simulations.
julia> cell_solve!(scft, config)
julia> cell_solve!(scft, updater, config)
julia> cell_solve!(opt, scft, config)Note that scft is being modified in place.
Stress-guided cell optimization (v0.11.0+)
This is an compromising approach which utilizes the stress information to guide the optimization of the simulation cell. The stress optimization and free energy optimization are performed alternately.
if neither scft.updater nor updater is a VariableCell updater, the stress-guided cell optimization is invoked with a default OptimStressGuidedCellOpt object.
Note that to use OptimStressGuidedCellOpt, the Optim.jl package should be installed and loaded first.
pkg> add Optim
julia> using Optim
# `scft.updater` should not be a `VariableCell` object.
julia> cell_solve!(scft)
# Provide an updater which should not be a `VariableCell` object.
julia> updater = SD(0.2)
julia> cell_solve!(scft, updater)
# A `OptimStressGuidedCellOpt` object can be provided explicitly to invoke the stress-guided cell optimization.
julia> algo = Optim.LBFGS() # or any other Optim algorithm
julia> options = (; ) # options supported by Optim to control the optimization
julia> opt = OptimStressGuidedCellOpt(algo, options)
julia> cell_solve!(opt, scft) # or `cell_solve!(opt, scft, updater)`
# Provide a `Polyorder.Config` object to further control the cell optimization and individual fixed-cell SCFT simulations.
julia> cell_solve!(scft, config)
julia> cell_solve!(scft, updater, config)
julia> cell_solve!(opt, scft, config)
julia> cell_solve!(opt, scft, updater, config)Gradient-free cell optimization
This is a legacy method for optimizing the simulation cell. It can only handle orthogonal simulation cells. The stress of the simulation cell is not computed.
Gradient-free cell optimization should be invoked explicitly with a GradientFreeCellOptAlgorithm object. Currently, There is only one type of such algorithm based on Optim.jl: OptimGradientFreeCellOpt.
Note that to use OptimGradientFreeCellOpt, the Optim.jl package should be installed and loaded first.
julia> using Optim
# Provide a `OptimGradientFreeCellOpt` object explicitly to invoke the gradient-free cell optimization.
julia> algo = Optim.Brent() # for 1D, or Optim.NelderMead() for higher dimensions
julia> options = (; ) # options supported by Optim to control the optimization
julia> opt = OptimGradientFreeCellOpt(algo, options)
# `scft.updater` should not be a `VariableCell` object.
julia> cell_solve!(opt, scft)
# An `updater` other than a `VariableCell` object can be provided explicitly.
julia> cell_solve!(opt, scft, updater)
# Provide a `Polyorder.Config` object to further control the cell optimization and individual fixed-cell SCFT simulations.
julia> cell_solve!(opt, scft, config)
julia> cell_solve!(opt, scft, updater, config)Good practice for cell optimization
A good practice is using solve! to obtain a well converged model, then serve this model to cell_solve!.
Since v0.12.0, there are 4 files: config.yml, summary.yml, tracecell, and tracesolve are saved to the path config.io.base_dir. By default, two additional files: fields.h5, densities.h5 are saved. It can be omitted by setting config.io.save_w and/or config.io.save_ϕ to be false. The available data format can be either HDF5 (built-in) or MAT (should install and load MAT.jl first) by setting config.io.data_format accordingly.
A progress meter for cell_solve! based on ProgressMeter.jl can be enabled by setting config.io.progress_cell to be true. (v0.12.0+)
Visualization
Polyorder provides a function to visualize the density fields of in a simulation cell of a SCFT model. Makie and one of its backend should be installed and loaded to use this function.
julia> using GLMakie
julia> densityplot(scft)Convergence traces of an SCFT simulation can be visualized by the traceplot function:
julia> traceplot(scft)
# or
julia> traceplot(updater)By default, the residual error is plotted. Other traces are also available. Please see documentation for details.
Serialization and configurations
We can save the full set of configurations for one scft model calculation, including the polymer system and the lattice (simulation cell) to the YAML file
# Serialize the scft model to a `Polyorder.Config` object.
julia> config = to_config(scft)
# Save the config to a YAML file.
julia> save_config("./config.yml", config)The YAML file can be loaded again as a Polyorder.Config object.
# Note that the second argument is mandatory.
julia> config_loaded = load_config("./config.yml", Polyorder.Config)We can construct PolymerSystem, BravaisLattice, and NoncyclicChainSCFT objects from PolymerSystemConfig, LatticeConfig and Polyorder.Config as:
julia> system = from_config(config_loaded.system)
julia> lattice = from_config(config_loaded.lattice)
julia> scft = from_config(config_loaded) # or NoncyclicChainSCFT(config_loaded)