Phase Behaviors of Block Copolymers via Polyorder.jl and Related Packages - Part 1

Quick Setup

Installing Packages

Currently, Polyorder.jl and its dependent packages are private repositories hosted at https://github.com/liuyxpp, which are not registered in the Julia General Registry:

• PolymerArchitechture.jl
• Scattering.jl
• PhaseDiagram.jl

It is recommended to install them via the private registry. First, you need to add the private registry to the Julia package manager in REPL:

(@v1.10) pkg> registry add git@github.com:liuyxpp/lab3083.git

or

(@v1.10) pkg> registry add https://github.com/liuyxpp/lab3083.git

This step only has to be done once. You can now install packages in both the General registry and the lab3083 registry. The dependent packages will be installed automatically just as installing packages normally.

Install all following packages to your project environment via:

(@v1.10) pkg> activate .
(your_project_dir) pkg> add Polymer PolymerArchitecture PolymerArchitectureMakie Scattering Polyorder PhaseDiagram

The environment your_project_dir can be any name you like. Make sure that

• Polymer.jl (v0.8.11+)
• PolymerArchitecture (v0.3.6+)
• Scattering.jl (v0.4.4+)
• Polyorder.jl (v0.23.0+)
• PhaseDiagram.jl (v0.8.2+)
• PolymerArchitectureMakie (v0.3.0+)

which can be checked via

(your_project_dir) pkg> st

To perform stress-guided cell optimization and compute stability limit via RPA, two extra packages are required:

(your_project_dir) pkg> add Optim Roots

Other packages required by this tutorial are registered in the General Registry, and can be installed via

(your_project_dir) pkg> add CairoMakie LaTeXStrings

Loading Packages

In the notebook, we first activate the environment reside in the folder you start Pluto. This step will disable the Pluto built-in package manager. If you are not in Pluto notebook, you should use an appropriate way to activate current environment.

import Pkg
Pkg.activate(".")

using Polymer
using Scattering
using Polyorder
import Polyorder as JP
using PhaseDiagram
using PolymerArchitectureMakie
using CairoMakie
using LaTeXStrings

  Activating project at `~/SynologyDrive/Develop/Polyorder.jl/docs/tutorial`

# Suppressing displaying source code file location for @info
using Logging
global_logger(ConsoleLogger())

ConsoleLogger(IJulia.IJuliaStdio{Base.PipeEndpoint}(IOContext(Base.PipeEndpoint(RawFD(43) open, 0 bytes waiting))), Info, Logging.default_metafmt, true, 0, Dict{Any, Int64}())

include("common.jl")

plot_density (generic function with 3 methods)

Create Polymer Systems

Polyorder performs computations on a PolymerSystem object which is created by the Polymer.jl package.

A PolymerSystem object contains BlockCopolymer and/or SmallMolecule objects.

Polymers and small molecules

Small molecules

Small molecules can be created by

SmallMolecule(:S)  # :S is the label/name of the samll molecule.

SmallMolecule S with b=1.0

Or you can use the handy method

solvent()  # the default label is :S

SmallMolecule S with b=1.0

Polymers

Polymer.jl currently only supports creating neat polymers and block copolymers, which are represented by a BlockCopolymer object. BlockCopolymer internally uses a graph to represent the chain architecture.

In Polymer.jl, a BlockCopolymer are consisted of PolymerBlocks, each of which is consisted of KuhnSegments, and they are connected by BranchPoints. They may also have FreeEnds.

A method for creating following polymer chain

[A1] - [B] - [A2] - [C1]
         |      |
         [C2]   [D]

Note that the two brackets around each block stand for two ends of that block.

function chainABCACD()
    # The chain has four distinct Kuhnsegments.
    sA = KuhnSegment(:A)
    sB = KuhnSegment(:B)
    sC = KuhnSegment(:C)
    sD = KuhnSegment(:D)
    # Among all block ends, there are three branch points and four free ends.
    eb1 = BranchPoint(:EB1)
    eb2 = BranchPoint(:EB2)
    eb3 = BranchPoint(:EB3)
    # The order of two ends is irrelevant 
    A1 = PolymerBlock(:A1, sA, 0.2, FreeEnd(:A1), eb1)
    B = PolymerBlock(:B, sB, 0.3, eb1, eb2)
    A2 = PolymerBlock(:A2, sA, 0.3, eb2, eb3)
    C1 = PolymerBlock(:C1, sC, 0.1, eb3, FreeEnd(:C1))
    C2 = PolymerBlock(:C2, sC, 0.05, eb2, FreeEnd(:C2))
    D = PolymerBlock(:D, sD, 0.05, eb3, FreeEnd(:D))
    # The order of these blocks is also irrelevant
    return BlockCopolymer(:ABCACD, [A1, B, A2, C1, C2, D])
end

chainABCACD (generic function with 1 method)

chainABCACD()

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/8h0bl/src/scenes.jl:227

Exercises

1. Write a function that creates a AB3 star block copolymer

 

  [B1]
  |
[A] - [B2]
  |
  [B3]

1. Write a function that creates the following block copolymer

 

    [A2]         [B2]
       |         |
[A3] - [B1] - [A1] - [B3]
       |         |
    [A4]         [B4]

Polymer systems

Before adding BlockCopolymer or SmallMolecule objects to PolymerSystem, we have to wrap each of them into a Component struct. A Component object carries extra information such as its volume fraction and its reference length.

AB_system()  # χN and fA can be specified by the keyword arguments

PolymerSystem{Float64}(Component[Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.5 of specie A
  * PolymerBlock B with f=0.5 of specie B
, 1.0, 1.0)], BulkConfinement(), Flory-Huggins interaction parameters betwen species:
    (A, B) => 20.0
, 1.0)

ABC_system()

PolymerSystem{Float64}(Component[Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer ABC with 3 blocks:
  * PolymerBlock A with f=0.3 of specie A
  * PolymerBlock B with f=0.4 of specie B
  * PolymerBlock C with f=0.29999999999999993 of specie C
, 1.0, 1.0)], BulkConfinement(), Flory-Huggins interaction parameters betwen species:
    (A, C) => 40.0
    (B, C) => 40.0
    (A, B) => 40.0
, 1.0)

For any block copolymer melts, we can use the PolymerSystem constructor:

let
    chain = chainABCACD()  # Create a polymer chain
    polymer = Component(chain)  # Wrap the chain into a component
    # Create an interaction map between different species
    χNmap = Dict([:A, :B]=>10.0,
                 [:A, :C]=>20.0,
                 [:A, :D]=>30.0,
                 [:B, :C]=>40.0,
                 [:B, :D]=>50.0,
                 [:C, :D]=>60.0
    )
    PolymerSystem([polymer], χNmap)  # Create a polymer system
end

PolymerSystem{Float64}(Component[Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer ABCACD with 6 blocks:
  * PolymerBlock A1 with f=0.2 of specie A
  * PolymerBlock B with f=0.3 of specie B
  * PolymerBlock A2 with f=0.3 of specie A
  * PolymerBlock C1 with f=0.1 of specie C
  * PolymerBlock C2 with f=0.05 of specie C
  * PolymerBlock D with f=0.05 of specie D
, 1.0, 1.0)], BulkConfinement(), Flory-Huggins interaction parameters betwen species:
    (D, B) => 50.0
    (A, C) => 20.0
    (B, C) => 40.0
    (A, B) => 10.0
    (A, D) => 30.0
    (D, C) => 60.0
, 1.0)

A/B binary blend

First, we create two homopolymers chainA and chainB via the convenient function homopolymer_chain of Polymer.jl.

chainA = homopolymer_chain(label=:A)
chainB = homopolymer_chain(label=:B)
nothing

Then we wrap them into two Component objects. Below, we make them equal length and equal volume fraction.

polymerA = Component(chainA, 1.0, 0.5)  # αA = 1.0, ϕA = 0.5
polymerB = Component(chainB, 1.0, 0.5)  # αB = 1.0, ϕB = 0.5

Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer B with 1 blocks:
  * PolymerBlock B with f=1.0 of specie B
, 1.0, 0.5)

Next, we will specify the interactions via a Dict object:

χNmap = Dict([:A, :B]=>20.0)  # requires Polymer.jl v0.8.5

Dict{Vector{Symbol}, Float64} with 1 entry:
  [:A, :B] => 20.0

Finally, we can construct a PolymerSystem object for a A/B binary blend:

PolymerSystem([polymerA, polymerB], χNmap)

PolymerSystem{Float64}(Component[Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer A with 1 blocks:
  * PolymerBlock A with f=1.0 of specie A
, 1.0, 0.5), Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer B with 1 blocks:
  * PolymerBlock B with f=1.0 of specie B
, 1.0, 0.5)], BulkConfinement(), Flory-Huggins interaction parameters betwen species:
    (A, B) => 20.0
, 1.0)

Note that a convenient function A_B_system has been provide by Polymer.jl to do exactly the above things.

A_B = A_B_system()

PolymerSystem{Float64}(Component[Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer A with 1 blocks:
  * PolymerBlock A with f=1.0 of specie A
, 1.0, 0.5), Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer B with 1 blocks:
  * PolymerBlock B with f=1.0 of specie B
, 1.0, 0.5)], BulkConfinement(), Flory-Huggins interaction parameters betwen species:
    (A, B) => 20.0
, 1.0)

AB/A binary blend

This can be done via the function AB_A_system provided by the package Polymer.jl.

AB_A_system accepts keyword arguments only. Its signature is

AB_A_system(; χN=20.0, ϕAB=0.5, fA=0.5, α=0.5)

The following is assumed:

• αAB = 1.0
• αA = α
• fA is the length fraction of A block within the AB block copolymer.

Below we create a binary blend which is the same as that in Fig. 3 of Ref. Mester, Z.; Lynd, N. A.; Fredrickson, G. H. Numerical Self-Consistent Field Theory of Multicomponent Polymer Blends in the Gibbs Ensemble. Soft Matter 2013, 9 (47), 11288.

AB_A = AB_A_system(; χN=10.0, fA=0.45, α=1.0)  # function from Polymer.jl

PolymerSystem{Float64}(Component[Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.45 of specie A
  * PolymerBlock B with f=0.55 of specie B
, 1.0, 0.5), Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer hA with 1 blocks:
  * PolymerBlock hA with f=1.0 of specie A
, 1.0, 0.5)], BulkConfinement(), Flory-Huggins interaction parameters betwen species:
    (A, B) => 10.0
, 1.0)

AB/A/B ternary blend

Based on Polymer.jl, arbitrary polymer blends including polymer solutions can be created. Below is an example on how to create an AB/A/B ternary blend.

function AB_A_B_system(; χN=20.0, ϕAB=0.5, fA=0.5, ϕA=0.25, αA=0.5, αB=0.5)
    AB = Component(diblock_chain(fA=fA), 1.0, ϕAB)
    A = Component(homopolymer_chain(label=:hA, segment=KuhnSegment(:A)), αA, ϕA)
    B = Component(homopolymer_chain(label=:hB, segment=KuhnSegment(:A)), αB, 1-ϕAB-ϕA)
    return PolymerSystem([A, B, AB], Dict([:A, :B]=>χN))
end

AB_A_B_system (generic function with 1 method)

Note that we put AB as the last component in the PolymerSystem object, which follows literature convention so that chemical related properties are reduced quantities using the last component as a reference.

AB_A_B = AB_A_B_system()

PolymerSystem{Float64}(Component[Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer hA with 1 blocks:
  * PolymerBlock hA with f=1.0 of specie A
, 0.5, 0.25), Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer hB with 1 blocks:
  * PolymerBlock hB with f=1.0 of specie A
, 0.5, 0.25), Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.5 of specie A
  * PolymerBlock B with f=0.5 of specie B
, 1.0, 0.5)], BulkConfinement(), Flory-Huggins interaction parameters betwen species:
    (A, B) => 20.0
, 1.0)

Exercises

1. Create a block copolymer melt system contains a miktoarm chain AB3. The interaction between A and B is 30.0.

2. Create a AB3/A binary blend. The volume fraction of homopolymer A is 0.1. The reference length is the length of AB3, and the length of homopolymer A is 0.2. The interaction between A and B is 30.0.

Create Simulation Cell

A polymer system has to be put in a simulation cell, which carries the information about the cell size, shape and symmetry. A simulation cell is a Scattering.BravaisLattice object.

Unit cells and crystal systems

Scattering.jl supports all crystal systems in 1D, 2D and 3D. A full list can be viewed by running:

values(Scattering.Crystal.ALL_CRYSTAL_SYSTEMS)

ValueIterator for a Dict{Tuple{Symbol, Symbol, Vararg{Symbol}}, DataType} with 14 entries. Values:
  Oblique
  Hexagonal
  Rectangular
  HexRect
  Line
  Hexagonal2D
  Triclinic
  Orthorhombic
  HexOrthorhombic
  Square
  Cubic
  Tetragonal
  Trigonal
  Monoclinic

Unit cell can be created by the UnitCell constructor function. One way is to specify crystal system and corresponding crystalline parameters.

UnitCell(Cubic(), 4.0)  # create a cubic cell

UnitCell
  * Crystal system: Cubic
  * Edges: [4.0, 4.0, 4.0]
  * Angles: [π/2, π/2, π/2]

Another way is to specify edge lengths and angles. If angles are ignored, they are assumed to be $\pi/2$. The crystal system is inferred from the edge lengths and angles.

let
    uc = UnitCell((2.0, 3.0, 4.0))
    crystalsystem(uc)
end

Orthorhombic crystal system with free lattice parameters [a, b, c]

Space Group

• For 1D crystals, there are two space groups.
• For 2D crystals, there are 17 space groups.
• For 3D crystals, there are 230 space groups.

These space groups can be described by the SpaceGroup struct defined in Crystalline.jl package. They can be created by specifying the space group index and the number of space dimension. For the name of each space group, please check out the list here.

spacegroup(70, 3)  # O70 (Fddd) space group

SpaceGroup{3} ⋕70 (Fddd) with 32 operations:
 1
 {2₀₀₁|¾,¾,0}
 {2₀₁₀|¾,0,¾}
 {2₁₀₀|0,¾,¾}
 -1
 {m₀₀₁|¼,¼,0}
 {m₀₁₀|¼,0,¼}
 {m₁₀₀|0,¼,¼}
 {1|0,½,½}
 {2₀₀₁|¾,¼,½}
 {2₀₁₀|¾,½,¼}
 {2₁₀₀|0,¼,¼}
 {-1|0,½,½}
 {m₀₀₁|¼,¾,½}
 {m₀₁₀|¼,½,¾}
 {m₁₀₀|0,¾,¾}
 {1|½,0,½}
 {2₀₀₁|¼,¾,½}
 {2₀₁₀|¼,0,¼}
 {2₁₀₀|½,¾,¼}
 {-1|½,0,½}
 {m₀₀₁|¾,¼,½}
 {m₀₁₀|¾,0,¾}
 {m₁₀₀|½,¼,¾}
 {1|½,½,0}
 {2₀₀₁|¼,¼,0}
 {2₀₁₀|¼,½,¾}
 {2₁₀₀|½,¼,¾}
 {-1|½,½,0}
 {m₀₀₁|¾,¾,0}
 {m₀₁₀|¾,½,¼}
 {m₁₀₀|½,¾,¼}

Lattice

A simulation cell can be created by specifying the unit cell only. The crystal system is inferred from the unit cell. The space group is chosen to be the first space group of that crystal system.

let
    uc = UnitCell((2.0, 3.0, 4.0))
    BravaisLattice(uc)
end

BravaisLattice
  * Centering: P
  * Space group: #16 (P222)
  * Crystal system: Orthorhombic
  * Unit cell: [2.0, 3.0, 4.0] [π/2, π/2, π/2]
  * Free lattice parameters: [a, b, c]

Alternatively, the space group can be specified explicitly as

let
    uc = UnitCell(Orthorhombic(), 2.0, 3.0, 4.0)
    BravaisLattice(uc, 40)  # No. 40 space group (3D) because uc is 3D.
end

BravaisLattice
  * Centering: A
  * Space group: #40 (Ama2)
  * Crystal system: Orthorhombic
  * Unit cell: [2.0, 3.0, 4.0] [π/2, π/2, π/2]
  * Free lattice parameters: [a, b, c]

An advanced usage of crystal system is to specify which unit cell parameters are fixed. The b edge length of following orthorhombic cell is fixed.

Note that UnitCell object does not explicitly carry any crystal system information. Thus control parameters cannot be set in UnitCell.

let
    uc = UnitCell((2.0, 3.0, 4.0))
    a = Scattering.Crystal.a
    c = Scattering.Crystal.c
    lat = BravaisLattice(uc, 40; freevars=[a, c])
end

BravaisLattice
  * Centering: A
  * Space group: #40 (Ama2)
  * Crystal system: Orthorhombic
  * Unit cell: [2.0, 3.0, 4.0] [π/2, π/2, π/2]
  * Free lattice parameters: [a, c]

Exercises

1. Create a cubic cell with space group 230 and side length 5.0.

2. Create a 2D hexagonal cell with default space group and side length 3.0.

Construct SCFT Models

An SCFT model holds all information needed for performing SCFT calculations. It typically includes:

• A PolymerSystem object: describing the polymer system
• A BravaisLattice object: describing the simulation cell
• An array of density fields
• An array of auxiliary fields
• An array of propagators
• An SCFT udpater
• An array of MDE solvers

Currently, Polyorder only supports non-cyclic chain architechures.

Constructor

The model can be created by the constructor of NoncyclicChainSCFT struct as follows,

let
    uc = UnitCell(4.0)  # Create a unit cell
    lat = BravaisLattice(uc)  # Create a simulation cell
    system = AB_system()  # Create a diblock copolymer system
    ds = 0.01  # global contour step size for all blocks
    model = NoncyclicChainSCFT(system, lat, ds)  # Create the SCFT model
end

Noncyclic Chain SCFT model:
* Free energy: NaN
* Residual: 0.0
* Stress norm: NaN
* Model type: SimpleFieldModel
* Model compressiblity: Incompressible
-----
PolymerSystem (AB) contains 1 components:

Component AB with ϕ=1.0 and α=1.0 contains BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.5 of specie A
  * PolymerBlock B with f=0.5 of specie B


with Flory-Huggins interaction parameters betwen species:
    (A, B) => 20.0
-----
Simulation Cell: BravaisLattice
  * Centering: p
  * Space group: #1 (p1)
  * Crystal system: Line
  * Unit cell: [4.0] [0]
  * Free lattice parameters: [a]
-----
* Spatial resolution: (27,)
* Contour steps: 0.01
* MDE solvers: Dict[Dict{Any, Any}((2 => 3) => OSF, (2 => 1) => OSF, (1 => 2) => OSF, (3 => 2) => OSF)]
-----
SCFT updater: PicardMann iteration with α=0.2.
-----
* First density field: 0.0 [0.0, 0.0]
* First auxiliary field: -0.002156 [-0.4471, 0.4528]

Configurations

Note that density fields, auxiliary fields, and propagators are created automatically according to the given simulation cell. Other configurable parameters, such as spatial resolution is given by the default Polyorder.Config() instance. You can check out the default settings by create a default Polyorder.Config instance:

JP.Config()

Polyorder.Config(PolymerSystemConfig(nothing, SpecieConfig[SpecieConfig(:A, :Segment, 1.0, nothing), SpecieConfig(:B, :Segment, 1.0, nothing)], Vector{Any}[[:A, :B, 20.0]], ComponentConfig[ComponentConfig(:BCP, :AB, 1.0, 1.0, BlockConfig[BlockConfig(:A, :A, 0.5, [:AB]), BlockConfig(:B, :B, 0.5, [:AB])])], 1.0), LatticeConfig(CrystalSystemConfig(:Line, Symbol[]), :P, UnitCellConfig(1, 1.0, nothing, nothing, nothing, nothing, nothing), SpaceGroupConfig(1, 1)), IOConfig(1, false, false, ".", "summary", "trace", "fields", "densities", "config", :HDF5, true, true, true, true, true, 100, 1000, 1000, 6, 0), SCFTConfig(:SD, false, :simple, false, [0.1, 0.5], 0.02, 0.9, 1, 100, 2000, 100, 1, 0.5, :vecnormInf, :vecnormInf, false, 0.15, false, :Residual, 1.0e-5, 0.001, 1.0, 100, 1.4901161193847656e-8, 100, 1.0, 100, 1.4901161193847656e-8, 10, 10, false, false, false, false, false), MDEConfig(:RQM4, Symbol[], :Fixf, [0.01], Int64[]), CellOptConfig(Symbol("Optim.Brent"), Symbol("Optim.NelderMead"), 0, true, 2.0, 0.0001, 1.0e-5, false, false, 50))

The default MDE solver is RQM4, and the default SCFT updater is SD. You can change them by passing them through keywords mde and updater, respectively:

let
    uc = UnitCell(4.0)  # Create a unit cell
    lat = BravaisLattice(uc)  # Create a simulation cell
    system = AB_system()  # Create a diblock copolymer system
    ds = 0.01  # global contour step size for all blocks
    # Create an SCFT model with the `OSF` MDE solver and the `SIS` updater.
    model = NoncyclicChainSCFT(system, lat, ds; mde=OSF, updater=SIS(0.1))
end

Noncyclic Chain SCFT model:
* Free energy: NaN
* Residual: 0.0
* Stress norm: NaN
* Model type: SimpleFieldModel
* Model compressiblity: Incompressible
-----
PolymerSystem (AB) contains 1 components:

Component AB with ϕ=1.0 and α=1.0 contains BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.5 of specie A
  * PolymerBlock B with f=0.5 of specie B


with Flory-Huggins interaction parameters betwen species:
    (A, B) => 20.0
-----
Simulation Cell: BravaisLattice
  * Centering: p
  * Space group: #1 (p1)
  * Crystal system: Line
  * Unit cell: [4.0] [0]
  * Free lattice parameters: [a]
-----
* Spatial resolution: (27,)
* Contour steps: 0.01
* MDE solvers: Dict[Dict{Any, Any}((2 => 3) => OSF, (2 => 1) => OSF, (1 => 2) => OSF, (3 => 2) => OSF)]
-----
SCFT updater: 1S semi-implicit method (SIS-1) with α=0.1.
-----
* First density field: 0.0 [0.0, 0.0]
* First auxiliary field: 0.0001273 [-0.4885, 0.5]

Or equivalently, create a Polymer.Config instance with configured mde and scft fields:

let
    mde = MDEConfig(algo=:OSF)
    scft = SCFTConfig(algo=:Euler)
    
    config = JP.Config(mde=mde, scft=scft)
    
    # Create an SCFT model with the `OSF` MDE solver and the `Euler` updater.
    model = NoncyclicChainSCFT(config)

    model.updater isa Euler
end

true

It is useful to serialize SCFT models as a Polyorder.Config. Polyorder provides a function to_config to do it. And save_config to save the Polyorder.Config instance as an YAML file.

let
    uc = UnitCell(4.0)  # Create a unit cell
    lat = BravaisLattice(uc)  # Create a simulation cell
    system = AB_system()  # Create a diblock copolymer system
    ds = 0.01  # global contour step size for all blocks
    # Create an SCFT model with the `OSF` MDE solver and the `SIS` updater.
    model = NoncyclicChainSCFT(system, lat, ds; mde=OSF, updater=SIS(0.1))
    
    config = to_config(model, JP.Config())  # serialize a model as a Config object
    save_config("config.yml", config)  # save the Config object as an YAML file
end

Exercises

1. Construct a SCFT model for a miktoarm block copolymer AB3 (fA=0.4, fB1=fB2=fB3=0.2), in a cubic unit cell with side length 4.0 Rg and space group No.230.

2. Construct a SCFT model for a AB3/A binary blend (fA=0.4, fB1=fB2=fB3=0.2, choose AB3 as reference length, the length of homopolymer A is 0.5), in an orthorhombic unit cell with side length 3.0, 4.0, 5.0, respectively. Use RQM4 as MDE solver, ETD as updater. Set contour step for all A blocks (A in AB3 and A in hompolymer) to be 0.02, contour step for all B blocks to be 0.01.

3. Explore the source file Polyorder.jl/src/config.jl, and construct a configuration instance that set symmetrize to false, minimum and maximum number of iterations to be 50 and 500, respectively in SCFTConfig.

Disordered (DIS) phase

The analytical solution at the mean-field level is available, which is the same as the Flory-Huggins lattice theory,

\[\tilde{F} = \frac{F}{CV} = \sum_{X,Y}\chi_{XY}N \bar\phi_X \bar\phi_Y + \sum_i \frac{\phi_i}{\alpha_i}\left( \ln\frac{C\phi_i}{\alpha_i} - 1 \right)\]

where $X, Y = \lbrace A, B, C, \dots \rbrace$ and $Y > X$ are the specie types in the polymer system, and $i = 1, 2, 3, \dots$ are the components in the polymer system. The averaging density for the specie X, $\bar\phi_X$, is

\[\bar\phi_X = \sum_i\sum_j \delta_{XS_{ij}} f_{ij} \phi_i\]

where $S_{ij}$ is the specie type in the $j$-th block of the $i$-th molecule, and $\delta_{XY}$ is the Kronecker delta function. Now we can compute its first-order partial derivative

\[\frac{\partial\bar\phi_X}{\partial\phi_i} = \sum_j \delta_{XS_{ij}} f_{ij}\]

Therefore,

\[\gamma_i = \frac{\partial\tilde{F}}{\partial\phi_i} = \sum_{X,Y}\sum_j ( \delta_{XS_{ij}}\bar\phi_Y + \delta_{YS_{ij}}\bar\phi_X)\chi_{XY}N f_{ij} + \frac{1}{\alpha_i}\ln\frac{C\phi_i}{\alpha_i}\]

and the chemical potential can be readily computed

\[\tilde\mu_i = \gamma_i - \gamma_{n_c}\]

for $i=1, 2, \dots, n_c$ with $n_c$ the number of components in the system.

Polyorder provides two ways to compute the disordered (homogeneous) phase: the analytical way and the numerical way.

Free energy

Create a NoncyclicChainSCFT instance.

begin
    uc_dis = UnitCell(1.0)  # 1D unit cell with length 1.0 Rg
    lattice_dis = BravaisLattice(uc_dis)
    # Space resolution `spacing` Δx is set to be equal to the length of the unit cell to ensure only one grid point is generated.
    nccscft_dis = NoncyclicChainSCFT(AB_A, lattice_dis, 0.01; spacing=1.0, mde=OSF)
end

Noncyclic Chain SCFT model:
* Free energy: Inf
* Residual: 0.0
* Stress norm: NaN
* Model type: SimpleFieldModel
* Model compressiblity: Incompressible
-----
PolymerSystem (AB + hA) contains 2 components:

Component AB with ϕ=0.5 and α=1.0 contains BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.45 of specie A
  * PolymerBlock B with f=0.55 of specie B

Component hA with ϕ=0.5 and α=1.0 contains BlockCopolymer hA with 1 blocks:
  * PolymerBlock hA with f=1.0 of specie A


with Flory-Huggins interaction parameters betwen species:
    (A, B) => 10.0
-----
Simulation Cell: BravaisLattice
  * Centering: p
  * Space group: #1 (p1)
  * Crystal system: Line
  * Unit cell: [1.0] [0]
  * Free lattice parameters: [a]
-----
* Spatial resolution: (1,)
* Contour steps: 0.01
* MDE solvers: Dict[Dict{Any, Any}((2 => 3) => OSF, (2 => 1) => OSF, (1 => 2) => OSF, (3 => 2) => OSF), Dict{Any, Any}((2 => 1) => OSF, (1 => 2) => OSF)]
-----
SCFT updater: PicardMann iteration with α=0.2.
-----
* First density field: 0.0 [0.0, 0.0]
* First auxiliary field: 0.1164 [0.1164, 0.1164]

It is not recommended to use solve! method to compute the DIS phase. But it is good for testing whether the implementation of a SCFT algorithm is correct. Note that, to ensure that the DIS phase is computed, the number of grid points must be set to 1 to suppress all possible inhomogeneous phases.

The numerical way only serves as a test tool to verify the correctness of a SCFT implementation.

JP.solve!(nccscft_dis)  # solve the SCFT model numerically

[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: Inf
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: NaN
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [1.0]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: final loss: 9.170760709498182e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 0.3006028194
[ Info: final loss: 9.17e-6
[ Info: final stress norm: 0.0
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [1.0]
└   * Angles: [0]
[ Info: iterations: 57
[ Info: time per iteration: 30 milliseconds, 167 microseconds, 769 nanoseconds

[ Info: Run time: 00:00:01.719562833
[ Info: =======================================

Polyorder.Successful()

We can simply obtain the free energy (which includes the contribution from the ideal gas terms) by calling the F function on the solved NoncyclicChainSCFT object.

JP.F(nccscft_dis)  # Evaluate the free energy numerically

0.3006028194400485

Compute the enthalpy part of the Flory-Huggins free energy.

JP.enthalpy(nccscft_dis)  # Evaluate the enthalpy numerically

1.9937500000000032

JP.enthalpy_FH(nccscft_dis)  # Evaluate the enthalpy analytically using FH theory.

1.9937500000000001

Compute the entropy part (i.e., the ideal gas terms) of the Flory-Huggins free energy,

\[\tilde{S} = -\tilde{F}_{ig} = \sum_i \frac{\phi_i}{\alpha_i}\left( \ln\frac{C\phi_i}{\alpha_i} - 1 \right)\]

JP.entropy(nccscft_dis)  # Evaluate the entropy numerically

1.6931471805599545

JP.entropy_ig(nccscft_dis)  # Evaluate the entropy analytically using FH theory

1.6931471805599454

JP.F_ig(nccscft_dis)  # shoud = -S

-1.6931471805599454

We can combine these two parts to obtain the full Flory-Huggins free energy: F = H - TS.

JP.enthalpy_FH(nccscft_dis) - JP.entropy_ig(nccscft_dis) # should be equal to JP.F(nccscft_dis)

0.30060281944005474

A convenient function is provided in the Polyorder package to compute the Flory-Huggins free energy, which is F_DIS or F_FH. Analytical expressions are utilized to carry out the computation.

We can verify the accuracy of the DIS free energy calculated from SCFT equations:

JP.F_DIS(nccscft_dis) - JP.F(nccscft_dis)

6.217248937900877e-15

It is not necessary to create an NoncyclicChainSCFT object to compute the disordered phases. Using PolymerSystem directly is the recommended way to compute the DIS free energy:

JP.F_DIS(AB_A)

0.30060281944005474

JP.F_DIS(AB_A_B)

1.210279229160082

JP.F_DIS(A_B)

3.3068528194400546

Chemical potential

Compute chemical potentials from the solved NoncyclicChainSCFT object.

JP.μ̃s(nccscft_dis)  # Evaluate the chemical potential numerically

1-element Vector{Float64}:
 2.474989912163239

Compute chemical potentials based on Flory Huggins lattice theory.

JP.μ̃s_DIS(nccscft_dis)  # Evaluate the chemical potential analytically

1-element Vector{Float64}:
 2.4750000000000005

Verify the accuracy of the SCFT computation.

JP.μ̃s(nccscft_dis) .- JP.μ̃s_DIS(nccscft_dis)

1-element Vector{Float64}:
 -1.0087836761485391e-5

Similarly, using PolymerSystem directly is the recommended way to compute the DIS chemical potential.

# the modified chemical potentials for each component except the last one.
# the order of component is the same as the array `PolymerSystem.components`.
JP.μ̃s_DIS(AB_A)

1-element Vector{Float64}:
 2.4750000000000005

AB_A.components

2-element Vector{Component}:
 Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.45 of specie A
  * PolymerBlock B with f=0.55 of specie B
, 1.0, 0.5)
 Component{BlockCopolymer{PolymerBlock{Float64}}, Float64}(BlockCopolymer hA with 1 blocks:
  * PolymerBlock hA with f=1.0 of specie A
, 1.0, 0.5)

Thus the chemical potential is for the AB diblock copolymer in the above case.

JP.γs_DIS(AB_A_B)  # the partial derivatives for each component.

3-element Vector{Float64}:
 3.613705638880109
 3.613705638880109
 9.306852819440055

JP.μ̃s_DIS(AB_A_B)  # the modified chemical potentials for each component (AB and A) except the last one (B).

2-element Vector{Float64}:
 -5.693147180559945
 -5.693147180559945

JP.μs_DIS(AB_A_B)  # the chemical potentials for each component.

3-element Vector{Float64}:
 -0.8181471805599454
 -0.8181471805599454
  4.056852819440055

Batch evaluations

Plot the free energy as a function of the volume fraction of the AB component for the AB/A binary blend.

begin
    system = AB_A_system(; χN=10.0, fA=0.45, α=1.0)
    ϕA_control = ϕControlParameter(:hA, system)  # a handle for modifying PolymerSystem
end

ϕControlParameter{Polymer.var"#30#32"}(2, ϕType{Float64}("volume fraction of a type of chain in a polymer system.", "ϕ", "phi", L"$\phi$", Float64, 0.0, 1.0, Float64[], Float64[]), Polymer.var"#30#32"())

begin
    ϕAs_dis = 0.001:0.001:0.999  # ϕA
    Fs_dis = similar(ϕAs_dis)
    μs_dis = similar(ϕAs_dis)
    for i in eachindex(ϕAs_dis)
        Polymer.update!(system, ϕAs_dis[i], ϕA_control)
        Fs_dis[i], μs_dis[i] = JP.F_DIS(system), JP.μ̃_DIS(system, 1)
    end
end

let
    f = Figure()
    ax = Axis(f[1, 1],
        xlabel = L"$\phi_{hA}$",
        ylabel = L"$F$ and $u$",
    )
    lines!(ax, ϕAs_dis, Fs_dis, label=L"$F$")
    f
end

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/8h0bl/src/scenes.jl:227

A straight line is subtracted from the actual free energy curve to visualize the double-well structure more easily.

let
    line = get_line(ϕAs_dis[1], Fs_dis[1], ϕAs_dis[end], Fs_dis[end])
    F_sub = subtract_line(line, ϕAs_dis, Fs_dis)
    f = Figure()
    ax = Axis(f[1, 1],
        xlabel = L"$\phi_{hA}$",
        ylabel = L"$F$",
    )
    lines!(ax, ϕAs_dis, F_sub, label="DIS")
    f
end

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/8h0bl/src/scenes.jl:227

Exercises

1. Draw DIS phase free energy curves of the A/B binary blend for χN={1.5, 1.8, 1.99, 2.01, 2.1, 3.0}, and other values which are helpful to understand the trend of the change of these free energy curves.

Ordered Phases

Block copolymers undergo a microphase separation when the interaction between species. Ordered phases emerges during the microphase separation. The ordered structures can be computed by SCFT numerically. Polyorder provides solve! and cell_solve! for fixed cell calculations and cell optimization, respectively.

Here, we use the most simple block copolymer system, AB diblock copolymer melt to demonstrate the usage of Polyorder.

LAM

By restricting in the 1D unit cell, the only possible phases are LAM and DIS.

Sometimes, when $\chi N$ is significantly close to the critical point or the cell size significantly deviates from the equilibrium value, the DIS is computed instead of LAM. Therefore, for a rigorous approach, we must check the density distribution of the final solution to make sure the LAM is actually computed. Or, simply compared its free energy with the DIS phase, if they are close within a threshold (say 1e-8), we can conclude that the DIS instead of LAM is computed.

If DIS has been computed when the target phase is an ordered phase, we can

• make sure that $\chi N$ is indeed above the binodal line.
• change the cell size
• provide an initial density distribution with a larger variation.

scft_AB_lam = let
    ds = 0.01
    uc = UnitCell(3.0)
    lat = BravaisLattice(uc)
    system = AB_system(χN=15.0, fA=0.36)
    NoncyclicChainSCFT(system, lat, ds)
end

Noncyclic Chain SCFT model:
* Free energy: 723.4313150479
* Residual: 0.0
* Stress norm: NaN
* Model type: SimpleFieldModel
* Model compressiblity: Incompressible
-----
PolymerSystem (AB) contains 1 components:

Component AB with ϕ=1.0 and α=1.0 contains BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.36 of specie A
  * PolymerBlock B with f=0.64 of specie B


with Flory-Huggins interaction parameters betwen species:
    (A, B) => 15.0
-----
Simulation Cell: BravaisLattice
  * Centering: p
  * Space group: #1 (p1)
  * Crystal system: Line
  * Unit cell: [3.0] [0]
  * Free lattice parameters: [a]
-----
* Spatial resolution: (20,)
* Contour steps: 0.01
* MDE solvers: Dict[Dict{Any, Any}((2 => 3) => OSF, (2 => 1) => OSF, (1 => 2) => OSF, (3 => 2) => OSF)]
-----
SCFT updater: PicardMann iteration with α=0.2.
-----
* First density field: 0.0 [0.0, 0.0]
* First auxiliary field: -0.035 [-0.3845, 0.4997]

Fixed cell calculation

The SCFT equations are solved for the current set of parameters using solve!. The model is updated in place. Note that the cell is not optimized.

JP.solve!(scft_AB_lam)

[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 723.4313150479
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: NaN
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.0]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >

[ Info: number: 100                                                                     
[ Info: F: 2.3537691835                                                                 
[ Info: resediual norm: 0.364                                                           
[ Info: number: 200                                                                     
[ Info: F: 2.3619065495                                                                 
[ Info: resediual norm: 0.00167                                                         
[ Info: number: 300                                                                     
[ Info: F: 2.3619038877                                                                 
[ Info: resediual norm: 0.000186                                                        
[ Info: number: 400                                                                     
[ Info: F: 2.3619039908                                                                 
[ Info: resediual norm: 2.83e-5                                                         
[ Info: final loss: 9.961214716192046e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3619039968
[ Info: final loss: 9.96e-6
[ Info: final stress norm: 0.09399202959385276
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.0]
└   * Angles: [0]
[ Info: iterations: 459
[ Info: time per iteration: 174 microseconds, 36 nanoseconds
[ Info: Run time: 00:00:00.079882625
[ Info: =======================================

Polyorder.Successful()

We can plot out the density fields to verify that the LAM phase is obtained.

plot_density(scft_AB_lam)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/8h0bl/src/scenes.jl:227

The free energy and other properties can be achieved from the solved model through various methods, such as F, μ, etc.

JP.F(scft_AB_lam)

2.361903996752715

JP.lattice(scft_AB_lam)

BravaisLattice
  * Centering: p
  * Space group: #1 (p1)
  * Crystal system: Line
  * Unit cell: [3.0] [0]
  * Free lattice parameters: [a]

Cell optimization

The cell should be optimized to have zero stress, i.e. the phase separation structure is stress-free. Polyorder provides three ways to optimize the cell with a unified interface function cell_solve!:

• Variable cell optimization: this is the most recommended way if you know a good initial guess. It optimizes the auxiliary fields and cell shape/size simultaneously. Therefore, it enjoys a fast convergence.
• Stress-guided cell optimization: using optimization methods to find the zero-stress unit cell with the aid of the gradient (stress) with respect to cell size and shape. This is a compromising approach to obtain the stress-free solution of the SCFT equations.
• Gradient-free cell optimization: it optimizes the unit cell sizes using gradient-free optimization methods. Note that only orthogonal unit cells are supported.

Variable cell optimization

To invoke the variable cell optimization mode via cell_solve! method, you need to either specify a VariableCell updater for the SCFT model or explicitly pass a VariableCell updater as an input argument to the cell_solve! method.

Note that the initial guess of the cell size is critical to find the equilibrium cell size. Here, by trial and error, we find the equilibrium cell size is about 7.0 $R_g$.

scft_vc = JP.clone(scft_AB_lam)

Noncyclic Chain SCFT model:
* Free energy: 721.1879402971
* Residual: 0.0
* Stress norm: NaN
* Model type: SimpleFieldModel
* Model compressiblity: Incompressible
-----
PolymerSystem (AB) contains 1 components:

Component AB with ϕ=1.0 and α=1.0 contains BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.36 of specie A
  * PolymerBlock B with f=0.64 of specie B


with Flory-Huggins interaction parameters betwen species:
    (A, B) => 15.0
-----
Simulation Cell: BravaisLattice
  * Centering: p
  * Space group: #1 (p1)
  * Crystal system: Line
  * Unit cell: [3.0] [0]
  * Free lattice parameters: [a]
-----
* Spatial resolution: (20,)
* Contour steps: 0.01
* MDE solvers: Dict[Dict{Any, Any}((2 => 3) => OSF, (2 => 1) => OSF, (1 => 2) => OSF, (3 => 2) => OSF)]
-----
SCFT updater: PicardMann iteration with α=0.2.
-----
* First density field: 0.0 [0.0, 0.0]
* First auxiliary field: 2.115 [-3.97, 8.235]

JP.cell_solve!(scft_vc, VariableCell(BB(1.0), SIS(1.0)))

[ Info: ###### Cell Optimization Start ######
┌ Info: Algorithm: Variable cell method
│ * cell updater: Barzilai-Borwein (BB2) method with max step size =1.0.
│ * fields updater: 1S semi-implicit method (SIS-1) with α=1.0.
└ Run fields updater 1 times per cell iteration.
[ Info: Cell optimization starts ...
[ Info: ******* SCFT Simulation Start *******
┌ Info: Algorithm: Variable cell method
│ * cell updater: Barzilai-Borwein (BB2) method with max step size =1.0.
│ * fields updater: 1S semi-implicit method (SIS-1) with α=1.0.
└ Run fields updater 1 times per cell iteration.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (stress norm)
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 721.1879402971
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: NaN
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.0]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: number: 100                                                                     

[ Info: F: 2.3352551579                                                                 
[ Info: resediual norm: 0.000136                                                        
[ Info: stress norm: 4.87e-7                                                            
[ Info: cell: [3.56891]                                                                 
[ Info: final loss: 9.638934645561648e-6
[ Info: Stop triggered by a `ThresholdStress` control. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.33525542
[ Info: final loss: 9.64e-6
[ Info: final stress norm: 1.43e-6
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.568909162030663]
└   * Angles: [0]
[ Info: iterations: 161
[ Info: time per iteration: 8 milliseconds, 454 microseconds, 567 nanoseconds
[ Info: Run time: 00:00:01.361185375
[ Info: =======================================
[ Info: Cell optimization finished.
[ Info: ------ Cell Optimization Summary ------

[ Info: Convergence: Polyorder.Successful()
┌ Info: Stress-free cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.568909162030663]
└   * Angles: [0]
[ Info: Final F: 2.335255420045568
[ Info: Final loss: 9.638934645561648e-6
[ Info: Final stress: 1.4250723279439873e-6
[ Info: Total solve! calls: 161
[ Info: Total SCFT iterations: 161
[ Info: Total time: 00:00:02

(Polyorder.Successful(), Noncyclic Chain SCFT model:
* Free energy: 2.33525542
* Residual: 9.639e-6
* Stress norm: 1.425e-6
* Model type: SimpleFieldModel
* Model compressiblity: Incompressible
-----
PolymerSystem (AB) contains 1 components:

Component AB with ϕ=1.0 and α=1.0 contains BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.36 of specie A
  * PolymerBlock B with f=0.64 of specie B


with Flory-Huggins interaction parameters betwen species:
    (A, B) => 15.0
-----
Simulation Cell: BravaisLattice
  * Centering: p
  * Space group: #1 (p1)
  * Crystal system: Line
  * Unit cell: [3.568909162030663] [0]
  * Free lattice parameters: [a]
-----
* Spatial resolution: (24,)
* Contour steps: 0.01
* MDE solvers: Dict[Dict{Any, Any}((2 => 3) => OSF, (2 => 1) => OSF, (1 => 2) => OSF, (3 => 2) => OSF)]
-----
SCFT updater: PicardMann iteration with α=0.2.
-----
* First density field: 0.36 [0.0617, 0.8146]
* First auxiliary field: 2.115 [-4.945, 8.1]
)

JP.F(scft_vc)

2.335255420045568

JP.lattice(scft_vc)

BravaisLattice
  * Centering: p
  * Space group: #1 (p1)
  * Crystal system: Line
  * Unit cell: [3.568909162030663] [0]
  * Free lattice parameters: [a]

Stress-guided cell optimization

If you call the cell_solve! method with only the SCFT model, i.e. no VariableCell updater is given as an input argument and the updater field of the SCFT model is not a VariableCell instance, then the stress-guided cell optimization is invoked.

scft_sgc = JP.clone(scft_AB_lam)

Noncyclic Chain SCFT model:
* Free energy: 743.4400719214
* Residual: 0.0
* Stress norm: NaN
* Model type: SimpleFieldModel
* Model compressiblity: Incompressible
-----
PolymerSystem (AB) contains 1 components:

Component AB with ϕ=1.0 and α=1.0 contains BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.36 of specie A
  * PolymerBlock B with f=0.64 of specie B


with Flory-Huggins interaction parameters betwen species:
    (A, B) => 15.0
-----
Simulation Cell: BravaisLattice
  * Centering: p
  * Space group: #1 (p1)
  * Crystal system: Line
  * Unit cell: [3.0] [0]
  * Free lattice parameters: [a]
-----
* Spatial resolution: (20,)
* Contour steps: 0.01
* MDE solvers: Dict[Dict{Any, Any}((2 => 3) => OSF, (2 => 1) => OSF, (1 => 2) => OSF, (3 => 2) => OSF)]
-----
SCFT updater: PicardMann iteration with α=0.2.
-----
* First density field: 0.0 [0.0, 0.0]
* First auxiliary field: 2.115 [-3.97, 8.235]

scft_sgc.updater  # updater is PicardMann (alias for SD) not VariableCell

PicardMann iteration with α=0.2.

JP.cell_solve!(scft_sgc)

[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 743.4400719214
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: NaN
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.0]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: final loss: 9.961214716192046e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3619039968
[ Info: final loss: 9.96e-6
[ Info: final stress norm: 0.09399202959385276
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.0]
└   * Angles: [0]
[ Info: iterations: 1
[ Info: time per iteration: 253 microseconds, 958 nanoseconds
[ Info: Run time: 00:00:00.000253958
[ Info: =======================================
[ Info: ###### Cell Optimization Start ######
[ Info: Algorithm: LBFGS
[ Info: Cell optimization starts ...

[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 2.3619039968
[ Info: initial residual norm: 9.96e-6
[ Info: initial stress norm: 0.094
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.0]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: final loss: 9.789049204123046e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3619039968
[ Info: final loss: 9.79e-6
[ Info: final stress norm: 0.093992029439695
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.0]
└   * Angles: [0]
[ Info: iterations: 1
[ Info: time per iteration: 280 microseconds, 666 nanoseconds
[ Info: Run time: 00:00:00.000280666
[ Info: =======================================
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 2.3619039968
[ Info: initial residual norm: 9.79e-6
[ Info: initial stress norm: 0.0857
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.093992029439695]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: number: 100                                                                     
[ Info: F: 2.3538079436                                                                 
[ Info: resediual norm: 0.000208                                                        
[ Info: number: 200                                                                     
[ Info: F: 2.3538082325                                                                 
[ Info: resediual norm: 1.12e-5                                                         
[ Info: final loss: 9.871820236595497e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3538082357
[ Info: final loss: 9.87e-6
[ Info: final stress norm: 0.07905519809631996
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.093992029439695]
└   * Angles: [0]
[ Info: iterations: 205
[ Info: time per iteration: 42 microseconds, 771 nanoseconds
[ Info: Run time: 00:00:00.008768083
[ Info: =======================================
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 720.9372703982
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: NaN
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.469960147198475]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: number: 100                                                                     
[ Info: F: 2.3360239927                                                                 
[ Info: resediual norm: 0.000591                                                        
[ Info: number: 200                                                                     
[ Info: F: 2.3360265282                                                                 
[ Info: resediual norm: 3.15e-5                                                         
[ Info: final loss: 9.780176487339531e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3360265871
[ Info: final loss: 9.78e-6
[ Info: final stress norm: 0.01612088457135104
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.469960147198475]
└   * Angles: [0]
[ Info: iterations: 248
[ Info: time per iteration: 51 microseconds, 122 nanoseconds
[ Info: Run time: 00:00:00.0126785
[ Info: =======================================
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 739.7765102753
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: NaN
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [5.349800735992375]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: number: 100                                                                     
[ Info: F: 2.4004884258                                                                 
[ Info: resediual norm: 0.187                                                           
[ Info: number: 200                                                                     
[ Info: F: 2.3993211334                                                                 
[ Info: resediual norm: 0.0122                                                          
[ Info: number: 300                                                                     
[ Info: F: 2.3993204423                                                                 
[ Info: resediual norm: 0.001                                                           
[ Info: number: 400                                                                     
[ Info: F: 2.399320443                                                                  
[ Info: resediual norm: 8.52e-5                                                         
[ Info: final loss: 9.945237756658898e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3993204447
[ Info: final loss: 9.95e-6
[ Info: final stress norm: 0.06627459198039426
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [5.349800735992375]
└   * Angles: [0]
[ Info: iterations: 487
[ Info: time per iteration: 78 microseconds, 296 nanoseconds
[ Info: Run time: 00:00:00.038130208
[ Info: =======================================
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: -537.4116114039
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: NaN
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [4.174900367996187]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >


[ Info: number: 100                                                                     
[ Info: F: 2.4574872526                                                                 
[ Info: resediual norm: 0.0547                                                          
[ Info: number: 200                                                                     
[ Info: F: 2.4563762019                                                                 
[ Info: resediual norm: 0.0254                                                          
[ Info: number: 300                                                                     
[ Info: F: 2.4546994387                                                                 
[ Info: resediual norm: 0.0997                                                          
[ Info: number: 400                                                                     
[ Info: F: 2.3586872282                                                                 
[ Info: resediual norm: 0.276                                                           
[ Info: number: 500                                                                     
[ Info: F: 2.3620967216                                                                 
[ Info: resediual norm: 0.000712                                                        
[ Info: number: 600                                                                     
[ Info: F: 2.3621019589                                                                 
[ Info: resediual norm: 3.41e-5                                                         
[ Info: final loss: 9.843467325687527e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3621020133
[ Info: final loss: 9.84e-6
[ Info: final stress norm: 0.08276625824473824
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [4.174900367996187]
└   * Angles: [0]
[ Info: iterations: 653
[ Info: time per iteration: 64 microseconds, 348 nanoseconds
[ Info: Run time: 00:00:00.042019458
[ Info: =======================================
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 720.3169816706
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: NaN
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.6247586544544896]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: number: 100                                                                     
[ Info: F: 2.3355178841                                                                 
[ Info: resediual norm: 0.000531                                                        
[ Info: number: 200                                                                     
[ Info: F: 2.3355153211                                                                 
[ Info: resediual norm: 3.74e-5                                                         
[ Info: final loss: 9.911253836314415e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3355152605
[ Info: final loss: 9.91e-6
[ Info: final stress norm: 0.008877477426241284
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.6247586544544896]
└   * Angles: [0]
[ Info: iterations: 258
[ Info: time per iteration: 51 microseconds, 391 nanoseconds
[ Info: Run time: 00:00:00.013258916
[ Info: =======================================
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 2.3355152605
[ Info: initial residual norm: 9.91e-6
[ Info: initial stress norm: 0.00929
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.570842960477296]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: number: 100                                                                     
[ Info: F: 2.3352564792                                                                 
[ Info: resediual norm: 6.67e-5                                                         
[ Info: final loss: 9.900442028420997e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3352561665
[ Info: final loss: 9.9e-6
[ Info: final stress norm: 0.0003117711104045113
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.570842960477296]
└   * Angles: [0]
[ Info: iterations: 167
[ Info: time per iteration: 53 microseconds, 728 nanoseconds
[ Info: Run time: 00:00:00.008972709
[ Info: =======================================
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 2.3352561665
[ Info: initial residual norm: 9.9e-6
[ Info: initial stress norm: 0.000376
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.355180184568523]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: number: 100                                                                     
[ Info: F: 2.338943139                                                                  
[ Info: resediual norm: 0.000318                                                        
[ Info: number: 200                                                                     
[ Info: F: 2.3389419407                                                                 
[ Info: resediual norm: 2.15e-5                                                         
[ Info: final loss: 9.890056571926209e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()





[ Info: final F: 2.3389419133
[ Info: final loss: 9.89e-6
[ Info: final stress norm: 0.035286904183002885
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.355180184568523]
└   * Angles: [0]
[ Info: iterations: 233
[ Info: time per iteration: 50 microseconds, 71 nanoseconds
[ Info: Run time: 00:00:00.011666666
[ Info: =======================================
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 2.3389419133
[ Info: initial residual norm: 9.89e-6
[ Info: initial stress norm: 0.0293
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.5689541985342235]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: number: 100                                                                     
[ Info: F: 2.3352539895                                                                 
[ Info: resediual norm: 0.000289                                                        
[ Info: number: 200                                                                     
[ Info: F: 2.3352554074                                                                 
[ Info: resediual norm: 1.61e-5                                                         
[ Info: final loss: 9.855689918403243e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3352554251
[ Info: final loss: 9.86e-6
[ Info: final stress norm: 8.649937644509866e-6
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.5689541985342235]
└   * Angles: [0]
[ Info: iterations: 220
[ Info: time per iteration: 50 microseconds, 736 nanoseconds
[ Info: Run time: 00:00:00.011161958
[ Info: =======================================
[ Info: Cell optimization finished.
[ Info: ------ Cell Optimization Summary ------
[ Info: Convergence: Unknown
┌ Info: Stress-free cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.5689541985342235]
└   * Angles: [0]
[ Info: Final F: 2.335255425070275
[ Info: Final loss: 9.855689918403243e-6
[ Info: Final stress: 8.649937644509866e-6
[ Info: Total solve! calls: 9
[ Info: Total SCFT iterations: 2472
[ Info: Total time: 00:00:01

( * Status: success

 * Candidate solution
    Final objective value:     2.335255e+00

 * Found with
    Algorithm:     L-BFGS

 * Convergence measures
    |x - x'|               = 5.58e-02 ≰ 1.0e-04
    |x - x'|/|x'|          = 1.56e-02 ≰ 0.0e+00
    |f(x) - f(x')|         = 2.60e-04 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 1.11e-04 ≰ 0.0e+00
    |g(x)|                 = 8.65e-06 ≤ 1.0e-05

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    2
    f(x) calls:    9
    ∇f(x) calls:   9
, Noncyclic Chain SCFT model:
* Free energy: 2.3352554251
* Residual: 9.856e-6
* Stress norm: 8.65e-6
* Model type: SimpleFieldModel
* Model compressiblity: Incompressible
-----
PolymerSystem (AB) contains 1 components:

Component AB with ϕ=1.0 and α=1.0 contains BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.36 of specie A
  * PolymerBlock B with f=0.64 of specie B


with Flory-Huggins interaction parameters betwen species:
    (A, B) => 15.0
-----
Simulation Cell: BravaisLattice
  * Centering: p
  * Space group: #1 (p1)
  * Crystal system: Line
  * Unit cell: [3.5689541985342235] [0]
  * Free lattice parameters: [a]
-----
* Spatial resolution: (24,)
* Contour steps: 0.01
* MDE solvers: Dict[Dict{Any, Any}((2 => 3) => OSF, (2 => 1) => OSF, (1 => 2) => OSF, (3 => 2) => OSF)]
-----
SCFT updater: PicardMann iteration with α=0.2.
Final state: n=220, #fevals=220, F=2.3352554251, residual=9.86e-6.
-----
* First density field: 0.36 [0.0617, 0.8146]
* First auxiliary field: 2.115 [-4.945, 8.1]
)

JP.F(scft_sgc) - JP.F(scft_vc)

5.024706872802653e-9

Scattering.edges(JP.unitcell(scft_sgc))[1] - Scattering.edges(JP.unitcell(scft_vc))[1]

4.5036503560247354e-5

Gradient-free cell optimization

scft_co = JP.clone(scft_AB_lam);

JP.cell_solve!(OptimGradientFreeCellOpt(), scft_co)

[ Info: UnitCell: 3.0
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 721.3580765177
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: NaN
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.0]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: final loss: 9.961214716192046e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3619039968
[ Info: final loss: 9.96e-6
[ Info: final stress norm: 0.09399202959385276
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.0]
└   * Angles: [0]
[ Info: iterations: 1
[ Info: time per iteration: 261 microseconds, 42 nanoseconds
[ Info: Run time: 00:00:00.000261042
[ Info: =======================================
[ Info: UnitCell: 3.4721359549995796
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: Inf
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: 465.0
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.4721359549995796]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: number: 100                                                                     
[ Info: F: 2.3359890624                                                                 
[ Info: resediual norm: 0.000764                                                        
[ Info: number: 200                                                                     
[ Info: F: 2.3359923615                                                                 
[ Info: resediual norm: 3.91e-5                                                         
[ Info: final loss: 9.933555568331087e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3359924412
[ Info: final loss: 9.93e-6
[ Info: final stress norm: 0.015761693116758988
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.4721359549995796]
└   * Angles: [0]
[ Info: iterations: 256
[ Info: time per iteration: 50 microseconds, 426 nanoseconds
[ Info: Run time: 00:00:00.012909083
[ Info: =======================================
[ Info: UnitCell: 3.76393202250021
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: -5.9720052156
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: NaN
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.76393202250021]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: number: 100                                                                     
[ Info: F: 2.3382653369                                                                 
[ Info: resediual norm: 0.000323                                                        
[ Info: number: 200                                                                     
[ Info: F: 2.3382672387                                                                 
[ Info: resediual norm: 1.68e-5                                                         
[ Info: final loss: 9.891976209708275e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3382672595
[ Info: final loss: 9.89e-6
[ Info: final stress norm: 0.030121219411065932
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.76393202250021]
└   * Angles: [0]
[ Info: iterations: 222
[ Info: time per iteration: 120 microseconds, 608 nanoseconds
[ Info: Run time: 00:00:00.026775042
[ Info: =======================================
[ Info: UnitCell: 3.570524482756406
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 720.2879273302
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: NaN
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.570524482756406]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: number: 100                                                                     
[ Info: F: 2.3352571357                                                                 
[ Info: resediual norm: 0.000225                                                        
[ Info: number: 200                                                                     
[ Info: F: 2.3352560159                                                                 
[ Info: resediual norm: 1.46e-5                                                         
[ Info: final loss: 9.904963834017622e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3352560037







[ Info: final loss: 9.9e-6
[ Info: final stress norm: 0.00026067525310166553
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.570524482756406]
└   * Angles: [0]
[ Info: iterations: 216
[ Info: time per iteration: 50 microseconds, 141 nanoseconds
[ Info: Run time: 00:00:00.010830542
[ Info: =======================================
[ Info: UnitCell: 3.5686982765044233
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 2.3352560037
[ Info: initial residual norm: 9.9e-6
[ Info: initial stress norm: 0.000261
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.5686982765044233]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: final loss: 9.483360800643936e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3352555406
[ Info: final loss: 9.48e-6
[ Info: final stress norm: 3.233796675738711e-5
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.5686982765044233]
└   * Angles: [0]
[ Info: iterations: 72
[ Info: time per iteration: 48 microseconds, 995 nanoseconds
[ Info: Run time: 00:00:00.003527709
[ Info: =======================================
[ Info: UnitCell: 3.5680292782592455
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 2.3352555406
[ Info: initial residual norm: 9.48e-6
[ Info: initial stress norm: 3.24e-5
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.5680292782592455]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: final loss: 9.661794706961047e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3352555876
[ Info: final loss: 9.66e-6
[ Info: final stress norm: 0.00013981526123044234
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.5680292782592455]
└   * Angles: [0]
[ Info: iterations: 56
[ Info: time per iteration: 49 microseconds, 924 nanoseconds
[ Info: Run time: 00:00:00.00279575
[ Info: =======================================
[ Info: UnitCell: 3.569055146332074
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 2.3352555876
[ Info: initial residual norm: 9.66e-6
[ Info: initial stress norm: 0.00014
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.569055146332074]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: final loss: 9.401436432909804e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3352551342
[ Info: final loss: 9.4e-6
[ Info: final stress norm: 2.4917020472009676e-5
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.569055146332074]
└   * Angles: [0]
[ Info: iterations: 59
[ Info: time per iteration: 50 microseconds, 322 nanoseconds
[ Info: Run time: 00:00:00.002969
[ Info: =======================================
[ Info: UnitCell: 3.5696163829052607
[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (Residual)
[ Info: initial F: 2.3352551342
[ Info: initial residual norm: 9.4e-6
[ Info: initial stress norm: 2.49e-5
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.5696163829052607]
└   * Angles: [0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: final loss: 9.920604750135453e-6
[ Info: Stop triggered by EarlyStopping.Threshold(1.0e-5) stopping criterion. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3352553413
[ Info: final loss: 9.92e-6
[ Info: final stress norm: 0.00011500648824438853
┌ Info: final unit cell: UnitCell
│   * Crystal system: Line
│   * Edges: [3.5696163829052607]
└   * Angles: [0]
[ Info: iterations: 50
[ Info: time per iteration: 49 microseconds, 62 nanoseconds
[ Info: Run time: 00:00:00.002453125
[ Info: =======================================

(Results of Optimization Algorithm
 * Algorithm: Brent's Method
 * Search Interval: [2.236068, 4.236068]
 * Minimizer: 3.569055e+00
 * Minimum: 2.335255e+00
 * Iterations: 7
 * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.0e-04*|x|+2.2e-16): true
 * Objective Function Calls: 8, Noncyclic Chain SCFT model:
* Free energy: 2.3352553413
* Residual: 9.921e-6
* Stress norm: 0.000115
* Model type: SimpleFieldModel
* Model compressiblity: Incompressible
-----
PolymerSystem (AB) contains 1 components:

Component AB with ϕ=1.0 and α=1.0 contains BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.36 of specie A
  * PolymerBlock B with f=0.64 of specie B


with Flory-Huggins interaction parameters betwen species:
    (A, B) => 15.0
-----
Simulation Cell: BravaisLattice
  * Centering: p
  * Space group: #1 (p1)
  * Crystal system: Line
  * Unit cell: [3.5696163829052607] [0]
  * Free lattice parameters: [a]
-----
* Spatial resolution: (24,)
* Contour steps: 0.01
* MDE solvers: Dict[Dict{Any, Any}((2 => 3) => OSF, (2 => 1) => OSF, (1 => 2) => OSF, (3 => 2) => OSF)]
-----
SCFT updater: PicardMann iteration with α=0.2.
Final state: n=50, #fevals=50, F=2.3352553413, residual=9.92e-6.
-----
* First density field: 0.36 [0.0617, 0.8147]
* First auxiliary field: 2.115 [-4.946, 8.1]
)

JP.F(scft_co) - JP.F(scft_vc)

-7.873114560297267e-8

Scattering.edges(JP.unitcell(scft_co))[1] - Scattering.edges(JP.unitcell(scft_vc))[1]

0.0007072208745975139

Exercises

1. Compute the free energies of the LAM phase for the same system with χN=20.0 and fA=[0.36, 0.38, 0.40, 0.42].

HEX

To compute the hexagonal cylinder phase, we construct a simulation cell with simple hexagonal unit cell an space group 17. The space group is used by symmetrizing functionality of Polyorder if enabled.

scft_AB_hex = let
    ds = 0.01
    uc = UnitCell(Hexagonal2D(), 3.6)  # A simple cell with side length 3.6 Rg
    lat = BravaisLattice(uc, 17)  # space group No.17
    system = AB_system(χN=15.0, fA=0.36)
    NoncyclicChainSCFT(system, lat, ds)
end;

let
    scftconfig = SCFTConfig(symmetrize=false, max_iter=500, tolmode=:F, tol=1e-8)
    config = Polyorder.Config(scft=scftconfig)
    JP.solve!(scft_AB_hex, config)
end

[ Info: ******* SCFT Simulation Start *******
[ Info: Algorithm: PicardMann iteration with α=0.2.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-8 (F)
[ Info: initial F: NaN
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: NaN
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Hexagonal2D
│   * Edges: [3.6, 3.6]
└   * Angles: [2π/3, 0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: number: 100                                                                     
[ Info: F: 2.3728691315                                                                 
[ Info: resediual norm: 1.31                                                            
[ Info: number: 200                                                                     
[ Info: F: 2.3276793157                                                                 
[ Info: resediual norm: 0.0207                                                          

[ Info: final loss: 0.0009996972397399162
[ Info: Stop triggered by a `ThresholdObjFun` control. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.3276901087
[ Info: final loss: 0.001
[ Info: final stress norm: 0.029823622728354825
┌ Info: final unit cell: UnitCell
│   * Crystal system: Hexagonal2D
│   * Edges: [3.6, 3.6]
└   * Angles: [2π/3, 0]
[ Info: iterations: 294
[ Info: time per iteration: 7 milliseconds, 310 microseconds, 967 nanoseconds
[ Info: Run time: 00:00:02.149424417
[ Info: =======================================

Polyorder.Successful()

Plot the density field of any component, and verify that the hexagonal cylinder phase is actually being computed. If not, you have to go back to the above procedure, change the cell size or disable symmetrizer or initialize by a perpared hexagonal cylinder phase.

plot_density(scft_AB_hex)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/8h0bl/src/scenes.jl:227

Optimize the simulation cell:

let
    scftconfig = SCFTConfig(symmetrize=false, max_iter=500, tolmode=:F, tol=1e-8)
    config = Polyorder.Config(scft=scftconfig)
    scft = JP.clone(scft_AB_hex)
    JP.cell_solve!(scft, VariableCell(BB(1.0), SD(0.2)), config)
end

[ Info: ###### Cell Optimization Start ######
┌ Info: Algorithm: Variable cell method
│ * cell updater: Barzilai-Borwein (BB2) method with max step size =1.0.
│ * fields updater: PicardMann iteration with α=0.2.
└ Run fields updater 1 times per cell iteration.
[ Info: Cell optimization starts ...
[ Info: ******* SCFT Simulation Start *******
┌ Info: Algorithm: Variable cell method
│ * cell updater: Barzilai-Borwein (BB2) method with max step size =1.0.
│ * fields updater: PicardMann iteration with α=0.2.
└ Run fields updater 1 times per cell iteration.
[ Info: MDE solvers: ["OSF"]
[ Info: tolerance: 1.0e-5 (stress norm)
[ Info: tolerance: 1.0e-8 (F)
[ Info: initial F: Inf
[ Info: initial residual norm: 0.0
[ Info: initial stress norm: NaN
┌ Info: initial unit cell: UnitCell
│   * Crystal system: Hexagonal2D
│   * Edges: [3.6, 3.6]
└   * Angles: [2π/3, 0]
[ Info: 
[ Info: Simulation starts > > > > > >
[ Info: number: 100                                                                     

[ Info: F: 2.3150912757                                                                 
[ Info: resediual norm: 0.000632                                                        
[ Info: stress norm: 6.54e-8                                                            
[ Info: cell: [4.0307]                                                                  
[ Info: final loss: 0.000610034134998827
[ Info: Stop triggered by a `ThresholdStress` control. 
[ Info: > > > > > > Simulation finished.
[ Info: 
[ Info: ------ SCFT Simulation Summary ------
[ Info: convergence: Polyorder.Successful()
[ Info: final F: 2.315091267
[ Info: final loss: 0.00061
[ Info: final stress norm: 1.86e-7
┌ Info: final unit cell: UnitCell
│   * Crystal system: Hexagonal2D
│   * Edges: [4.0306960068361795, 4.0306960068361795]
└   * Angles: [2π/3, 0]
[ Info: iterations: 101
[ Info: time per iteration: 20 milliseconds, 242 microseconds, 231 nanoseconds
[ Info: Run time: 00:00:02.044465334
[ Info: =======================================
[ Info: Cell optimization finished.
[ Info: ------ Cell Optimization Summary ------
[ Info: Convergence: Polyorder.Successful()
┌ Info: Stress-free cell: UnitCell
│   * Crystal system: Hexagonal2D
│   * Edges: [4.0306960068361795, 4.0306960068361795]
└   * Angles: [2π/3, 0]

[ Info: Final F: 2.315091266973812
[ Info: Final loss: 0.000610034134998827
[ Info: Final stress: 1.859505600132121e-7
[ Info: Total solve! calls: 101
[ Info: Total SCFT iterations: 101
[ Info: Total time: 00:00:03

(Polyorder.Successful(), Noncyclic Chain SCFT model:
* Free energy: 2.315091267
* Residual: 0.00061
* Stress norm: 1.86e-7
* Model type: SimpleFieldModel
* Model compressiblity: Incompressible
-----
PolymerSystem (AB) contains 1 components:

Component AB with ϕ=1.0 and α=1.0 contains BlockCopolymer AB with 2 blocks:
  * PolymerBlock A with f=0.36 of specie A
  * PolymerBlock B with f=0.64 of specie B


with Flory-Huggins interaction parameters betwen species:
    (A, B) => 15.0
-----
Simulation Cell: BravaisLattice
  * Centering: p
  * Space group: #17 (p6mm)
  * Crystal system: Hexagonal2D
  * Unit cell: [4.0306960068361795, 4.0306960068361795] [2π/3, 0]
  * Free lattice parameters: [a]
-----
* Spatial resolution: (27, 27)
* Contour steps: 0.01
* MDE solvers: Dict[Dict{Any, Any}((2 => 3) => OSF, (2 => 1) => OSF, (1 => 2) => OSF, (3 => 2) => OSF)]
-----
SCFT updater: PicardMann iteration with α=0.2.
-----
* First density field: 0.36 [0.0829, 0.9316]
* First auxiliary field: 2.1 [-6.024, 7.431]
)

Exercises

1. Compute the free energies of the HEX phase for the same system with χN=20.0 and fA=[0.36, 0.38, 0.40, 0.42].