mutable struct Wizard

holds lots of HNC data.

Important fields for input

• Also refer to immutable arguments in Wizard()

Immutable fields set by Wizard()

• ncomp: number of chemical components
• nfnc : number of functions = ncomp (ncomp + 1) / 2
• r : real space grid
• k : reciprocal space grid (computed)
• deltak : reciprocal space grid spacing (computed)

Mutable fields, which may be modified before invoking a solver.

• rho : density array
• z: valence array
• arep :repulsion amplitude array
• dd: hard core diameter array
• diam : hard core diameter array, pairwise
• sigma : used both for the long-range Coulomb smearing length for the soft potentials, and for hard core diameter for RPM models.
• sigmap : +- long-range Coulomb smearing length (URPM)
• kappa : +- long-range Coulomb smoothing parameter (RPM)
• rgroot : linear charge smearing range (Groot)
• lbda : Slater charge smearing range (exact)
• beta : Slater charge smearing range (approx)
• alpha: Picard method, fraction of new solution
• rc : short-range DPD repulsion range
• lb : long-range Coulomb coupling length

Computed and allocated fields, calculated by solvers.

• c : direct correlation functions (dcfs)
• e : indirect correlation functions (icfs)
• ck : transform of dcfs
• ek : transform of icfs
• sk : partial structure factors
• h0 : reference total correlation functions
• hr : current total correlation functions
• hc : current contact values
• ushort : short range potential in real space
• dushort : derivative of ushort
• ulong : long range potential in real space
• dulong : derivative of ulong
• ulongk : long range potential in reciprocal space
• expnegus : exp(-ushort) (includes hard cores)
• u0 : r = 0 of long range potential (diagonal)
• muex: chemical potential array
• tp : mean field pressure contribution
• tu : mean field energy contribution
• tl : mean field long range potential

Important fields for output

• ERROR : difference between current and previous solutions
• return_code : return (error) code, any of return_code
• closure_type : any of closure_type
• model_type : which potential was (last) chosen
• nsteps : actual number of steps taken

Working fields for FFT

• planFFT
• fftw

Some fields are calculated by following methods

• make_pair_functions!
• make_structure_factors!
• make_thermodynamics!

Wizard(; kwargs...)

Construct Wizard object.

Immutable arguments, each used as the struct's field as it is

• ncomp::Int64=1 >= 1: number of chemical components
• ng::Int64=2^8 > 0: grid size
• nps::Int64=6 >= 2: number of previous states used in Ng method
• npic::Int64=6 >= 2: number of Picard steps
• deltar::Float64=1e-2 > 0 : real space grid spacing

Mutable arguments, used to control solver

• tol::Float64=1e-12 > 0 : Error tolerance for claiming convergence
• maxsteps::Int64=100 : max number of steps to take for convergence
• verbose::Bool=false : Print solver diagnostics
• silent::Bool=false : Prevent printing warning/error messages
• suppress_msgs::Bool=false : In case of severe numerical instability
• auto_fns:Bool=true : Calculate stuff at end
• cold_start : Force a cold start
• start_type : how to initialise in a cold start

dpd_potential!(w::Wizard, charge_type=DPD_GAUSSIAN_CHARGES)

Build the potential arrays.

Arguments

• charge_type : Use the defined integer constants for these.
  • DPD_GAUSSIAN_CHARGES=1, Gaussian (default)
  • DPD_BESSEL_CHARGES=2, Bessel
  • DPD_LINEAR_CHARGES=3, Groot (linear)
  • DPD_SLATER_APPROX_CHARGES=4, Slater (approx)
  • DPD_SLATER_EXACT_CHARGES=5, Slater (exact)

Referred fields in Wizard object

A factor beta = 1/kT is implicit in the following definitions.

Fields referred by all charge_type

• arep[:,:]
• z

rpm_potential!(w::Wizard, use_ushort=false)

Build the potential arrays for the softened RPM (charged hard spheres).

Using kappa < 0 implies the pure RPM case (kappa -> infinity).

The case ncomp = 3 corresponds to the RPM in the presence of a neutral hard sphere solvent.

Arguments

• use_ushort::Bool=false: whether ushort is used or not.

Referred fields for input

• diam
• sigma
• lb
• kappa

urpm_potential!(w::Wizard, use_ushort=false)

Build the potential arrays for the softened URPM (Gaussian charges), This expects ncomp = 2, and will set z[1] = 1, z[2] = -1.

Arguments

• use_ushort::Bool=false: whether ushort is used or not.

Referred fields for input

• sigma
• sigmap
• lb

hs_potential!(w::Wizard)

Build the potential arrays for hard spheres with diameter sigma. This works for arbitrary numbers of components.

Referred fields for input

• diam
• sigma

NO_MODEL_TYPE - constant for charge_type

DPD_GAUSSIAN_CHARGES - constant for charge_type

DPD_BESSEL_CHARGES - constant for charge_type

DPD_LINEAR_CHARGES - constant for charge_type

DPD_SLATERAPPROX\CHARGES - constant for charge_type

DPD_SLATEREXACT\CHARGES - constant for charge_type

URPM_WITHOUT_USHORT - constant for charge_type

URPM_WITH_USHORT - constant for charge_type

RPM_WITHOUT_USHORT - constant for charge_type

RPM_WITH_USHORT - constant for charge_type

HARD_SPHERES - constant for charge_type

oz_solve!(w::Wizard)

Solve the Ornstein-Zernicke equation to determine e = h - c, given c. We re-partition the long range part of the potential so that the routine actually calculates c' = c + Ulong and e' = e - Ulong. This is because the Fourier transform of Ulong can be computed in closed form. Note h = e + c = e' + c'.

oz_solve2!(w::Wizard)

Solve an alternate version of the Ornstein-Zernicke equation to determine c and e from the reference h0. In practice as always we actually calculate c' = c + Ulong and e' = e - Ulong. Note h = e + c = e' + c'. The result is saved to position 1 in the history trajectory. Note that the offset ulongk in.

hnc_ng!(w::Wizard)

This method implements the Ng method K-C Ng, J. Chem. Phys. v61, 2680 (1974), doi:10.1063/1.1682399 as an accelerated solver for the HNC closure. As we iterate we move forward in the history trajectory, cyclically.

hnc_picard!(w::Wizard)

This method implements the HNC closure expressed as c = exp( -beta v + e) - e - 1 where e = h - c is the indirect correlation function, c is the direct correlation function from the Ornstein-Zernicke relation, h = g - 1, and g is the pair distribution function. One can show this is equivalent to g = exp(-v + h - c) in Hansen + McDonald. As above, the routine actually works with c' = c + Ulong and e' = e - Ulong where Ulong is the long-range part of the potential for which the Fourier transform is simple. This means that v in the above expression is the short-range part of the potential only. As we iterate we move forward in the history trajectory, cyclically.

hnc_solve!(w::Wizard)

Basic driver routine for solving HNC: take a number of Picard iterations to pump-prime the Ng method. Stop when error is less than tolerance, or when exceed maximum number of iterations. The flag cold_start indicates whether the direct correlation function should be re-initialised.

The initial guess to the direct correlation function is either

• 0 (zero) (start_type = 1), or
• c = - Ushort (start_type = 2), or
• c = e^(-Ushort)-1 (start_type = 3).

Any of these should do in principle, but the initial convergence may be different. Note from above that c is actually defined c' = c + Ulong, ie with the long-range part of the potential added.

Note that we always start from position 1 in the history trajectory, and the history trajectory is pump-primed by Picard steps before attempting the Ng accelerator. At the end, the final solution is copied back to position 1 in the history trajectory.

msa_ng!(w::Wizard)

This method implements the Ng method K-C Ng, J. Chem. Phys. v61, 2680 (1974), doi: 10.1063/1.1682399 as an accelerated solver for the MSA closure (cf HNC above). As we iterate we move forward in the history trajectory, cyclically.

msa_picard!(w::Wizard)

This method implements the MSA closure expressed as c' = - e' - 1 within the hard core, in similar terms to the HNC closure. Outwith the hard core, c' = - beta v' is left untouched, presuming it is set correctly in the MSA initialisation step. As we iterate we move forward in the history trajectory, cyclically.

msa_solve!(w::Wizard)

Basic driver routine for solving MSA: take a number of Picard iterations to pump-prime the Ng method. Stop when error is less than tolerance, or when exceed maximum number of iterations. The flag cold_start indicates whether the direct correlation function should be re-initialised.

For a cold start, the initial guess to the direct correlation function inside the hard core is c' = -1.

Irrespective of this, we need to make sure c' is correctly initialised outwith the hard core.

Note that we always start from position 1 in the history trajectory, and the history trajectory is pump-primed by Picard steps before attempting the Ng accelerator. At the end, the final solution is copied back to position 1 in the history trajectory.

rpa_solve!(w::Wizard)

Given the HNC machinery, the implementation of the RPA is almost completely trivial and corresponds to one iteration through the Ornstein-Zernike solver given the choice c = - Ushort. We save the result to position 1 in the history trajectory.

NO_ERROR – constant for return_code

CONVERGENCE_ERROR - constant for return_code

AXEQB_ERROR - constant for return_code

DSYSV_ERROR - constant for return_code

DGEEV_ERROR - constant for return_code

MISMATCH_ERROR - constant for return_code

NO_CLOSURE - constant for closure_type

HNC_CLOSURE - constant for closure_type

RPA_CLOSURE - constant for closure_type

MSA_CLOSURE - constant for closure_type

EXP_CLOSURE - constant for closure_type

conv_test!(w::Wizard)

Calculate the difference between the direct correlation functions for the current and previous iteration, used as a convergence test; return answer in variable 'error'.

hnc_kspace_integrals!(w::Wizard)

Calculate the reciprocal space contributions to the HNC free energy.

Note that c = c' - Ulong, and ck is available after a call to OZ solver.

exp_refine!(w::Wizard)

The EXP approximation refines the current RPA/MSA solution by using the current solution (h = c + e) and a reference solution (h0) to send h0 --> (1 + h0) exp(h - h0) - 1. A round trip through the alternate version of the Ornstein-Zernike relation re-calculates the direct and indirect correlation functions. We assume the current solution is in position 1 in the history trajectory. The reference state is reset to h0 = 0 after a call to this function, for safety!

make_pair_functions!(w::Wizard)

Construct the total correlation functions out of the direct correlation functions. Note that the above routines actually works with c' = c + Ulong and e' = e - Ulong where Ulong is the long-range part of the potential, but h = g - 1 = e + c = e' + c'. The pair correlation functions are g = 1 + h - the addition of 1 is left for the user to implement. We assume the c and e functions are those in the position 1 in the history trajectory.

Fields calculated by this method

• hc : current contact values
• hr : current total correlation functions

make_structure_factors!(w::Wizard)

Construct the structure factors out of the transform of the total correlation function. Note that ck and ek are available after a call to the OZ solver.

Fields calculated by this method

• hk : fourier transformed pair functions
• sk : fourier transformed partial structure factors

make_thermodynamics!(w::Wizard)

Calculate various thermodynamics properties by spatial integration (as contrasted to thermodynamic integration). We use the trapezium rule, taking account where necessary the end-point values (see intro), h_0 = 2 h_1 - h_2 at r = 0 (i = 1), and h_n = 0 at r = L (i = ng). Also we have r = i * deltar for i = 1 to ng-1. Note that the above routines actually work with c' = c + Ulong and e' = e - Ulong where Ulong is the long-range part of the potential, so we have h = g - 1 = e + c = e' + c'. See also Vrbka et al, J. Chem. Phys. 131, 154109 (2009), doi:10.1063/1.3248218. We assume the c and e functions are those in the position 1 in the history trajectory.

The mean-field thermodynamic expressions can often be obtained analytically from the potential. In this routine they are calculated from species pair contributions, which are themselves calculated in the potential routines (which do not have access to the densities).

Fields calculated by this method

• cf_mf : Compressibility factor: mean field contribution
• cf_xc : Compressibility factor: contact contribution
• cf_gc : Compressibility factor: correlation contribution
• pex : Excess pressure, virial route
• press : Pressure, virial route
• comp : Compressibility: total
• comp_xc : Compressibility: correlation contribution
• uex : Internal energy density
• uex_mf : Internal energy density: mean field contribution
• uex_xc : Internal energy density: correlation contribution
• aex_rl : HNC excess free energy density, real space (long)
• aex_rs : HNC excess free energy density, real space (short)
• aex_kd : HNC excess free energy density, k-space log(det)
• aex_kl : HNC excess free energy density, k-space tr (long)
• aex_ks : HNC excess free energy density, k-space tr (short)
• aex : HNC total free energy density
• deficit : HNC excess free energy deficit (see docs)

save_reference(w::Wizard)

Save the reference state, assuming the c and e functions are those in the position 1 in the history trajectory. The corresponding c and e functions can be restored by a call to oz_solve2.

write_params(w::Wizard)

Write out parameters for system and potentials.

write_thermodynamics(w::Wizard)

Write out parameters for thermodynamics.

upperTriangularIndices(n::Int)

Traverse every element [i,j] of the upper triangular matrix

upperTriangularIndices(w::Wizard)

Equivalent to upperTriangularIndices(w.ncomp).

upperTriangularIndices(n::Int, nlimit::Int)

Traverse every element [i,j] of the upper triangular matrix. Restrict i and j <= nlimit

upperTriangularIndices(w::Wizard, nlimit::Int)

Equivalent to upperTriangularIndices(w.ncomp, nlimit).

diagonalUpperTriangularIndices(n::Int)

Traverse diagonal element [i,i] of the upper triangular matrix

diagonalUpperTriangularIndices(w::Wizard)

Equivalen to diagonalUpperTriangularIndices(w.ncomp).