PSCF v1.1
Pspc::AmIterator

The AmIterator algorithm used by the pscf_pc programs solves the SCFT equations for a periodic system using the Anderson Mixing algorithm.

The variant of this algorithm implemented here is almost identical to that described by Arora et al. [Arora, Morse, Bates, and Dorfman, J. Chem. Phys. vol. 146, 244902 (2017)], which described the implementation of this algorithm in the earlier Fortran version of PSCF. The version of this algorithm used here adjusts the fields by adjusting the compoments of the field in a basis of symmetry adapated basis functions, thus preserving the space-group symmetry of the fields throughout the iteration process. The algorithm can be used to either solve the SCFT equations for a rigid unit cell or to solve the SCFT equations and also optimize some or all of the unit cell parameters of a flexible unit cell so as to minimize the free energy density, depending on the value of "isFlexible" and "flexibleParams", as described below.

Parameter File

A typical example of the parameter file format for this iterator is shown below:

AmIterator{
epsilon 1e-8
maxItr 100
maxHist 50
isFlexible 1
}

The format of this block is described more formally below:

AmIterator{
epsilon real
maxItr* int (200 by default)
maxHist* int (50 by default)
verbose* int (0-2, 0 by default)
outputTime* bool (false by default)
errorType* string ("norm", "rms", "max", or "relNorm", "relNorm" by default)
isFlexible* bool (0 or 1, 1/true by default)
flexibleParams* Array [ bool ] (nParameters elements)
scaleStress* real (10.0 by default)
}

Here, as elsewhere, labels followed by an asterisk (*) represent optional parameters. The meaning of the various parameters are described below:

Label Description
epsilon Desired tolerance for convergence - iteration stops if the magnitude of the error drops below epsilon.
maxItr* Maximum number of iterations that will be attempted (200 by default)
maxHist* Maximum number of previous trial solultions that will be retained in memory for used by the AM algorithm (50 by default)
verbose* Integer level 0, 1, 2 for verbosity of log output during iteration. Optional, and 0 (concise) by default.
showTime* Set to 1 (true) to report wall clock time components in log file. Optional, and 0 (false) by default.
errorType* Identifer for the type of variable used to define scalar error . The only allowed values are "norm", "rms" "max", and "relNorm", as discussed below. Optional, and equal to "relNorm" by default.
isFlexible* Set isFlexible true to enable or false to disable iterative optimization of unit cell parameters so as to minimize the free energy density. Optional and true by default.
flexibleParams* Use flexibleParams to declare only some lattice parameters as flexible. This optional entry is allowed only if isFlexible is true. This array must contain one entry for each lattice parameter, with a value 1 (true) if the parameter is flexible or 0 (false) if the parameter value is fixed. If isFlexible is true but flexibleParams is absent, then all all lattice parameters will be flexible.
scaleStress* Constant factor by which stress components are multiplied in the definition of the residual attempted if isFlexible is true (optional).

The iterative loop exits if the number of iterations has reached maxItr or if the magnitude of the scalar error drops below epsilon.

errorType : Several different definitions may be used for the scalar error, depending on the value of the identifier errorType.
All of them ar based on computation of different types of norm of the residual vector, defined as discussed below.

  • If errorType == norm, then the scalar error is take to be the L2 norm of the residual vector (square-root of sum of squares of the elements) of the residual vector as defined below.
  • If errorType == rms, then the scalar error is take to be the root-mean-square magnitude of elements of the residual vector (the L2 norm divided by the square root of the number of elements).
  • If errorType == max, then the scalar error is take to be the maximum of the absolute magnitude of the elements of the residual vector (the L infinity norm of the residual vector).
  • If errorType == relNorm, then the scalar error is take to be the ratio of the L2 norm of the residual vector to the L2 norm of the w field, as in Stasiak and Matsen, Eur. Phys. J. E 34 , 110 (2011).

scaleStress : If isFlexible is true, the choice of a value for the parameter scaleStress determines how strongly the definition of the scalar error weights errors that arise from nonzero derivatives of the free energy with respect to the unit cell parameters, relative to errors arising from errors in the SCFT equations for the w fields in a fixed unit cell. This variable is irrelevant if isFlexible is false, and should normally be ommitted in this case.

Residual Definition

The vector of residuals used in this algorithm is described by Eqs. (10-13) of the article by Arora et al. (J. Chem. Phys. 2017). To describe the residuals here, we use zero based indexing for monomer types and basis functions, us use \( N_{b} \) to denote the number of basis functions (referred to in the source code as Nbasis), and use \( N_{m} \) to denote the number of monomer types (referred to elsewhere as nMonomer).

The vector of residuals for a calculation for a rigid unit cell (i.e., one with isFlexible set false) has \( N_{b}N_{m} \) (or nBasis*nMonomer) scalar residuals.
Let \( w_{ia} \) and \( \phi_{ai} \) denote the coefficients of basis function a in the symmetry-adapted Fourier expansions of the chemical potential field \( w_{i} \) and volume fraction \( \phi_{i} \) for monomer type i, respectively. The residual component \( R_{ai}\) associated with basis function a and monomer type i is given, as in Eq. (10) of Arora et al. (JCP 2017), as

\[ R_{ai} = \sum_{j=0}^{N_{m}-1} \left ( \chi_{ij} \phi_{aj} - P_{ij}w_{aj} \right ) \quad. \]

Here, \( P_{ij} \) is an element of an idempotent \( N_{m} \times N_{m} \) matrix \( P \) that can be expressed in matrix notation as

\[ P = I - \frac{\epsilon \epsilon^{T} \chi^{-1} } { \epsilon^{T} \chi^{-1} \epsilon } \quad. \]

where \( \chi^{-1} \) is the matrix inverse of the \(\chi\) matrix, and \(\epsilon\) is a \( N_{m} \)-component column vector with elements given by \( \epsilon_{i} = 1 \) for all \( i = 0, \ldots, N_{m}-1 \).

The vector of residuals for a calculation with a flexible unit cell (with isFleixble set true) has an additional stress residual associated with each flexible unit cell parameter. The residual \( R_{m} \) associated with flexible unit cell parameter \( \theta_{m} \) is given by

\[ R_{m} = -S\frac{\partial f}{\partial \theta_{m}} \]

where \( f \) is the free energy per monomer (i.e., the free energy density times the monomer reference volume) in units in which \( kT = 1 \), while \(S\) is an arbitrary prefactor whose value is given by the parameter scaleStress. (The factor \(S\) is set to 10.0 by default).

Closed Systems (Canonical Ensemble)

A slight modification of the residual definition is required in the case case when a value of phi rather than mu is specified for every molecular species, giving a closed system, or a canonical statistical ensemble. In this case, the solution of the SCF problem is not unique, becase the values of the residuals defined above can be shown to be invariant under a shift in all of the chemical potential fields at all grid point points by any spatially homgeneous constant. To obtain a unique solution, in this case the spatial average of the Laplace multiplier field \( \xi({\bf r}) \) is implicitly set to zero by convention.