Parameter Ramps (Prev) Developer Information (Next)

A perturbation is an optional feature that allows a user to define a modification of the form of the standard field theoretic Hamiltonian. Thus far, it has only been used as a framework for implementing the Einstein crystal method of computing free energies.

Perturbation Concept

Any perturbation of the Hamiltonian must be implemented by a subclass of a base class named Perturbation. A Perturbation defines a functional of the fields, denoted here by \( H_{1} \), that is added to the standard Hamiltonian \( H_{0} \) to obtain a total Hamiltonian \( H_{0} + H_{1} \). The perturbation Hamiltonian must depend upon a single real parameter, denoted by \( \lambda \), giving a modified Hamiltonian of the form

\[ H(\lambda) = H_{0} + H_{1}(\lambda) \]

A Perturbation class must provide methods to compute \( H_{1} \), and to compute its functional derivatives with respect to the fields (which are needed to compute forces in Brownian dynamics simulation) and the the partial derivative of \( H_{1} \) with respect to the parameter \( \lambda \) at fixed field values (which is needed for thermodynamic integration).

A perturbation is an optional element of a simulation that may be added by including a Perturbation subblock of the BdSimulator or McSimulator block of a parameter file. If present, the Perturbation block must appear immediately after the Compressor block (if any). Full descriptions of formats of for the BdSimulator and McSimulator parameter file blocks are given here and here, respectively.

The PerturbationDerivative analyzer class is designed to enable thermodynamic integration with respect to the \( \lambda \) parameter associated with a perturbation. This analyzer is designed to work with any subclass of Perturbation, and relies on the fact that any such Perturbation is required to define a functions to compute the value of the partial derivative \( \partial H_{1} / \partial \lambda \) . The PerturbationDerivative class provides tools to periodically evaluate this partial derivative, compute its time average, and write values or block averages of this quantity to a file.

Continuous thermodynamic integration with respect to the parameter lambda asssociated with a perturbation may be performed by using the parameter type lambda_pert in the parameter block for the ramp to perform. A ramp that uses this parameter type performs a linear variation of the lambda parameter over the course of a simulation, as discussed here . The implementation of this feature of the LinearRamp class only uses the interface defined by the Perturbation base class, and so is designed to work with any type of perturbation (i.e., any subclass of this base class).

Einstein Crystal Perturbation

EinsteinCrystalPerturbation : The only type of Perturbation that is provided with PSCF v1.2 is a class named EinsteinCrystalPerturbation that enables calculation of free energies by the Einstein crystal method. The resulting perturbation is a functional of the form

\[ H_{1}(\lambda) = (1 - \lambda)( H_{{\rm EC}} - H_{0} ) \]

in which \( H_{0} \) is the standard unperturbed Hamiltonian functional, \( H_{{\rm EC}} \) is the Hamiltonian for an Einstein crystal model (described below), and \( \lambda \) is a parameter that must be in the range \( 0 \leq f \leq 1 \).

In the common case of a standard AB system with \( \chi > 0 \), the EinsteinCrystalPerturbation class defines a harmonic functional

\[ H_{{\rm EC}} = \frac{1}{v\epsilon} \int \! d {\bf r} \; [ W_{-}({\bf r}) - W^{({\rm r})}_{-}({\bf r}) ]^2 \]

in which \( W^{({\rm r})}_{-}({\bf r}) \) is a reference exchange field provided by the user and \( \epsilon \) is a scalar parameter that controls the magnitude of fluctuations. The \( \epsilon \) parameter for a standard AB system is set equal to \( \chi \) by default, and this default setting is used in the example given below. Generalization of the Einstein crystal Hamiltonian to cases with more than two monomer types, and the full parameter file format for the class are discussed in manual page for the EinsteinCrystalPerturbation class.

When a perturbation of this form is added to the unperturbed Hamiltonian \( H_{0} \), it yields a total Hamiltonian

\[ H(\lambda) = \lambda H_{0} + (1-\lambda) H_{{\rm EC}} \]

By construction, \( H(\lambda) \) is thus equal to the Einstein crystal Hamiltonian \( H_{{\rm EC}} \) when \( \lambda = 0 \) and equal to the unperturbed Hamiltonian \( H_{0} \) when \( \lambda = 1 \). The derivative of \( H(\lambda) \) with respect to \( \lambda \), which is used in thermodynamic integration, is given by the difference

\[ \frac{\partial H_{1}(\lambda)}{\partial \lambda} = H_{0} - H_{{\rm EC}} \]

where the partial derivative is evaluated at fixed field values.

The parameter file block for an EinsteinCrystalPerturbation requires a user to enter a value for the perturbation parameter \( \lambda \), and the name of a file that contains the reference field configuration in r-grid format, as described here. The reference field configuration is usually taken to be a SCFT solution for the phase of interest. When applied to a disordered phase, this reference field normally contains spatially homogeneous fields.

Thermodynamic Integration : The EinsteinCrystalPerturbation and PerturbationDerivative classes must be used together to compute free energies by the Einstein crystal method. When applied to a standard AB system, with chi(0,0) = chi(1,1) = 0.0 and \( \chi = \) chi(0,1) = chi(1,0), with the default choice of \( \epsilon = \chi \), the free energy \( F \) of the interacting system with \( \lambda = 1 \) can be expressed as a sum

\[ F = H_{0}^{({\rm r})} + \int\limits_{0}^{1}\! d \lambda \; \left \langle \frac{\partial H(\lambda)}{\partial \lambda} \right \rangle \quad. \]

in which \( H_{0}^{({\rm r})} \) denotes the value of the standard unperturbed Hamiltonian in the reference state used by the Einstein crystal Hamiltonian. In this expression, the value of the integrand at each value of \( \lambda \in [0, 1]\) is an equilibrium average for a system governed by a perturbed Hamiltonian \( H(\lambda) \rangle \).

The required integral with respect to \( \lambda \) can be evaluated by either of two methods:

Multiple equilibrium simulations : Perform separate equilibrium simulations with an EinsteinCrystalPerturbation at multiple values of \( \lambda \in [0,1] \) and use a standard numerical integration scheme to approximate the integral. For example, one can obtain the integral over \( [0, 1] \) or any subinterval thereof by performing simulations at an odd number of evenly spaced values of \( \lambda \) (including the endpoints) that divide the interval into an even number of subintervals, and then use Simpson's rule for integration.
Continuous integration : Perform a very slow linear ramp of \( \lambda \) from 0 to 1 by using the LinearRamp class with a parameter of type lambda_pert. The integral can then be approximated by the time average of \( \partial H(\lambda)/\partial \lambda \) over the course of this ramp.

Ramps can also be applied to subdomains of the domain \( [0 1] \) in order to use different ramp rates over different ranges. When computing the free energy of an ordered phase using an inhomogeneous reference field, special attention to accuracy is required for a range of values near the final value \( \lambda = 1 \). In either method, the required values of \( \partial H(\lambda)/\partial \lambda \) can be sampled by the PerturbationDerivative analyzer.

Example: Continuous Einstein Crystal Integration

Below, we show a complete parameter file for a calculation that uses continuous Einstein crystal integration to compute a free energy difference.

System{
  Mixture{
    nMonomer  2
    monomers[
             1.0  
             1.0 
    ]
    nPolymer  1
    Polymer{
      type    linear
      nBlock  2
      blocks[
              0   0.5
              1   0.5
      ]
      phi     1.0
    }
    vMonomer  0.01
    ds   0.02
  }
  Interaction{
    chi(
        1   0   14.0
    )
  }
  Domain{
    mesh         40    40    40
    lattice      cubic   
  }
  BdSimulator{
    seed           486893701
    BdStep{
      mobility  5.0E-2
    }
    Compressor{
      epsilon      1.0e-4
      maxHist      40
    }
    EinsteinCrystalPerturbation{
      lambda           0.0
      fieldFileName    in/ecReference.rf
    }
    Ramp{
      nParameters  1 
      parameters[
         lambda_pert    1.0
      ]
    }
    AnalyzerManager{
 
      StepLogger{
        interval         50
      }
 
      PerturbationDerivative{
        interval          5
        outputFileName    out/derivative
      }
 
    }
  }
}

Discussion : The main new elements in this example, compared to simpler previous examples, are the presence of an EinsteinCrystalPerturbation block and a PerturbationDerivative analyzer, as well as a linear Ramp that is used to perform continuous integration. Here, we have used generic block labels BdStep, Compressor, and Ramp to select default algorithms for these blocks, but have used the specific class name EinsteinCrystalPerturbation for the Perturbation block. Because EinsteinCrystalPerturbation is the default perturbation algorithm, we could instead have used a generic label Perturbation for this block without changing the result.

In the EinsteinCrystalPerturbation block, we set an initial value lambda = 0.0 for the perturbation parameter, corresponding to a pure Einstein crystal Hamiltonian. The epsilon parameter has been set equal to the value 14.0 that is assigned to chi(0,1) in the interaction block, using the convention alpha = chi(0,1) recommended above. The reference field is read from an r-grid field file named ecReference.rf in subdirectory in/ of the current working directory (i.e., the directory from which the program was invoked).

The Ramp block uses the parameter identifier lambda_pert to indicate that the lambda parameter of the associated Perturbation should be varied, and gives a range 1.0 to creat a ramp from the initial value of 0 to a final value of 1.0.

In the PerturbationDerivative block of the AnalyzerManager, this analyzer is programmed to sample the partial derivative \( \partial H/\partial \lambda \) every 5 simulation steps, and to output every sampled value to a file. The file name out/derivative will be used as a base file name for the file containing these sampled values and for another file that contains information about the time average.

Perturbation Concept

Einstein Crystal Perturbation

Example: Continuous Einstein Crystal Integration

See Also: